laura is hosted by Hepforge, IPPP Durham
Laura++  v1r2
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauAbsFitModel.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2004 - 2013.
3 // Distributed under the Boost Software License, Version 1.0.
4 // (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5 
6 // Authors:
7 // Thomas Latham
8 // John Back
9 // Paul Harrison
10 
15 #include <iostream>
16 #include <vector>
17 using std::cout;
18 using std::cerr;
19 using std::endl;
20 using std::vector;
21 
22 #include "TMessage.h"
23 #include "TMonitor.h"
24 #include "TServerSocket.h"
25 #include "TSocket.h"
26 #include "TSystem.h"
27 #include "TVectorD.h"
28 #include "TVirtualFitter.h"
29 
30 #include "LauAbsFitModel.hh"
31 #include "LauAbsPdf.hh"
32 #include "LauComplex.hh"
33 #include "LauFitter.hh"
34 #include "LauFitDataTree.hh"
35 #include "LauFitNtuple.hh"
36 #include "LauGenNtuple.hh"
37 #include "LauParameter.hh"
38 #include "LauParamFixed.hh"
39 #include "LauPrint.hh"
40 #include "LauSPlot.hh"
41 
42 ClassImp(LauAbsFitModel)
43 
44 
46  useAsymmFitErrors_(kFALSE),
47  compareFitData_(kFALSE),
48  writeLatexTable_(kFALSE),
49  writeSPlotData_(kFALSE),
50  storeDPEff_(kFALSE),
51  randomFit_(kFALSE),
52  emlFit_(kFALSE),
53  poissonSmear_(kFALSE),
54  enableEmbedding_(kFALSE),
55  twoStageFit_(kFALSE),
56  usingDP_(kTRUE),
57  pdfsDependOnDP_(kFALSE),
58  firstExpt_(0),
59  nExpt_(0),
60  evtsPerExpt_(0),
61  iExpt_(0),
62  inputFitData_(0),
63  fitNtuple_(0),
64  genNtuple_(0),
65  sPlotNtuple_(0),
66  fitStatus_(0),
67  NLL_(0),
68  numberOKFits_(0),
69  numberBadFits_(0),
70  nParams_(0),
71  nFreeParams_(0),
72  worstLogLike_(DBL_MAX),
73  withinMINOS_(kFALSE),
74  nullString_(""),
75  doSFit_(kFALSE),
76  sWeightBranchName_(""),
77  sWeightScaleFactor_(1.0),
78  fitToyMCFileName_("fitToyMC.root"),
79  fitToyMCTableName_("fitToyMCTable"),
80  fitToyMCScale_(10),
81  fitToyMCPoissonSmear_(kFALSE),
82  sPlotFileName_(""),
83  sPlotTreeName_(""),
84  sPlotVerbosity_(""),
85  socketMonitor_(0),
86  sMaster_(0),
87  vectorPar_(0),
88  vectorRes_(0),
89  messageFromSlave_(0),
90  messageFromMaster_(0),
91  slaveId_(-1),
92  nSlaves_(0)
93 {
94 }
95 
96 
98 {
99  delete inputFitData_; inputFitData_ = 0;
100  delete fitNtuple_; fitNtuple_ = 0;
101  delete genNtuple_; genNtuple_ = 0;
102  delete sPlotNtuple_; sPlotNtuple_ = 0;
103  delete socketMonitor_; socketMonitor_ = 0;
104  delete sMaster_; sMaster_ = 0;
105  for ( std::vector<TSocket*>::iterator iter = sSlaves_.begin(); iter != sSlaves_.end(); ++iter ) {
106  delete (*iter);
107  }
108  sSlaves_.clear();
109  delete vectorPar_; vectorPar_ = 0;
110  delete vectorRes_; vectorRes_ = 0;
113 }
114 
115 // It's necessary to define an external function that specifies the address of the function
116 // that Minuit needs to minimise. Minuit doesn't know about any classes - therefore
117 // use gMinuit->SetFCN(external_function), gMinuit->SetObjectFit(this).
118 // Here, we use TVirtualFitter* fitter instead of gMinuit, defined below.
119 // Then, within the external function, invoke an object from this class (LauAllModel),
120 // and use the member functions to access the parameters/variables.
121 extern void logLikeFun(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
122 
123 void LauAbsFitModel::run(const TString& applicationCode, const TString& dataFileName, const TString& dataTreeName,
124  const TString& histFileName, const TString& tableFileName)
125 {
126  // Chose whether you want to generate or fit events in the Dalitz plot.
127  // To generate events choose applicationCode = "gen", to fit events choose
128  // applicationCode = "fit".
129 
130  TString runCode(applicationCode);
131  runCode.ToLower();
132 
133  TString histFileNameCopy(histFileName);
134  TString tableFileNameCopy(tableFileName);
135  TString dataFileNameCopy(dataFileName);
136  TString dataTreeNameCopy(dataTreeName);
137 
138  // Initialise the fit par vectors. Each class that inherits from this one
139  // must implement this sensibly for all vectors specified in clearFitParVectors,
140  // i.e. specify parameter names, initial, min, max and fixed values
141  this->initialise();
142 
143  if (dataFileNameCopy == "") {dataFileNameCopy = "data.root";}
144  if (dataTreeNameCopy == "") {dataTreeNameCopy = "genResults";}
145 
146  if (runCode.Contains("gen")) {
147 
148  if (histFileNameCopy == "") {histFileNameCopy = "parInfo.root";}
149  if (tableFileNameCopy == "") {tableFileNameCopy = "genResults";}
150 
151  this->setGenValues();
152  this->generate(dataFileNameCopy, dataTreeNameCopy, histFileNameCopy, tableFileNameCopy);
153 
154  } else if (runCode.Contains("fit")) {
155 
156  if (histFileNameCopy == "") {histFileNameCopy = "parInfo.root";}
157  if (tableFileNameCopy == "") {tableFileNameCopy = "fitResults";}
158 
159  this->fit(dataFileNameCopy, dataTreeNameCopy, histFileNameCopy, tableFileNameCopy);
160 
161  } else if (runCode.Contains("reweight")) {
162 
163  this->weightEvents(dataFileNameCopy, dataTreeNameCopy);
164  }
165 }
166 
167 void LauAbsFitModel::runMaster(const TString& applicationCode, const TString& dataFileName, const TString& dataTreeName,
168  const TString& histFileName, const TString& tableFileName, const UInt_t nSlaves)
169 {
170  // Chose whether you want to generate or fit events in the Dalitz plot.
171  // To generate events choose applicationCode = "gen", to fit events choose
172  // applicationCode = "fit".
173 
174  if ( nSlaves < 1 ) {
175  return this->run( applicationCode, dataFileName, dataTreeName, histFileName, tableFileName );
176  }
177 
178  TString runCode(applicationCode);
179  runCode.ToLower();
180 
181  TString histFileNameCopy(histFileName);
182  TString tableFileNameCopy(tableFileName);
183  TString dataFileNameCopy(dataFileName);
184  TString dataTreeNameCopy(dataTreeName);
185 
186  // Initialise the fit par vectors. Each class that inherits from this one
187  // must implement this sensibly for all vectors specified in clearFitParVectors,
188  // i.e. specify parameter names, initial, min, max and fixed values
189  this->initialise();
190 
191  if (dataFileNameCopy == "") {dataFileNameCopy = "data.root";}
192  if (dataTreeNameCopy == "") {dataTreeNameCopy = "genResults";}
193 
194  if (runCode.Contains("gen")) {
195 
196  if (histFileNameCopy == "") {histFileNameCopy = "parInfo.root";}
197  if (tableFileNameCopy == "") {tableFileNameCopy = "genResults";}
198 
199  this->setGenValues();
200  this->generate(dataFileNameCopy, dataTreeNameCopy, histFileNameCopy, tableFileNameCopy);
201 
202  } else if (runCode.Contains("fit")) {
203 
204  //initialize connections with slaves
205  nSlaves_ = nSlaves;
206  this->initSockets();
207  vectorPar_ = new TVectorD(200);
208  vectorRes_ = new TVectorD(200);
209  messageFromSlave_ = new TMessage(kMESS_OBJECT);
210 
211  if (histFileNameCopy == "") {histFileNameCopy = "parInfo.root";}
212  if (tableFileNameCopy == "") {tableFileNameCopy = "fitResults";}
213 
214  this->fitMaster(dataFileNameCopy, dataTreeNameCopy, histFileNameCopy, tableFileNameCopy);
215 
216  // Close the socket.
217  for ( UInt_t i(0); i<nSlaves_; ++i ) {
218  TString msg = "Finish";
219  sSlaves_[i]->Send(msg);
220  sSlaves_[i]->Close();
221  }
222  }
223 }
224 
225 void LauAbsFitModel::runSlave(const TString& applicationCode, const TString& dataFileName, const TString& dataTreeName,
226  const TString& histFileName, const TString& tableFileName, const TString& addressMaster)
227 {
228  // Chose whether you want to generate or fit events in the Dalitz plot.
229  // To generate events choose applicationCode = "gen", to fit events choose
230  // applicationCode = "fit".
231 
232  TString runCode(applicationCode);
233  runCode.ToLower();
234 
235  if ( runCode != "fit" ) {
236  std::cerr << "Not doing a fit - nothing for me to do!" << std::endl;
237  return;
238  }
239 
240  TString histFileNameCopy(histFileName);
241  TString tableFileNameCopy(tableFileName);
242  TString dataFileNameCopy(dataFileName);
243  TString dataTreeNameCopy(dataTreeName);
244 
245  // Open connection to server
246  char str[64];
247  sMaster_ = new TSocket(addressMaster, 9090);
248  sMaster_->Recv(str, 32);
249  TString msg( str );
250  std::cout << "Socket received message: " << msg << std::endl;
251  msg.Remove(0, 6);
252  TString slaveID = msg(0,msg.Index("/"));
253  TString nSlaves = msg(msg.Index("/")+1,msg.Length()-1);
254  slaveId_ = atoi( slaveID );
255  nSlaves_ = atoi( nSlaves );
256 
257 
258  // Initialise the fit par vectors. Each class that inherits from this one
259  // must implement this sensibly for all vectors specified in clearFitParVectors,
260  // i.e. specify parameter names, initial, min, max and fixed values
261  this->initialise();
262 
263  if (dataFileNameCopy == "") {dataFileNameCopy = "data.root";}
264  if (dataTreeNameCopy == "") {dataTreeNameCopy = "genResults";}
265 
266  if (histFileNameCopy == "") {histFileNameCopy = "parInfo.root";}
267  if (tableFileNameCopy == "") {tableFileNameCopy = "fitResults";}
268 
269  this->fitSlave(dataFileNameCopy, dataTreeNameCopy);
270 
271  cout << "Fit slave " << slaveId_ << " has finished successfully" << std::endl;
272 }
273 
275 {
276  cout << "\n\nWaiting for connection with " << nSlaves_ << " workers...\n\n" << std::endl;
277 
278  //initialize socket connection, then accept a connection and return a full-duplex communication socket.
279  socketMonitor_ = new TMonitor();
280 
281  TServerSocket *ss = new TServerSocket(9090, kTRUE);
282  sSlaves_.resize(nSlaves_);
283  for ( UInt_t i(0); i<nSlaves_; ++i ) {
284  sSlaves_[i] = ss->Accept();
285  }
286 
287  // tell the clients to start
288  for ( UInt_t i(0); i<nSlaves_; ++i ) {
289  TString msg = "Slave ";
290  msg += i;
291  msg += "/";
292  msg += nSlaves_;
293  sSlaves_[i]->Send(msg);
294  socketMonitor_->Add(sSlaves_[i]);
295  cout << "Added slave " << i<<endl;
296  }
297 
298  cout << "Now start" << endl;
299 
300  ss->Close();
301  delete ss;
302 }
303 
304 void LauAbsFitModel::doSFit( const TString& sWeightBranchName, Double_t scaleFactor )
305 {
306  if ( sWeightBranchName == "" ) {
307  cerr << "WARNING in LauAbsFitModel::doSFit : sWeight branch name is empty string, not setting-up sFit." << endl;
308  return;
309  }
310 
311  doSFit_ = kTRUE;
312  sWeightBranchName_ = sWeightBranchName;
313  sWeightScaleFactor_ = scaleFactor;
314 }
315 
316 void LauAbsFitModel::setBkgndClassNames( const std::vector<TString>& names )
317 {
318  if ( !bkgndClassNames_.empty() ) {
319  cerr << "WARNING in LauAbsFitModel::setBkgndClassNames : Names already stored, not changing them." << endl;
320  return;
321  }
322 
323  UInt_t nBkgnds = names.size();
324  for ( UInt_t i(0); i < nBkgnds; ++i ) {
325  bkgndClassNames_.insert( std::make_pair( i, names[i] ) );
326  }
327 
328  this->setupBkgndVectors();
329 }
330 
331 Bool_t LauAbsFitModel::validBkgndClass( const TString& className ) const
332 {
333  if ( bkgndClassNames_.empty() ) {
334  return kFALSE;
335  }
336 
337  Bool_t found(kFALSE);
338  for ( LauBkgndClassMap::const_iterator iter = bkgndClassNames_.begin(); iter != bkgndClassNames_.end(); ++iter ) {
339  if ( iter->second == className ) {
340  found = kTRUE;
341  break;
342  }
343  }
344 
345  return found;
346 }
347 
348 UInt_t LauAbsFitModel::bkgndClassID( const TString& className ) const
349 {
350  if ( ! this->validBkgndClass( className ) ) {
351  cerr << "ERROR in LauAbsFitModel::bkgndClassID : Request for ID for invalid background class \"" << className << "\"." << endl;
352  return (bkgndClassNames_.size() + 1);
353  }
354 
355  UInt_t bgID(0);
356  for ( LauBkgndClassMap::const_iterator iter = bkgndClassNames_.begin(); iter != bkgndClassNames_.end(); ++iter ) {
357  if ( iter->second == className ) {
358  bgID = iter->first;
359  break;
360  }
361  }
362 
363  return bgID;
364 }
365 
366 const TString& LauAbsFitModel::bkgndClassName( UInt_t classID ) const
367 {
368  LauBkgndClassMap::const_iterator iter = bkgndClassNames_.find( classID );
369 
370  if ( iter == bkgndClassNames_.end() ) {
371  cerr << "ERROR in LauAbsFitModel::bkgndClassName : Request for name of invalid background class ID " << classID << "." << endl;
372  return nullString_;
373  }
374 
375  return iter->second;
376 }
377 
379 {
380  cout << "Clearing fit vectors" << endl;
381  fitVars_.clear();
382 }
383 
385 {
386  cout << "Clearing extra ntuple variable vectors" << endl;
387  extraVars_.clear();
388 }
389 
391 {
392  // makes sure each parameter holds its genValue as its current value
393  for (LauParameterPList::iterator iter = fitVars_.begin(); iter != fitVars_.end(); ++iter) {
394  (*iter)->value((*iter)->genValue());
395  }
396  this->propagateParUpdates();
397 }
398 
399 void LauAbsFitModel::writeSPlotData(const TString& fileName, const TString& treeName, Bool_t storeDPEfficiency, const TString& verbosity)
400 {
401  if (this->writeSPlotData()) {
402  cerr << "ERROR in LauAbsFitModel::writeSPlotData : Already have an sPlot ntuple setup, not doing it again." << endl;
403  return;
404  }
405  writeSPlotData_ = kTRUE;
406  sPlotFileName_ = fileName;
407  sPlotTreeName_ = treeName;
408  sPlotVerbosity_ = verbosity;
409  storeDPEff_ = storeDPEfficiency;
410 }
411 
412 // TODO : histFileName isn't used here at the moment but could be used for
413 // storing the values of the parameters used in the generation.
414 // These could then be read and used for setting the "true" values
415 // in a subsequent fit.
416 void LauAbsFitModel::generate(const TString& dataFileName, const TString& dataTreeName, const TString& /*histFileName*/, const TString& tableFileNameBase)
417 {
418  // Create the ntuple for storing the results
419  cout << "Creating generation ntuple." << endl;
420  if (genNtuple_ != 0) {delete genNtuple_; genNtuple_ = 0;}
421  genNtuple_ = new LauGenNtuple(dataFileName,dataTreeName);
422 
423  // add branches for storing the experiment number and the number of
424  // the event within the current experiment
425  this->addGenNtupleIntegerBranch("iExpt");
426  this->addGenNtupleIntegerBranch("iEvtWithinExpt");
427  this->setupGenNtupleBranches();
428 
429  // Start the cumulative timer
430  cumulTimer_.Start();
431 
432  Bool_t genOK(kTRUE);
433  do {
434  // Loop over the number of experiments
435  for (iExpt_ = firstExpt_; iExpt_ < (firstExpt_+nExpt_); ++iExpt_) {
436 
437  // Start the timer to see how long each experiment takes to generate
438  timer_.Start();
439 
440  // Store the experiment number in the ntuple
441  this->setGenNtupleIntegerBranchValue("iExpt",iExpt_);
442 
443  // Do the generation for this experiment
444  cout << "Generating experiment number " << iExpt_ << endl;
445  genOK = this->genExpt();
446 
447  // Stop the timer and see how long the program took so far
448  timer_.Stop();
449  timer_.Print();
450 
451  if (!genOK) {
452  // delete and recreate an empty tree
454 
455  // then break out of the experiment loop
456  cerr << "ERROR in LauAbsFitModel::generate : Problem in toy MC generation. Starting again with updated parameters..." << endl;
457  break;
458  }
459 
460  if (this->writeLatexTable()) {
461  TString tableFileName(tableFileNameBase);
462  tableFileName += "_";
463  tableFileName += iExpt_;
464  tableFileName += ".tex";
465  this->writeOutTable(tableFileName);
466  }
467 
468  } // Loop over number of experiments
469  } while (!genOK);
470 
471  // Print out total timing info.
472  cumulTimer_.Stop();
473  cout << "Finished generating all experiments." << endl;
474  cout << "Cumulative timing:" << endl;
475  cumulTimer_.Print();
476 
477  // Build the event index
478  cout << "Building experiment:event index." << endl;
479  // TODO - can test this return value?
480  //Int_t nIndexEntries =
481  genNtuple_->buildIndex("iExpt","iEvtWithinExpt");
482 
483  // Write out toy MC ntuple
484  cout << "Writing data to file " << dataFileName << "." << endl;
486 }
487 
489 {
491 }
492 
494 {
496 }
497 
499 {
500  genNtuple_->setIntegerBranchValue(name,value);
501 }
502 
504 {
505  genNtuple_->setDoubleBranchValue(name,value);
506 }
507 
509 {
510  return genNtuple_->getIntegerBranchValue(name);
511 }
512 
514 {
515  return genNtuple_->getDoubleBranchValue(name);
516 }
517 
519 {
521 }
522 
524 {
526 }
527 
529 {
531 }
532 
534 {
535  sPlotNtuple_->setIntegerBranchValue(name,value);
536 }
537 
539 {
540  sPlotNtuple_->setDoubleBranchValue(name,value);
541 }
542 
544 {
546 }
547 
548 void LauAbsFitModel::fit(const TString& dataFileName, const TString& dataTreeName, const TString& histFileName, const TString& tableFileNameBase)
549 {
550  // Routine to perform the total fit.
551 
552  // Check whether we're going to use asymmetric errors
553  // This boolean can be changed with the useAsymmFitErrors(Bool_t) function
554  if (useAsymmFitErrors_ == kTRUE) {
555  cout << "We are going to use MINOS to calculate the asymmetric fit errors." << endl;
556  cout << "This will in general significantly increase the CPU time required for fitting." << endl;
557  cout << "Use setCalcAsymmFitErrors(kFALSE) if you want disable this feature." << endl;
558  }
559 
560  cout << "First experiment = " << firstExpt_ << endl;
561  cout << "Number of experiments = " << nExpt_ << endl;
562 
563  // Start the cumulative timer
564  cumulTimer_.Start();
565 
566  numberOKFits_ = 0, numberBadFits_ = 0;
567  fitStatus_ = -1;
568 
569  // Create and setup the fit results ntuple
570  cout << "Creating fit ntuple." << endl;
571  if (fitNtuple_ != 0) {delete fitNtuple_; fitNtuple_ = 0;}
572  fitNtuple_ = new LauFitNtuple(histFileName);
573 
574  // Create and setup the sPlot ntuple
575  if (this->writeSPlotData()) {
576  cout << "Creating sPlot ntuple." << endl;
577  if (sPlotNtuple_ != 0) {delete sPlotNtuple_; sPlotNtuple_ = 0;}
579  this->setupSPlotNtupleBranches();
580  }
581 
582  // This reads in the given dataFile and creates an input
583  // fit data tree that stores them for all events and experiments.
584  Bool_t dataOK = this->cacheFitData(dataFileName,dataTreeName);
585  if (!dataOK) {
586  cerr << "ERROR in LauAbsFitModel::fit : Problem caching the fit data." << endl;
587  gSystem->Exit(EXIT_FAILURE);
588  }
589 
590  // Loop over the number of experiments
591  for (iExpt_ = firstExpt_; iExpt_ < (firstExpt_+nExpt_); ++iExpt_) {
592 
593  // Start the timer to see how long each fit takes
594  timer_.Start();
595 
598 
599  if (this->eventsPerExpt() < 1) {
600  cerr << "ERROR in LauAbsFitModel::fit : Zero events in experiment " << iExpt_ << ", skipping..." << endl;
601  timer_.Stop();
602  continue;
603  }
604 
605  // Now the sub-classes must implement whatever they need to do
606  // to cache any more input fit data they need in order to
607  // calculate the likelihoods during the fit.
608  // They need to use the inputFitData_ tree as input. For example,
609  // inputFitData_ contains m13Sq and m23Sq. The appropriate fit model
610  // then caches the resonance dynamics for the signal model, as
611  // well as the background likelihood values in the Dalitz plot
612  this->cacheInputFitVars();
613  if ( this->doSFit() ) {
614  this->cacheInputSWeights();
615  }
616 
617  // Do the fit for this experiment
618  this->fitExpt();
619 
620  // Write the results into the ntuple
621  this->finaliseFitResults(tableFileNameBase);
622 
623  // Store the per-event likelihood values
624  if ( this->writeSPlotData() ) {
625  this->storePerEvtLlhds();
626  }
627 
628  // Create a toy MC sample for the 1st successful experiment
629  // using the fitted parameters so that the user can compare
630  // the fit to the actual data. The toy MC stats are a
631  // factor of 10 higher than the number of events specified
632  // via setNEvGen. This is to reduce the statistical
633  // fluctuations for the toy MC data.
634  if (compareFitData_ == kTRUE && fitStatus_ == 3) {
636  compareFitData_ = kFALSE; // only do this for the first successful experiment
637  }
638 
639  // Stop the timer and see how long the program took so far
640  timer_.Stop();
641  timer_.Print();
642 
643  // Keep track of how many fits worked or failed
644  // NB values of fitStatus_ now indicate the status of the error matrix:
645  // 0= not calculated at all
646  // 1= approximation only, not accurate
647  // 2= full matrix, but forced positive-definite
648  // 3= full accurate covariance matrix
649  if (fitStatus_ == 3) {
650  numberOKFits_++;
651  } else {
652  numberBadFits_++;
653  }
654 
655  } // Loop over number of experiments
656 
657  // Print out total timing info.
658  cumulTimer_.Stop();
659  cout << " Cumulative timing:" << endl;
660  cumulTimer_.Print();
661 
662  // Print out stats on OK fits.
663  cout << "Number of OK Fits = " << numberOKFits_ << endl;
664  cout << "Number of Failed Fits = " << numberBadFits_ << endl;
665  Double_t fitEff(0.0);
666  if (nExpt_ != 0) {fitEff = numberOKFits_/(1.0*nExpt_);}
667  cout << "Fit efficiency = " << fitEff*100.0 << "%." << endl;
668 
669  // Write out any fit results (ntuples etc...).
670  this->writeOutAllFitResults();
671  if ( this->writeSPlotData() ) {
672  this->calculateSPlotData();
673  }
674 }
675 
676 void LauAbsFitModel::fitMaster(const TString& dataFileName, const TString& dataTreeName, const TString& histFileName, const TString& tableFileNameBase)
677 {
678  this->fit( dataFileName, dataTreeName, histFileName, tableFileNameBase );
679 }
680 
681 void LauAbsFitModel::fitSlave(const TString& dataFileName, const TString& dataTreeName)
682 {
683 
684  // This reads in the given dataFile and creates an input
685  // fit data tree that stores them for all events and experiments.
686  Bool_t dataOK = this->cacheFitData(dataFileName,dataTreeName);
687  if (!dataOK) {
688  cerr << "ERROR in LauAbsFitModel::fit : Problem caching the fit data." << endl;
689  gSystem->Exit(EXIT_FAILURE);
690  }
691 
694 
695  if (this->eventsPerExpt() < 1) {
696  cerr << "ERROR in LauAbsFitModel::fit : Zero events in experiment " << firstExpt_ << ", aborting..." << endl;
697  gSystem->Exit(EXIT_FAILURE);
698  }
699 
700  this->cacheInputFitVars();
701  if ( this->doSFit() ) {
702  this->cacheInputSWeights();
703  }
704  std::cout << "\nSlave " << slaveId_ << " ready...\n" << std::endl;
705 
706 
707  // infinite loop waiting messages from server with new parameters
708  TMessage messageToMaster_(kMESS_OBJECT);
709  char str[128];
710 
711  while (1) {
713  if (!messageFromMaster_) break;
714 
715  if (messageFromMaster_->What() == kMESS_STRING) {
716  messageFromMaster_->ReadString(str, 128);
717  TString msg( str );
718  std::cout << "Message from master: " << msg << std::endl;
719  if (msg == "Finish") {
720  break;
721  }
722  } else if (messageFromMaster_->What() == kMESS_OBJECT) {
723  TVectorD *vecPars_ = dynamic_cast<TVectorD*>( messageFromMaster_->ReadObject( messageFromMaster_->GetClass() ) );
724  if (vecPars_) {
725  Double_t * params = vecPars_->GetMatrixArray();
726  UInt_t npar = static_cast<UInt_t>( params[0] );
727  //Int_t iflag = (Int_t)params[1];
728 
729  // We split the data among the slaves.
730  UInt_t nEventsPerSlave = this->eventsPerExpt() / nSlaves_;
731  UInt_t start( nEventsPerSlave*slaveId_ );
732  UInt_t end( nEventsPerSlave*(slaveId_+1) );
733  if ( slaveId_ == (nSlaves_-1) ) {
734  end = this->eventsPerExpt();
735  }
736 
737  for (UInt_t i=2;i<npar+2;i++) {
738  if (!fitVars_[i-2]->fixed()) {
739  fitVars_[i-2]->value(params[i]);
740  }
741  }
742  this->propagateParUpdates();
743 
744  Double_t func = this->getLogLikelihood(start,end);
745 
746  // send the function in the same vector
747  vecPars_[0] = func;
748  messageToMaster_.Reset();
749  messageToMaster_.WriteObject(vecPars_);
750  sMaster_->Send(messageToMaster_);
751  delete vecPars_;
752  }
753  } else {
754  std::cerr << "ERROR in LauAbsFitModel::fitSlave : *** Unexpected type of message from master ***" << std::endl;
755  gSystem->Exit(EXIT_FAILURE);
756  }
757 
758  }
759 }
760 
761 Bool_t LauAbsFitModel::cacheFitData(const TString& dataFileName, const TString& dataTreeName)
762 {
763  // From the input data stream, store the variables into the
764  // internal tree inputFitData_ that can be used by the sub-classes
765  // in calculating their likelihood functions for the fit
766  delete inputFitData_;
767  inputFitData_ = new LauFitDataTree(dataFileName,dataTreeName);
768  Bool_t dataOK = inputFitData_->findBranches();
769 
770  if (!dataOK) {
771  delete inputFitData_; inputFitData_ = 0;
772  }
773 
774  return dataOK;
775 }
776 
778 {
779  Bool_t hasBranch = inputFitData_->haveBranch( sWeightBranchName_ );
780  if ( ! hasBranch ) {
781  std::cerr << "ERROR in LauAbsFitModel::cacheInputSWeights : Input data does not contain variable \"" << sWeightBranchName_ << "\".\n";
782  std::cerr << " : Turning off sFit!" << std::endl;
783  doSFit_ = kFALSE;
784  return;
785  }
786 
787  UInt_t nEvents = this->eventsPerExpt();
788  sWeights_.clear();
789  sWeights_.reserve( nEvents );
790 
791  for (UInt_t iEvt = 0; iEvt < nEvents; ++iEvt) {
792 
793  const LauFitData& dataValues = inputFitData_->getData(iEvt);
794 
795  LauFitData::const_iterator iter = dataValues.find( sWeightBranchName_ );
796 
797  sWeights_.push_back( iter->second * sWeightScaleFactor_ );
798  }
799 }
800 
802 {
803  // Routine to perform the actual fit for the given experiment, and store
804  // the results in the histograms if required (doHist == kTRUE).
805 
806  // Reset the worst likelihood found to its catch-all value
807  worstLogLike_ = DBL_MAX;
808 
809  // Hook the external likelihood function to this LauFitter::fitter() and this class.
810  LauFitter::fitter()->SetFCN(logLikeFun);
811  LauFitter::fitter()->SetObjectFit(this);
812 
813  // Clear any stored parameters etc... before using
814  LauFitter::fitter()->Clear();
815 
816  // Define the default relative error
817  const Double_t defaultError(0.01);
818 
819  // Update initial fit parameters if required (e.g. if using random numbers).
820  this->checkInitFitParams();
821 
822  nParams_ = fitVars_.size();
823  cout << "Total number of parameters = " << nParams_ << endl;
824  cout << "Setting fit parameters" << endl;
825 
826  // Set-up the parameters to be fit
827  for (UInt_t i = 0; i < nParams_; i++) {
828  TString name = fitVars_[i]->name();
829  Double_t initVal = fitVars_[i]->initValue();
830  Double_t initErr = fitVars_[i]->error();
831  if ( initVal == 0.0 ) {
832  initErr = defaultError;
833  } else if ( TMath::Abs(initErr/initVal) < 1e-6 ) {
834  initErr = TMath::Abs(defaultError * initVal);
835  }
836  Double_t minVal = fitVars_[i]->minValue();
837  Double_t maxVal = fitVars_[i]->maxValue();
838  Bool_t secondStage = fitVars_[i]->secondStage();
839  if (this->twoStageFit() && secondStage == kTRUE) {
840  fitVars_[i]->fixed(kTRUE);
841  }
842  Bool_t fixVar = fitVars_[i]->fixed();
843 
844  cout << "Setting parameter " << i << " called " << name << " to have initial value " << initVal << ", error " << initErr << " and range " << minVal << " to " << maxVal << endl;
845  LauFitter::fitter()->SetParameter(i, name, initVal, initErr, minVal, maxVal);
846 
847  // Fix parameter if required
848  if (fixVar == kTRUE) {
849  cout << "Fixing parameter " << i << std::endl;
850  LauFitter::fitter()->FixParameter(i);
851  }
852  }
853 
854  LauParamFixed pred;
855  nFreeParams_ = nParams_ - std::count_if(fitVars_.begin(),fitVars_.end(),pred);
856 
857  // Need to set the "SET ERR" command to +0.5 for +/-1 sigma errors
858  // for maximum likelihood fit. Very important command, otherwise all
859  // extracted errors will be too big, and pull distributions will be too narrow!
860  // TODO - The alternative to this is to make FCN = -2log(L) rather than -log(L)
861  Double_t argL[2];
862  argL[0] = 0.5;
863  fitStatus_ = LauFitter::fitter()->ExecuteCommand("SET ERR", argL, 1);
864 
865  //argL[0] = 0;
866  //fitStatus_ = LauFitter::fitter()->ExecuteCommand("SET STRATEGY", argL, 1);
867 
868  // Now ready for minimization step
869  cout << "\n...Start minimization...\n";
870  Bool_t ok = this->runMinimisation();
871 
872  // If we're doing a two stage fit we can now release (i.e. float)
873  // the 2nd stage parameters and re-fit
874  if (this->twoStageFit()) {
875 
876  if (!ok) {
877  cerr << "ERROR in LauAbsFitModel:fitExpt : Not running second stage fit since first stage failed." << endl;
878  // Need to mark the second stage parameters as free
879  // now so that the fit ntuple doesn't get confused
880  for (UInt_t i = 0; i < nParams_; i++) {
881  Bool_t secondStage = fitVars_[i]->secondStage();
882  if (secondStage == kTRUE) {
883  fitVars_[i]->fixed(kFALSE);
884  }
885  }
886  } else {
887  for (UInt_t i = 0; i < nParams_; i++) {
888  // Release "secondStage" parameters
889  Bool_t secondStage = fitVars_[i]->secondStage();
890  if (secondStage == kTRUE) {
891  fitVars_[i]->fixed(kFALSE);
892  LauFitter::fitter()->ReleaseParameter(i);
893  }
894  // Fix "firstStage" parameters
895  Bool_t firstStage = fitVars_[i]->firstStage();
896  if (firstStage == kTRUE) {
897  fitVars_[i]->fixed(kTRUE);
898  LauFitter::fitter()->FixParameter(i);
899  }
900  }
901  nFreeParams_ = nParams_ - std::count_if(fitVars_.begin(),fitVars_.end(),pred);
902  this->runMinimisation();
903  }
904  }
905 
906  // Store the final fit results and errors into protected internal vectors that
907  // all sub-classes can use within their own finalFitResults implementation
908  // used below (e.g. putting them into an ntuple in a root file)
909  for (UInt_t i = 0; i < nParams_; i++) {
910  // Get the value and errors from MINUIT
911  Double_t value = LauFitter::fitter()->GetParameter(i);
912  Double_t error(0.0);
913  Double_t negError(0.0);
914  Double_t posError(0.0);
915  Double_t globalcc(0.0);
916  if (useAsymmFitErrors_) {
917  LauFitter::fitter()->GetErrors(i, posError, negError, error, globalcc);
918  } else {
919  error = LauFitter::fitter()->GetParError(i);
920  }
921  fitVars_[i]->valueAndErrors(value, error, negError, posError);
922 
923  // release "firstStage" parameters so fit results are stored
924  Bool_t firstStage = fitVars_[i]->firstStage();
925  if (firstStage == kTRUE) {
926  fitVars_[i]->fixed(kFALSE);
927  }
928  }
929 }
930 
932 {
933  Double_t arglist[2];
934  arglist[0] = 1000*nParams_; // maximum iterations
935  arglist[1] = 0.05; // tolerance -> min EDM = 0.001*tolerance (0.05)
936  fitStatus_ = LauFitter::fitter()->ExecuteCommand("MIGRAD", arglist, 2);
937 
938  // Dummy variables - need to feed them to the function
939  // used for getting NLL and error matrix status
940  Double_t edm, errdef;
941  Int_t nvpar, nparx;
942 
943  if (fitStatus_ != 0) {
944 
945  cerr << "Error in minimising loglike." << endl;
946 
947  } else {
948 
949  // Check that the error matrix is ok
950  NLL_ = 0.0;
951  fitStatus_ = LauFitter::fitter()->GetStats(NLL_, edm, errdef, nvpar, nparx);
952  cout << "Error matrix status after MIGRAD is: " << fitStatus_ << endl;
953  // 0= not calculated at all
954  // 1= approximation only, not accurate
955  // 2= full matrix, but forced positive-definite
956  // 3= full accurate covariance matrix
957 
958  // Fit result was OK. Now get the more precise errors.
959  fitStatus_ = LauFitter::fitter()->ExecuteCommand("HESSE", arglist, 1);
960 
961  if (fitStatus_ != 0) {
962 
963  cerr << "Error in Hesse routine." << endl;
964 
965  } else {
966 
967  // Check that the error matrix is ok
968  NLL_ = 0.0;
969  fitStatus_ = LauFitter::fitter()->GetStats(NLL_, edm, errdef, nvpar, nparx);
970  cout << "Error matrix status after HESSE is: " << fitStatus_ << endl;
971  // 0= not calculated at all
972  // 1= approximation only, not accurate
973  // 2= full matrix, but forced positive-definite
974  // 3= full accurate covariance matrix
975 
976  // Symmetric errors and eror matrix were OK.
977  // Get asymmetric errors if asked for.
978  if (useAsymmFitErrors_ == kTRUE) {
979  withinMINOS_ = kTRUE;
980  fitStatus_ = LauFitter::fitter()->ExecuteCommand("MINOS", arglist, 1);
981  withinMINOS_ = kFALSE;
982  if (fitStatus_ != 0) {
983  cerr << "Error in Minos routine." << endl;
984  }
985  }
986  }
987  }
988 
989  // Print results
990  NLL_ = 0.0;
991  fitStatus_ = LauFitter::fitter()->GetStats(NLL_, edm, errdef, nvpar, nparx);
992  // If the error matrix is not accurate, fail the fit
993  cout << "Final error matrix status is: " << fitStatus_ << endl;
994  // 0= not calculated at all
995  // 1= approximation only, not accurate
996  // 2= full matrix, but forced positive-definite
997  // 3= full accurate covariance matrix
998 
999  LauFitter::fitter()->PrintResults(3, NLL_);
1000 
1001  return (fitStatus_ == 3);
1002 }
1003 
1005 {
1006  // Write out histograms at end
1007  if (fitNtuple_ != 0) {
1009  }
1010 }
1011 
1013 {
1014  if (sPlotNtuple_ != 0) {
1017  LauSPlot splot(sPlotNtuple_->fileName(), sPlotNtuple_->treeName(), this->firstExpt(), this->nExpt(),
1018  this->variableNames(), this->freeSpeciesNames(), this->fixdSpeciesNames(), this->twodimPDFs(),
1019  this->splitSignal(), this->scfDPSmear());
1021  splot.writeOutResults();
1022  }
1023 }
1024 
1025 void LauAbsFitModel::compareFitData(UInt_t toyMCScale, const TString& mcFileName, const TString& tableFileName, Bool_t poissonSmearing)
1026 {
1027  compareFitData_ = kTRUE;
1028  fitToyMCScale_ = toyMCScale;
1029  fitToyMCFileName_ = mcFileName;
1030  fitToyMCTableName_ = tableFileName;
1031  fitToyMCPoissonSmear_ = poissonSmearing;
1032 }
1033 
1034 void LauAbsFitModel::createFitToyMC(const TString& mcFileName, const TString& tableFileName)
1035 {
1036  // Create a toy MC sample so that the user can eventually
1037  // compare the "fitted" result with the data
1038  // Generate more toy MC to reduce statistical fluctuations. Use the
1039  // rescaling value fitToyMCScale_.
1040 
1041  // Store the info on the number of experiments, first expt and current expt
1042  UInt_t oldNExpt(this->nExpt());
1043  UInt_t oldFirstExpt(this->firstExpt());
1044  UInt_t oldIExpt(iExpt_);
1045 
1046  // Turn off Poisson smearing if required
1047  Bool_t poissonSmearing(this->doPoissonSmearing());
1049 
1050  // Turn off embedding, since we need toy MC, not reco'ed events
1051  Bool_t enableEmbeddingOrig(this->enableEmbedding());
1052  this->enableEmbedding(kFALSE);
1053 
1054  // Need to make sure that the generation of the DP co-ordinates is
1055  // switched on if any of our PDFs depend on it
1056  Bool_t origUseDP = this->useDP();
1057  if ( this->pdfsDependOnDP() && !origUseDP ) {
1058  this->useDP( kTRUE );
1059  this->initialiseDPModels();
1060  }
1061 
1062  // Generate the toy MC
1063  cout << "Generating toy MC in " << mcFileName << " to compare fit with data..." << endl;
1064  cout << "Number of experiments to generate = " << this->nExpt() << ", which is a factor of "
1065  <<fitToyMCScale_ << " higher than that originally specified. "
1066  <<"This is to allow the toy MC to be made with reduced statistical "
1067  <<"fluctuations." << endl;
1068 
1069  // Set the genValue of each parameter to its current (fitted) value
1070  // but first store the original genValues for restoring later
1071  vector<Double_t> origGenValues; origGenValues.reserve(nParams_);
1072  for (LauParameterPList::iterator iter = fitVars_.begin(); iter != fitVars_.end(); ++iter) {
1073  origGenValues.push_back((*iter)->genValue());
1074  (*iter)->genValue((*iter)->value());
1075  }
1076 
1077  // If we're asked to generate more than 100 experiments then split it
1078  // up into multiple files since otherwise can run into memory issues
1079  // when building the index
1080 
1081  UInt_t totalExpts = oldNExpt * fitToyMCScale_;
1082  if ( totalExpts > 100 ) {
1083  UInt_t nFiles = totalExpts/100;
1084  if ( totalExpts%100 ) {
1085  nFiles += 1;
1086  }
1087  for ( UInt_t iFile(0); iFile < nFiles; ++iFile ) {
1088 
1089  UInt_t firstExp( iFile*100 );
1090 
1091  // Set number of experiments and first experiment to generate
1092  UInt_t nExp = ((firstExp + 100)>totalExpts) ? totalExpts-firstExp : 100;
1093  this->setNExpts(nExp, firstExp);
1094 
1095  // Create a unique filename and generate the events
1096  TString extraname = "_file";
1097  extraname += iFile;
1098  TString filename(mcFileName);
1099  filename.Insert( filename.Index("."), extraname );
1100  this->generate(filename, "genResults", "dummy.root", tableFileName);
1101  }
1102  } else {
1103  // Set number of experiments to new value
1104  this->setNExpts(oldNExpt*fitToyMCScale_, 0);
1105  // Generate the toy
1106  this->generate(mcFileName, "genResults", "dummy.root", tableFileName);
1107  }
1108 
1109  // Reset number of experiments to original value
1110  iExpt_ = oldIExpt;
1111  this->setNExpts(oldNExpt, oldFirstExpt);
1112 
1113  // Restore the Poisson smearing to its former value
1114  this->doPoissonSmearing(poissonSmearing);
1115 
1116  // Restore the embedding status
1117  this->enableEmbedding(enableEmbeddingOrig);
1118 
1119  // Restore "useDP" to its former status
1120  this->useDP( origUseDP );
1121 
1122  // Restore the original genValue to each parameter
1123  for (UInt_t i(0); i<nParams_; ++i) {
1124  fitVars_[i]->genValue(origGenValues[i]);
1125  }
1126 
1127  cout << "Finished in createFitToyMC." << endl;
1128 }
1129 
1131 {
1132  // Calculate the total negative log-likelihood over all events.
1133  // This function assumes that the fit parameters and data tree have
1134  // already been set-up correctly.
1135 
1136  // Loop over the data points to calculate the log likelihood
1137  Double_t logLike(0.0);
1138 
1139  // Loop over the number of events in this experiment
1140  Bool_t ok(kTRUE);
1141  for (UInt_t iEvt = 0; iEvt < this->eventsPerExpt(); ++iEvt) {
1142 
1143  Double_t likelihood = this->getTotEvtLikelihood(iEvt);
1144 
1145  if (likelihood > DBL_MIN) { // Is the likelihood zero?
1146  Double_t evtLogLike = TMath::Log(likelihood);
1147  if ( doSFit_ ) {
1148  evtLogLike *= sWeights_[iEvt];
1149  }
1150  logLike += evtLogLike;
1151  } else {
1152  ok = kFALSE;
1153  cerr << "WARNING in LauAbsFitModel::getTotNegLogLikelihood : Strange likelihood value for event " << iEvt << ": " << likelihood << "\n";
1154  this->printEventInfo(iEvt);
1155  this->printVarsInfo(); //Write the values of the floated variables for which the likelihood is zero
1156  break;
1157  }
1158  }
1159 
1160  // Include the Poisson term in the extended likelihood if required
1161  if (this->doEMLFit()) {
1162  logLike -= this->getEventSum();
1163  }
1164 
1165  if (!ok) {
1166  cerr << " : Returning worst NLL found so far to force MINUIT out of this region." << endl;
1167  logLike = worstLogLike_;
1168  } else if (logLike < worstLogLike_) {
1169  worstLogLike_ = logLike;
1170  }
1171 
1172  Double_t totNegLogLike = -logLike;
1173  return totNegLogLike;
1174 }
1175 
1176 Double_t LauAbsFitModel::getLogLikelihood( UInt_t iStart, UInt_t iEnd )
1177 {
1178  // Calculate the total negative log-likelihood over all events.
1179  // This function assumes that the fit parameters and data tree have
1180  // already been set-up correctly.
1181 
1182  // Loop over the data points to calculate the log likelihood
1183  Double_t logLike(0.0);
1184 
1185  // Loop over the number of events in this experiment
1186  Bool_t ok(kTRUE);
1187  for (UInt_t iEvt = iStart; iEvt < iEnd; ++iEvt) {
1188 
1189  Double_t likelihood = this->getTotEvtLikelihood(iEvt);
1190 
1191  if (likelihood > DBL_MIN) { // Is the likelihood zero?
1192  Double_t evtLogLike = TMath::Log(likelihood);
1193  if ( doSFit_ ) {
1194  evtLogLike *= sWeights_[iEvt];
1195  }
1196  logLike += evtLogLike;
1197  } else {
1198  ok = kFALSE;
1199  cerr << "WARNING in LauAbsFitModel::getLogLikelihood : Strange likelihood value for event " << iEvt << ": " << likelihood << "\n";
1200  this->printEventInfo(iEvt);
1201  this->printVarsInfo(); //Write the values of the floated variables for which the likelihood is zero
1202  break;
1203  }
1204  }
1205 
1206  if (!ok) {
1207  cerr << " : Returning worst NLL found so far to force MINUIT out of this region." << endl;
1208  logLike = worstLogLike_;
1209  } else if (logLike < worstLogLike_) {
1210  worstLogLike_ = logLike;
1211  }
1212 
1213  return logLike;
1214 }
1215 
1216 void LauAbsFitModel::setParsFromMinuit(Double_t* par, Int_t npar)
1217 {
1218  // This function sets the internal parameters based on the values
1219  // that Minuit is using when trying to minimise the total likelihood function.
1220 
1221  // MINOS reports one fewer free parameter than there actually is,
1222  // so increment npar so the following check doesn't fail
1223  if ( withinMINOS_ ) {
1224  ++npar;
1225  }
1226 
1227  if (static_cast<UInt_t>(npar) != nFreeParams_) {
1228  cerr << "ERROR in LauAbsFitModel::setParsFromMinuit : Unexpected number of free parameters: " << npar << ".\n";
1229  cerr << " Expected: " << nFreeParams_ << ".\n" << endl;
1230  gSystem->Exit(EXIT_FAILURE);
1231  }
1232 
1233  // Despite npar being the number of free parameters
1234  // the par array actually contains all the parameters,
1235  // free and floating...
1236  // Update all the floating ones with their new values.
1237  for (UInt_t i(0); i<nParams_; ++i) {
1238  if (!fitVars_[i]->fixed()) {
1239  fitVars_[i]->value(par[i]);
1240  }
1241  }
1242 
1243  this->propagateParUpdates();
1244 }
1245 
1247 {
1248  UInt_t nParsAdded(0);
1249  for (LauPdfList::iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) {
1250  LauAbsPdf* pdf = (*pdf_iter);
1251  if ( pdf->isDPDependent() ) {
1252  this->pdfsDependOnDP( kTRUE );
1253  }
1254  LauParameterPList& pars = pdf->getParameters();
1255  for (LauParameterPList::iterator pars_iter = pars.begin(); pars_iter != pars.end(); ++pars_iter) {
1256  if ( !(*pars_iter)->clone() && ( !(*pars_iter)->fixed() ||
1257  (this->twoStageFit() && ( (*pars_iter)->secondStage() || (*pars_iter)->firstStage())) ) ) {
1258  fitVars_.push_back(*pars_iter);
1259  ++nParsAdded;
1260  }
1261  }
1262  }
1263  return nParsAdded;
1264 }
1265 
1267 {
1268  for (LauPdfList::iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) {
1269  (*pdf_iter)->updatePulls();
1270  }
1271 }
1272 
1273 void LauAbsFitModel::printFitParameters(const LauPdfList& pdfList, std::ostream& fout) const
1274 {
1275  LauPrint print;
1276  for (LauPdfList::const_iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) {
1277  const LauParameterPList& pars = (*pdf_iter)->getParameters();
1278  for (LauParameterPList::const_iterator pars_iter = pars.begin(); pars_iter != pars.end(); ++pars_iter) {
1279  if (!(*pars_iter)->clone()) {
1280  fout << (*pars_iter)->name() << " & $";
1281  print.printFormat(fout, (*pars_iter)->value());
1282  if ((*pars_iter)->fixed() == kTRUE) {
1283  fout << "$ (fixed) \\\\";
1284  } else {
1285  fout << " \\pm ";
1286  print.printFormat(fout, (*pars_iter)->error());
1287  fout << "$ \\\\" << endl;
1288  }
1289  }
1290  }
1291  }
1292 }
1293 
1295 {
1296  for (LauPdfList::iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) {
1297  (*pdf_iter)->cacheInfo(theData);
1298  }
1299 }
1300 
1301 Double_t LauAbsFitModel::prodPdfValue(LauPdfList& pdfList, UInt_t iEvt)
1302 {
1303  Double_t pdfVal = 1.0;
1304  for (LauPdfList::iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) {
1305  (*pdf_iter)->calcLikelihoodInfo(iEvt);
1306  pdfVal *= (*pdf_iter)->getLikelihood();
1307  }
1308  return pdfVal;
1309 }
1310 
1311 void LauAbsFitModel::printEventInfo(UInt_t iEvt) const
1312 {
1313  const LauFitData& data = inputFitData_->getData(iEvt);
1314  cerr << " : Input data values for this event:" << endl;
1315  for (LauFitData::const_iterator iter = data.begin(); iter != data.end(); ++iter) {
1316  cerr << " " << iter->first << " = " << iter->second << endl;
1317  }
1318 }
1319 
1321 {
1322  cerr << " : Current values of fit parameters:" << endl;
1323  for (UInt_t i(0); i<nParams_; ++i) {
1324  cerr << " " << (fitVars_[i]->name()).Data() << " = " << fitVars_[i]->value() << endl;
1325  }
1326 }
1327 
1328 Double_t LauAbsFitModel::calculateLogLike(Int_t npar, Double_t* par, Int_t iflag)
1329 {
1330  // check if we're running a parallel setup or not
1331  if ( socketMonitor_ == 0 ) {
1332  // we're not, so just do things normally
1333  return this->getTotNegLogLikelihood();
1334  }
1335 
1336  TMessage messageToSlave(kMESS_OBJECT);
1337  TSocket *sActual(0);
1338 
1339  //send current parameters to the clients. Must send all the par vector
1340  //In fact, the number npar = number of free parameters, while par has all parameters
1341  if ((UInt_t)npar > nParams_){
1342  std::cerr << "ERROR in LauAbsFitModel::calculateLogLike : Oops, too many parameters, abort" << std::endl;
1343  gSystem->Exit(EXIT_FAILURE);
1344  }
1345 
1346  (*vectorPar_)[0] = nParams_;
1347  (*vectorPar_)[1] = iflag;
1348  for (UInt_t i(0);i<nParams_;i++) (*vectorPar_)[i+2] = par[i];
1349 
1350  messageToSlave.WriteObject(vectorPar_); //write vector of parameters in message buffer
1351  for ( UInt_t i(0); i<nSlaves_; ++i ) {
1352  sSlaves_[i]->Send(messageToSlave);
1353  }
1354 
1355  Double_t func(0.0);
1356  UInt_t msgsReceived(0);
1357  while (1) {
1358  sActual = socketMonitor_->Select();
1359  sActual->Recv(messageFromSlave_);
1360 
1361  vectorRes_ = dynamic_cast<TVectorD*>( messageFromSlave_->ReadObject( messageFromSlave_->GetClass() ) );
1362  if (!vectorRes_) {
1363  std::cerr << "ERROR in LauAbsFitModel::calculateLogLike : vector parameter is NULL !" << std::endl;
1364  return 0.0;
1365  }
1366 
1367  func += (*vectorRes_)[0];
1368  ++msgsReceived;
1369 
1370  if (msgsReceived == nSlaves_) {
1371  break;
1372  }
1373  }
1374 
1375  // Include the Poisson term in the extended likelihood if required
1376  if (this->doEMLFit()) {
1377  func -= this->getEventSum();
1378  }
1379 
1380  Double_t negLogLike = -func;
1381  return negLogLike;
1382 }
1383 
1384 // Definition of the fitting function for Minuit
1385 void logLikeFun(Int_t& npar, Double_t* /*first-derivatives*/, Double_t& f, Double_t* par, Int_t iflag)
1386 {
1387  // Routine that specifies the negative log-likelihood function for the fit.
1388  // Used by the MINUIT minimising code.
1389 
1390  LauAbsFitModel* theModel = dynamic_cast<LauAbsFitModel*>(LauFitter::fitter()->GetObjectFit());
1391 
1392  // Set the internal parameters for this model using parameters from Minuit (pars):
1393  theModel->setParsFromMinuit( par, npar );
1394 
1395  // Set the value of f to be the total negative log-likelihood for the data sample.
1396  f = theModel->calculateLogLike( npar, par, iflag );
1397 }
1398 
Class for defining the SPlot technique.
Definition: LauSPlot.hh:55
TString fitToyMCTableName_
The output table name for Toy MC.
virtual void addSPlotNtupleIntegerBranch(const TString &name)
Add a branch to the sPlot tree for storing an integer.
Bool_t findBranches()
Find all of the branches in the tree.
virtual Bool_t genExpt()=0
The method that actually generates the toy MC events for the given experiment.
void createFitToyMC(const TString &mcFileName, const TString &tableFileName)
Create a toy MC sample from the fitted parameters.
UInt_t firstExpt() const
Obtain the number of the first experiment.
File containing declaration of LauFitDataTree class.
virtual Double_t getEventSum() const =0
Returns the sum of the expected events over all hypotheses; used in the EML fit scenario.
Bool_t fixed() const
Check whether the parameter is fixed or floated.
void run(const TString &applicationCode, const TString &dataFileName="", const TString &dataTreeName="", const TString &histFileName="", const TString &tableFileName="")
Start the toy generation / fitting.
UInt_t eventsPerExpt() const
Obtain the total number of events in the current experiment.
void setGenValues()
Make sure all parameters hold their genValue as the current value.
TString sPlotFileName_
The name of the sPlot file.
virtual void writeOutTable(const TString &outputFile)=0
Write the latex table.
Bool_t writeLatexTable() const
Determine whether writing out of the latex table is enabled.
Bool_t doSFit_
Option to perfom the sFit.
Bool_t writeSPlotData() const
Determine whether the sPlot data is to be written out.
TStopwatch timer_
The fit timer.
UInt_t nParams_
The number of fit parameters.
File containing declaration of LauParamFixed class.
Double_t getTotNegLogLikelihood()
Retrieves the total negative log-likelihood.
Bool_t storeDPEff_
Option to store DP efficiencies in the sPlot ntuple.
void cacheInfo(LauPdfList &pdfList, const LauFitDataTree &theData)
Have all PDFs in the list cache the data.
Bool_t twoStageFit() const
Determine whether the two-stage fit is enabled.
void runCalculations(const TString &option="q")
Method to calculate the sWeights and cN coeffs.
Definition: LauSPlot.cc:626
std::vector< TSocket * > sSlaves_
Sockets for each of the slaves to communicate with the master.
virtual void setGenNtupleIntegerBranchValue(const TString &name, Int_t value)
Set the value of an integer branch in the gen tree.
virtual void setupBkgndVectors()=0
Method to set up the storage for background-related quantities called by setBkgndClassNames.
TVectorD * vectorPar_
Parameters to send to the slaves.
LauGenNtuple * genNtuple_
The generated ntuple.
File containing declaration of LauAbsFitModel class.
const TString & name() const
The parameter name.
virtual void setupGenNtupleBranches()=0
Setup the generation ntuple branches.
UInt_t addFitParameters(LauPdfList &pdfList)
Add parameters of the PDFs in the list to the list of all fit parameters.
virtual void generate(const TString &dataFileName, const TString &dataTreeName, const TString &histFileName, const TString &tableFileNameBase)
Generate toy MC.
virtual void setBkgndClassNames(const std::vector< TString > &names)
Setup the background class names.
Bool_t useDP() const
Is the Dalitz plot term in the likelihood.
TVectorD * vectorRes_
Parameters returned from the slaves.
Int_t buildIndex(const TString &majorName, const TString &minorName="0")
Create an index table using leaves of the tree.
void compareFitData(UInt_t toyMCScale=10, const TString &mcFileName="fitToyMC.root", const TString &tableFileName="fitToyMCTable.tex", Bool_t poissonSmearing=kTRUE)
Specify that a toy MC sample should be created for a successful fit to an experiment.
const TString & treeName() const
Ntuple tree name.
Definition: LauGenNtuple.hh:56
Int_t getIntegerBranchValue(const TString &name) const
Get value of an integer branch.
virtual void cacheInputSWeights()
Cache the value of the sWeights to be used in the sFit.
virtual void setGenNtupleDoubleBranchValue(const TString &name, Double_t value)
Set the value of a double branch in the gen tree.
void deleteAndRecreateTree()
Delete and recreate tree.
Double_t negError() const
The lower error on the parameter.
TSocket * sMaster_
A socket to enable parallel setup.
virtual const std::vector< LauParameter * > & getParameters() const
Retrieve the parameters of the PDF, e.g. so that they can be loaded into a fit.
Definition: LauAbsPdf.hh:321
void addFriendTree(const TString &rootFileName, const TString &rootTreeName)
Add a friend tree.
virtual void storePerEvtLlhds()=0
Store the per-event likelihood values.
virtual void printEventInfo(UInt_t iEvt) const
Prints the values of all the fit variables for the specified event - useful for diagnostics.
Bool_t pdfsDependOnDP() const
Do any of the PDFs have a dependence on the DP?
virtual ~LauAbsFitModel()
Destructor.
Bool_t compareFitData_
Option to make toy from 1st successful experiment.
File containing declaration of LauPrint class.
Double_t prodPdfValue(LauPdfList &pdfList, UInt_t iEvt)
Calculate the product of the per-event likelihoods of the PDFs in the list.
void runMaster(const TString &applicationCode, const TString &dataFileName="", const TString &dataTreeName="", const TString &histFileName="", const TString &tableFileName="", const UInt_t nSlaves=0)
Start the master process for the generation / fitting.
Double_t getDoubleBranchValue(const TString &name) const
Get value of a double branch.
std::vector< LauAbsPdf * > LauPdfList
List of Pdfs.
std::map< TString, Double_t > LauFitData
Type for holding event data.
void fillBranches()
Fill branches in the ntuple.
virtual Int_t getGenNtupleIntegerBranchValue(const TString &name) const
Get the value of an integer branch in the gen tree.
UInt_t nSlaves_
The number of slaves.
virtual LauSPlot::NumbMap fixdSpeciesNames() const =0
Returns the names and yields of species that are fixed in the fit.
LauBkgndClassMap bkgndClassNames_
The background class names.
TMessage * messageFromMaster_
Message from master to the slaves.
Class to define various output print commands.
Definition: LauPrint.hh:29
virtual void setupSPlotNtupleBranches()=0
Setup the branches of the sPlot tuple.
TMonitor * socketMonitor_
Parallel setup monitor.
virtual void fillGenNtupleBranches()
Fill the gen tuple branches.
Bool_t useAsymmFitErrors_
Option to use asymmetric errors.
virtual void setSPlotNtupleDoubleBranchValue(const TString &name, Double_t value)
Set the value of a double branch in the sPlot tree.
void logLikeFun(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
void printFitParameters(const LauPdfList &pdfList, std::ostream &fout) const
Print the fit parameters for all PDFs in the list.
void readExperimentData(UInt_t iExpt)
Read events only for the given experiment.
Double_t posError() const
The upper error on the parameter.
UInt_t slaveId_
Slave id number.
virtual Double_t getTotEvtLikelihood(UInt_t iEvt)=0
Calculates the likelihood for a given event.
TStopwatch cumulTimer_
The total fit timer.
File containing LauFitter namespace.
virtual void propagateParUpdates()=0
This function (specific to each model) calculates anything that depends on the fit parameter values...
void initSockets()
Initialise socket connections for the slaves.
Bool_t withinMINOS_
Flag to indicate if MINOS is currently running.
LauFitDataTree * inputFitData_
The input data.
Bool_t doSFit() const
Return the flag to store the status of using an sFit or not.
std::vector< LauParameter * > LauParameterPList
List of parameter pointers.
const TString & treeName() const
Retrieve the tree name.
LauParameterList extraVars_
Extra variables that aren&#39;t in the fit but are stored in the ntuple.
Bool_t firstStage() const
Check whether the parameter should be floated only in the first stage of a two stage fit...
void writeOutFitResults()
Write out fit results.
void addIntegerBranch(const TString &name)
Add integer branch to tree.
Definition: LauGenNtuple.cc:73
const TString & fileName() const
Ntuple file name.
Definition: LauGenNtuple.hh:50
Abstract interface to the fitting and toy MC model.
Double_t calculateLogLike(Int_t npar, Double_t *par, Int_t iflag)
Calculate the new value of the negative log likelihood.
UInt_t numberBadFits_
The number of bad fits.
virtual void addGenNtupleIntegerBranch(const TString &name)
Add a branch to the gen tree for storing an integer.
UInt_t nExpt() const
Obtain the number of experiments.
File containing declaration of LauParameter class.
Bool_t doPoissonSmearing() const
Determine whether Poisson smearing is enabled for the toy MC generation.
Predicate to allow counting of the number of fixed parameters.
virtual Double_t getGenNtupleDoubleBranchValue(const TString &name) const
Get the value of a double branch in the gen tree.
void setNExpts(UInt_t nExperiments, UInt_t firstExperiment=0)
Set the number of experiments and the first experiment.
virtual void weightEvents(const TString &dataFileName, const TString &dataTreeName)=0
Reweighting - allows e.g. MC events to be weighted by the DP model.
Bool_t doEMLFit() const
Determine whether an extended maximum likelihood fit it being performed.
Bool_t secondStage() const
Check whether the parameter should be floated only in the second stage of a two stage fit...
File containing declaration of LauComplex class.
UInt_t numberOKFits_
The number of good fits.
const TString & fileName() const
Retrieve the file name.
virtual Bool_t isDPDependent() const
Specifies whether or not the PDF is DP dependent.
Definition: LauAbsPdf.hh:110
virtual LauSPlot::NumbMap freeSpeciesNames() const =0
Returns the names and yields of species that are free in the fit.
Class to store the results from the toy MC generation into an ntuple.
Definition: LauGenNtuple.hh:32
File containing declaration of LauSPlot class.
virtual Bool_t splitSignal() const =0
Check if the signal is split into well-reconstructed and mis-reconstructed types. ...
Double_t error() const
The error on the parameter.
Double_t worstLogLike_
The worst LL value found so far.
TString sPlotTreeName_
The name of the sPlot tree.
virtual void initialiseDPModels()=0
Initialise the DP models.
LauGenNtuple * sPlotNtuple_
The sPlot ntuple.
TString sPlotVerbosity_
Control the verbosity of the sFit.
const TString nullString_
An empty string.
void clearFitParVectors()
Clear the vectors containing fit parameters.
virtual void addSPlotNtupleDoubleBranch(const TString &name)
Add a branch to the sPlot tree for storing a double.
void clearExtraVarVectors()
Clear the vectors containing extra ntuple variables.
virtual void cacheInputFitVars()=0
Cache the input data values to calculate the likelihood during the fit.
void fitExpt()
Routine to perform the actual fit for a given experiment.
Double_t sWeightScaleFactor_
The sWeight scaling factor.
File containing declaration of LauAbsPdf class.
UInt_t bkgndClassID(const TString &className) const
The number assigned to a background class.
const LauFitData & getData(UInt_t iEvt) const
Retrieve the data for a given event.
virtual void writeOutAllFitResults()
Write out any fit results.
virtual void finaliseFitResults(const TString &tableFileName)=0
Write the results of the fit into the ntuple.
void fitSlave(const TString &dataFileName, const TString &dataTreeName)
Slaves required when performing a simultaneous fit.
UInt_t firstExpt_
The number of the first experiment to consider.
void setIntegerBranchValue(const TString &name, Int_t value)
Set value of an integer branch.
Definition: LauGenNtuple.cc:91
virtual void checkInitFitParams()=0
Update initial fit parameters if required.
virtual void printVarsInfo() const
Same as printEventInfo, but printing out the values of the variables in the fit.
Bool_t writeSPlotData_
Option to write sPlot data.
Bool_t validBkgndClass(const TString &className) const
Check if the given background class is in the list.
Bool_t fitToyMCPoissonSmear_
Option to perform Poisson smearing.
TString sWeightBranchName_
The name of the sWeight branch.
UInt_t iExpt_
The number of the current experiment.
TMessage * messageFromSlave_
Message from slaves to the master.
LauParameterPList fitVars_
Internal vectors of fit parameters.
virtual Bool_t scfDPSmear() const =0
Check if the mis-reconstructed signal is to be smeared in the DP.
Bool_t cacheFitData(const TString &dataFileName, const TString &dataTreeName)
Store variables from the input file into the internal data storage.
void runSlave(const TString &applicationCode, const TString &dataFileName="", const TString &dataTreeName="", const TString &histFileName="", const TString &tableFileName="", const TString &addressMaster="")
Start the slave process for the generation / fitting.
void fit(const TString &dataFileName, const TString &dataTreeName, const TString &histFileName, const TString &tableFileNameBase)
Perform the total fit.
virtual void setSPlotNtupleIntegerBranchValue(const TString &name, Int_t value)
Set the value of an integer branch in the sPlot tree.
void fitMaster(const TString &dataFileName, const TString &dataTreeName, const TString &histFileName, const TString &tableFileNameBase)
Master fit used when simultaneous fit is required.
UInt_t fitToyMCScale_
The scaling factor (toy vs data statistics)
Bool_t runMinimisation()
Routine to perform the minimisation.
virtual void fillSPlotNtupleBranches()
Fill the sPlot tuple.
TString fitToyMCFileName_
The output file name for Toy MC.
void writeOutGenResults()
Write out the results from the generation.
void printFormat(std::ostream &stream, Double_t value) const
Method to choose the printing format to a specified level of precision.
Definition: LauPrint.cc:32
virtual LauSPlot::TwoDMap twodimPDFs() const =0
Returns the species and variables for all 2D PDFs in the fit.
Bool_t haveBranch(const TString &name) const
Check if the named branch is stored.
const TString & bkgndClassName(UInt_t classID) const
Get the name of a background class from the number.
Bool_t enableEmbedding() const
Determine whether embedding of events is enabled in the generation.
Double_t NLL_
The negative log-likelihood.
Class for defining the abstract interface for PDF classes.
Definition: LauAbsPdf.hh:40
Double_t getLogLikelihood(UInt_t iStart, UInt_t iEnd)
Calculate the sum of the log-likelihood over the specified events.
UInt_t nExpt_
The number of experiments to consider.
LauFitNtuple * fitNtuple_
The fit ntuple.
virtual LauSPlot::NameSet variableNames() const =0
Returns the names of all variables in the fit.
Double_t value() const
The value of the parameter.
UInt_t nEvents() const
Retrieve the number of events.
virtual void addGenNtupleDoubleBranch(const TString &name)
Add a branch to the gen tree for storing a double.
friend class LauFitNtuple
LauFitNtuple is a friend class.
void updateFitParameters(LauPdfList &pdfList)
Update the fit parameters for all PDFs in the list.
void addDoubleBranch(const TString &name)
Add double branch to tree.
Definition: LauGenNtuple.cc:82
virtual void setParsFromMinuit(Double_t *par, Int_t npar)
This function sets the parameter values from Minuit.
Int_t fitStatus_
The status of the fit.
virtual void initialise()=0
Initialise the fit par vectors.
Class to store the input fit variables.
File containing declaration of LauFitNtuple class.
virtual void calculateSPlotData()
Calculate the sPlot data.
File containing declaration of LauGenNtuple class.
void setDoubleBranchValue(const TString &name, Double_t value)
Set value of a double branch.
Definition: LauGenNtuple.cc:96
std::vector< Double_t > sWeights_
The vector of sWeights.
UInt_t nFreeParams_
The number of free fit parameters.
TVirtualFitter * fitter(const TString &fitterString="Minuit", Int_t maxPar=100)
Method that provides a pointer to a TVirtualFitter object.
Definition: LauFitter.cc:20