laura is hosted by Hepforge, IPPP Durham
Laura++  v2r2
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 - 2014.
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 <limits>
17 #include <vector>
18 
19 #include "TMessage.h"
20 #include "TMonitor.h"
21 #include "TServerSocket.h"
22 #include "TSocket.h"
23 #include "TSystem.h"
24 #include "TVectorD.h"
25 #include "TVirtualFitter.h"
26 
27 #include "LauAbsFitModel.hh"
28 #include "LauAbsFitter.hh"
29 #include "LauAbsPdf.hh"
30 #include "LauComplex.hh"
31 #include "LauFitter.hh"
32 #include "LauFitDataTree.hh"
33 #include "LauFitNtuple.hh"
34 #include "LauGenNtuple.hh"
35 #include "LauParameter.hh"
36 #include "LauParamFixed.hh"
37 #include "LauPrint.hh"
38 #include "LauSPlot.hh"
39 
41 
42 
44  storeCon_(0),
45  twoStageFit_(kFALSE),
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  usingDP_(kTRUE),
56  pdfsDependOnDP_(kFALSE),
57  firstExpt_(0),
58  nExpt_(0),
59  evtsPerExpt_(0),
60  iExpt_(0),
61  inputFitData_(0),
62  fitNtuple_(0),
63  genNtuple_(0),
64  sPlotNtuple_(0),
65  fitStatus_(0),
66  NLL_(0),
67  numberOKFits_(0),
68  numberBadFits_(0),
69  nParams_(0),
70  nFreeParams_(0),
71  worstLogLike_(std::numeric_limits<Double_t>::max()),
72  withinAsymErrorCalc_(kFALSE),
73  nullString_(""),
74  doSFit_(kFALSE),
75  sWeightBranchName_(""),
76  sWeightScaleFactor_(1.0),
77  fitToyMCFileName_("fitToyMC.root"),
78  fitToyMCTableName_("fitToyMCTable"),
79  fitToyMCScale_(10),
80  fitToyMCPoissonSmear_(kFALSE),
81  sPlotFileName_(""),
82  sPlotTreeName_(""),
83  sPlotVerbosity_(""),
84  sMaster_(0),
85  messageFromMaster_(0),
86  slaveId_(-1),
87  nSlaves_(0),
88  parValues_(0)
89 {
90 }
91 
92 
94 {
95  delete inputFitData_; inputFitData_ = 0;
96  delete fitNtuple_; fitNtuple_ = 0;
97  delete genNtuple_; genNtuple_ = 0;
98  delete sPlotNtuple_; sPlotNtuple_ = 0;
99  delete sMaster_; sMaster_ = 0;
100  delete[] parValues_; parValues_ = 0;
101 
102  // Remove the components created to apply constraints to fit parameters
103  for (std::vector<LauAbsRValue*>::iterator iter = conVars_.begin(); iter != conVars_.end(); ++iter){
104  if ( !(*iter)->isLValue() ){
105  delete (*iter);
106  (*iter) = 0;
107  }
108  }
109 }
110 
111 void LauAbsFitModel::run(const TString& applicationCode, const TString& dataFileName, const TString& dataTreeName,
112  const TString& histFileName, const TString& tableFileName)
113 {
114  // Chose whether you want to generate or fit events in the Dalitz plot.
115  // To generate events choose applicationCode = "gen", to fit events choose
116  // applicationCode = "fit".
117 
118  TString runCode(applicationCode);
119  runCode.ToLower();
120 
121  TString histFileNameCopy(histFileName);
122  TString tableFileNameCopy(tableFileName);
123  TString dataFileNameCopy(dataFileName);
124  TString dataTreeNameCopy(dataTreeName);
125 
126  // Initialise the fit par vectors. Each class that inherits from this one
127  // must implement this sensibly for all vectors specified in clearFitParVectors,
128  // i.e. specify parameter names, initial, min, max and fixed values
129  this->initialise();
130 
131  // Add variables to Gaussian constrain to a list
132  this->addConParameters();
133 
134  if (dataFileNameCopy == "") {dataFileNameCopy = "data.root";}
135  if (dataTreeNameCopy == "") {dataTreeNameCopy = "genResults";}
136 
137  if (runCode.Contains("gen")) {
138 
139  if (histFileNameCopy == "") {histFileNameCopy = "parInfo.root";}
140  if (tableFileNameCopy == "") {tableFileNameCopy = "genResults";}
141 
142  this->setGenValues();
143  this->generate(dataFileNameCopy, dataTreeNameCopy, histFileNameCopy, tableFileNameCopy);
144 
145  } else if (runCode.Contains("fit")) {
146 
147  if (histFileNameCopy == "") {histFileNameCopy = "parInfo.root";}
148  if (tableFileNameCopy == "") {tableFileNameCopy = "fitResults";}
149 
150  this->fit(dataFileNameCopy, dataTreeNameCopy, histFileNameCopy, tableFileNameCopy);
151 
152  } else if (runCode.Contains("reweight")) {
153 
154  this->weightEvents(dataFileNameCopy, dataTreeNameCopy);
155  }
156 }
157 
158 void LauAbsFitModel::runSlave(const TString& dataFileName, const TString& dataTreeName,
159  const TString& histFileName, const TString& tableFileName,
160  const TString& addressMaster, const UInt_t portMaster)
161 {
162  if ( sMaster_ != 0 ) {
163  std::cerr << "ERROR in LauAbsFitModel::runSlave : master socket already present" << std::endl;
164  return;
165  }
166 
167  // Open connection to master
168  sMaster_ = new TSocket(addressMaster, portMaster);
169  sMaster_->Recv( messageFromMaster_ );
170  messageFromMaster_->ReadUInt( slaveId_ );
171  messageFromMaster_->ReadUInt( nSlaves_ );
172 
173  delete messageFromMaster_;
174  messageFromMaster_ = 0;
175 
176  std::cout << "INFO in LauAbsFitModel::runSlave : Established connection to master on port " << portMaster << std::endl;
177  std::cout << " : We are slave " << slaveId_ << " of " << nSlaves_ << std::endl;
178 
179  // Initialise the fit par vectors. Each class that inherits from this one
180  // must implement this sensibly for all vectors specified in clearFitParVectors,
181  // i.e. specify parameter names, initial, min, max and fixed values
182  this->initialise();
183 
184  // NB call to addConParameters() is intentionally not included here cf.
185  // run() since this has to be dealt with by the master to avoid
186  // multiple inclusions of each penalty term
187 
188  // Create array to efficiently exchange parameter values with master
189  nParams_ = fitVars_.size();
190  parValues_ = new Double_t[nParams_];
191  for ( UInt_t iPar(0); iPar < nParams_; ++iPar ) {
192  parValues_[iPar] = fitVars_[iPar]->initValue();
193  }
194 
195  TString dataFileNameCopy(dataFileName);
196  TString dataTreeNameCopy(dataTreeName);
197  TString histFileNameCopy(histFileName);
198  TString tableFileNameCopy(tableFileName);
199 
200  if (dataFileNameCopy == "") {dataFileNameCopy = "data.root";}
201  if (dataTreeNameCopy == "") {dataTreeNameCopy = "genResults";}
202  if (histFileNameCopy == "") {histFileNameCopy = "parInfo.root";}
203  if (tableFileNameCopy == "") {tableFileNameCopy = "fitResults";}
204 
205  this->fitSlave(dataFileNameCopy, dataTreeNameCopy, histFileNameCopy, tableFileNameCopy);
206 
207  std::cout << "INFO in LauAbsFitModel::runSlave : Fit slave " << slaveId_ << " has finished successfully" << std::endl;
208 }
209 
210 void LauAbsFitModel::doSFit( const TString& sWeightBranchName, Double_t scaleFactor )
211 {
212  if ( sWeightBranchName == "" ) {
213  std::cerr << "WARNING in LauAbsFitModel::doSFit : sWeight branch name is empty string, not setting-up sFit." << std::endl;
214  return;
215  }
216 
217  doSFit_ = kTRUE;
218  sWeightBranchName_ = sWeightBranchName;
219  sWeightScaleFactor_ = scaleFactor;
220 }
221 
222 void LauAbsFitModel::setBkgndClassNames( const std::vector<TString>& names )
223 {
224  if ( !bkgndClassNames_.empty() ) {
225  std::cerr << "WARNING in LauAbsFitModel::setBkgndClassNames : Names already stored, not changing them." << std::endl;
226  return;
227  }
228 
229  UInt_t nBkgnds = names.size();
230  for ( UInt_t i(0); i < nBkgnds; ++i ) {
231  bkgndClassNames_.insert( std::make_pair( i, names[i] ) );
232  }
233 
234  this->setupBkgndVectors();
235 }
236 
237 Bool_t LauAbsFitModel::validBkgndClass( const TString& className ) const
238 {
239  if ( bkgndClassNames_.empty() ) {
240  return kFALSE;
241  }
242 
243  Bool_t found(kFALSE);
244  for ( LauBkgndClassMap::const_iterator iter = bkgndClassNames_.begin(); iter != bkgndClassNames_.end(); ++iter ) {
245  if ( iter->second == className ) {
246  found = kTRUE;
247  break;
248  }
249  }
250 
251  return found;
252 }
253 
254 UInt_t LauAbsFitModel::bkgndClassID( const TString& className ) const
255 {
256  if ( ! this->validBkgndClass( className ) ) {
257  std::cerr << "ERROR in LauAbsFitModel::bkgndClassID : Request for ID for invalid background class \"" << className << "\"." << std::endl;
258  return (bkgndClassNames_.size() + 1);
259  }
260 
261  UInt_t bgID(0);
262  for ( LauBkgndClassMap::const_iterator iter = bkgndClassNames_.begin(); iter != bkgndClassNames_.end(); ++iter ) {
263  if ( iter->second == className ) {
264  bgID = iter->first;
265  break;
266  }
267  }
268 
269  return bgID;
270 }
271 
272 const TString& LauAbsFitModel::bkgndClassName( UInt_t classID ) const
273 {
274  LauBkgndClassMap::const_iterator iter = bkgndClassNames_.find( classID );
275 
276  if ( iter == bkgndClassNames_.end() ) {
277  std::cerr << "ERROR in LauAbsFitModel::bkgndClassName : Request for name of invalid background class ID " << classID << "." << std::endl;
278  return nullString_;
279  }
280 
281  return iter->second;
282 }
283 
285 {
286  std::cout << "INFO in LauAbsFitModel::clearFitParVectors : Clearing fit variable vectors" << std::endl;
287 
288  // Remove the components created to apply constraints to fit parameters
289  for (std::vector<LauAbsRValue*>::iterator iter = conVars_.begin(); iter != conVars_.end(); ++iter){
290  if ( !(*iter)->isLValue() ){
291  delete (*iter);
292  (*iter) = 0;
293  }
294  }
295  conVars_.clear();
296  fitVars_.clear();
297 }
298 
300 {
301  std::cout << "INFO in LauAbsFitModel::clearExtraVarVectors : Clearing extra ntuple variable vectors" << std::endl;
302  extraVars_.clear();
303 }
304 
306 {
307  // makes sure each parameter holds its genValue as its current value
308  for (LauParameterPList::iterator iter = fitVars_.begin(); iter != fitVars_.end(); ++iter) {
309  (*iter)->value((*iter)->genValue());
310  }
311  this->propagateParUpdates();
312 }
313 
314 void LauAbsFitModel::writeSPlotData(const TString& fileName, const TString& treeName, Bool_t storeDPEfficiency, const TString& verbosity)
315 {
316  if (this->writeSPlotData()) {
317  std::cerr << "ERROR in LauAbsFitModel::writeSPlotData : Already have an sPlot ntuple setup, not doing it again." << std::endl;
318  return;
319  }
320  writeSPlotData_ = kTRUE;
321  sPlotFileName_ = fileName;
322  sPlotTreeName_ = treeName;
323  sPlotVerbosity_ = verbosity;
324  storeDPEff_ = storeDPEfficiency;
325 }
326 
327 // TODO : histFileName isn't used here at the moment but could be used for
328 // storing the values of the parameters used in the generation.
329 // These could then be read and used for setting the "true" values
330 // in a subsequent fit.
331 void LauAbsFitModel::generate(const TString& dataFileName, const TString& dataTreeName, const TString& /*histFileName*/, const TString& tableFileNameBase)
332 {
333  // Create the ntuple for storing the results
334  std::cout << "INFO in LauAbsFitModel::generate : Creating generation ntuple." << std::endl;
335  if (genNtuple_ != 0) {delete genNtuple_; genNtuple_ = 0;}
336  genNtuple_ = new LauGenNtuple(dataFileName,dataTreeName);
337 
338  // add branches for storing the experiment number and the number of
339  // the event within the current experiment
340  this->addGenNtupleIntegerBranch("iExpt");
341  this->addGenNtupleIntegerBranch("iEvtWithinExpt");
342  this->setupGenNtupleBranches();
343 
344  // Start the cumulative timer
345  cumulTimer_.Start();
346 
347  Bool_t genOK(kTRUE);
348  do {
349  // Loop over the number of experiments
350  for (iExpt_ = firstExpt_; iExpt_ < (firstExpt_+nExpt_); ++iExpt_) {
351 
352  // Start the timer to see how long each experiment takes to generate
353  timer_.Start();
354 
355  // Store the experiment number in the ntuple
356  this->setGenNtupleIntegerBranchValue("iExpt",iExpt_);
357 
358  // Do the generation for this experiment
359  std::cout << "INFO in LauAbsFitModel::generate : Generating experiment number " << iExpt_ << std::endl;
360  genOK = this->genExpt();
361 
362  // Stop the timer and see how long the program took so far
363  timer_.Stop();
364  timer_.Print();
365 
366  if (!genOK) {
367  // delete and recreate an empty tree
369 
370  // then break out of the experiment loop
371  std::cerr << "ERROR in LauAbsFitModel::generate : Problem in toy MC generation. Starting again with updated parameters..." << std::endl;
372  break;
373  }
374 
375  if (this->writeLatexTable()) {
376  TString tableFileName(tableFileNameBase);
377  tableFileName += "_";
378  tableFileName += iExpt_;
379  tableFileName += ".tex";
380  this->writeOutTable(tableFileName);
381  }
382 
383  } // Loop over number of experiments
384  } while (!genOK);
385 
386  // Print out total timing info.
387  cumulTimer_.Stop();
388  std::cout << "INFO in LauAbsFitModel::generate : Finished generating all experiments." << std::endl;
389  std::cout << "INFO in LauAbsFitModel::generate : Cumulative timing:" << std::endl;
390  cumulTimer_.Print();
391 
392  // Build the event index
393  std::cout << "INFO in LauAbsFitModel::generate : Building experiment:event index." << std::endl;
394  // TODO - can test this return value?
395  //Int_t nIndexEntries =
396  genNtuple_->buildIndex("iExpt","iEvtWithinExpt");
397 
398  // Write out toy MC ntuple
399  std::cout << "INFO in LauAbsFitModel::generate : Writing data to file " << dataFileName << "." << std::endl;
401 }
402 
404 {
406 }
407 
409 {
411 }
412 
414 {
415  genNtuple_->setIntegerBranchValue(name,value);
416 }
417 
419 {
420  genNtuple_->setDoubleBranchValue(name,value);
421 }
422 
424 {
425  return genNtuple_->getIntegerBranchValue(name);
426 }
427 
429 {
430  return genNtuple_->getDoubleBranchValue(name);
431 }
432 
434 {
436 }
437 
439 {
441 }
442 
444 {
446 }
447 
449 {
450  sPlotNtuple_->setIntegerBranchValue(name,value);
451 }
452 
454 {
455  sPlotNtuple_->setDoubleBranchValue(name,value);
456 }
457 
459 {
461 }
462 
463 void LauAbsFitModel::fit(const TString& dataFileName, const TString& dataTreeName, const TString& histFileName, const TString& tableFileNameBase)
464 {
465  // Routine to perform the total fit.
466 
467  std::cout << "INFO in LauAbsFitModel::fit : First experiment = " << firstExpt_ << std::endl;
468  std::cout << "INFO in LauAbsFitModel::fit : Number of experiments = " << nExpt_ << std::endl;
469 
470  // Start the cumulative timer
471  cumulTimer_.Start();
472 
473  numberOKFits_ = 0, numberBadFits_ = 0;
474  fitStatus_ = -1;
475 
476  // Create and setup the fit results ntuple
477  std::cout << "INFO in LauAbsFitModel::fit : Creating fit ntuple." << std::endl;
478  if (fitNtuple_ != 0) {delete fitNtuple_; fitNtuple_ = 0;}
479  fitNtuple_ = new LauFitNtuple(histFileName, this->useAsymmFitErrors());
480 
481  // Create and setup the sPlot ntuple
482  if (this->writeSPlotData()) {
483  std::cout << "INFO in LauAbsFitModel::fit : Creating sPlot ntuple." << std::endl;
484  if (sPlotNtuple_ != 0) {delete sPlotNtuple_; sPlotNtuple_ = 0;}
486  this->setupSPlotNtupleBranches();
487  }
488 
489  // This reads in the given dataFile and creates an input
490  // fit data tree that stores them for all events and experiments.
491  Bool_t dataOK = this->cacheFitData(dataFileName,dataTreeName);
492  if (!dataOK) {
493  std::cerr << "ERROR in LauAbsFitModel::fit : Problem caching the fit data." << std::endl;
494  gSystem->Exit(EXIT_FAILURE);
495  }
496 
497  // Loop over the number of experiments
498  for (iExpt_ = firstExpt_; iExpt_ < (firstExpt_+nExpt_); ++iExpt_) {
499 
500  // Start the timer to see how long each fit takes
501  timer_.Start();
502 
505 
506  if (this->eventsPerExpt() < 1) {
507  std::cerr << "ERROR in LauAbsFitModel::fit : Zero events in experiment " << iExpt_ << ", skipping..." << std::endl;
508  timer_.Stop();
509  continue;
510  }
511 
512  // Now the sub-classes must implement whatever they need to do
513  // to cache any more input fit data they need in order to
514  // calculate the likelihoods during the fit.
515  // They need to use the inputFitData_ tree as input. For example,
516  // inputFitData_ contains m13Sq and m23Sq. The appropriate fit model
517  // then caches the resonance dynamics for the signal model, as
518  // well as the background likelihood values in the Dalitz plot
519  this->cacheInputFitVars();
520  if ( this->doSFit() ) {
521  this->cacheInputSWeights();
522  }
523 
524  // Do the fit for this experiment
525  this->fitExpt();
526 
527  // Write the results into the ntuple
528  this->finaliseFitResults(tableFileNameBase);
529 
530  // Stop the timer and see how long the program took so far
531  timer_.Stop();
532  timer_.Print();
533 
534  // Store the per-event likelihood values
535  if ( this->writeSPlotData() ) {
536  this->storePerEvtLlhds();
537  }
538 
539  // Create a toy MC sample using the fitted parameters so that
540  // the user can compare the fit to the data.
541  if (compareFitData_ == kTRUE && fitStatus_ == 3) {
543  }
544 
545  // Keep track of how many fits worked or failed
546  // NB values of fitStatus_ now indicate the status of the error matrix:
547  // 0= not calculated at all
548  // 1= approximation only, not accurate
549  // 2= full matrix, but forced positive-definite
550  // 3= full accurate covariance matrix
551  if (fitStatus_ == 3) {
552  numberOKFits_++;
553  } else {
554  numberBadFits_++;
555  }
556 
557  } // Loop over number of experiments
558 
559  // Print out total timing info.
560  cumulTimer_.Stop();
561  std::cout << "INFO in LauAbsFitModel::fit : Cumulative timing:" << std::endl;
562  cumulTimer_.Print();
563 
564  // Print out stats on OK fits.
565  std::cout << "INFO in LauAbsFitModel::fit : Number of OK Fits = " << numberOKFits_ << std::endl;
566  std::cout << "INFO in LauAbsFitModel::fit : Number of Failed Fits = " << numberBadFits_ << std::endl;
567  Double_t fitEff(0.0);
568  if (nExpt_ != 0) {fitEff = numberOKFits_/(1.0*nExpt_);}
569  std::cout << "INFO in LauAbsFitModel::fit : Fit efficiency = " << fitEff*100.0 << "%." << std::endl;
570 
571  // Write out any fit results (ntuples etc...).
572  this->writeOutAllFitResults();
573  if ( this->writeSPlotData() ) {
574  this->calculateSPlotData();
575  }
576 }
577 
578 void LauAbsFitModel::fitSlave(const TString& dataFileName, const TString& dataTreeName, const TString& histFileName, const TString& tableFileNameBase)
579 {
580  // Create and setup the fit results ntuple
581  std::cout << "INFO in LauAbsFitModel::fitSlave : Creating fit ntuple." << std::endl;
582  if (fitNtuple_ != 0) {delete fitNtuple_; fitNtuple_ = 0;}
583  fitNtuple_ = new LauFitNtuple(histFileName, this->useAsymmFitErrors());
584 
585  // This reads in the given dataFile and creates an input
586  // fit data tree that stores them for all events and experiments.
587  Bool_t dataOK = this->cacheFitData(dataFileName,dataTreeName);
588  if (!dataOK) {
589  std::cerr << "ERROR in LauAbsFitModel::fitSlave : Problem caching the fit data." << std::endl;
590  gSystem->Exit(EXIT_FAILURE);
591  }
592 
593  // Now process the various requests from the master
594 
595  TMessage messageToMaster(kMESS_ANY);
596 
597  while ( kTRUE ) {
598 
599  sMaster_->Recv( messageFromMaster_ );
600 
601  if ( messageFromMaster_->What() == kMESS_STRING ) {
602 
603  TString msgStr;
604  messageFromMaster_->ReadTString( msgStr );
605 
606  std::cout << "INFO in LauAbsFitModel::fitSlave : Received message from master: " << msgStr << std::endl;
607 
608  if ( msgStr == "Send Parameters" ) {
609 
610  // Update initial fit parameters if required (e.g. if using random numbers).
611  this->checkInitFitParams();
612 
613  // Send the fit parameters
614  TObjArray array;
615  for ( LauParameterPList::iterator iter = fitVars_.begin(); iter != fitVars_.end(); ++iter ) {
616  array.Add( *iter );
617  }
618 
619  messageToMaster.Reset( kMESS_OBJECT );
620  messageToMaster.WriteObject( &array );
621  sMaster_->Send( messageToMaster );
622 
623  } else if ( msgStr == "Read Expt" ) {
624 
625  // Read the data for this experiment
626  messageFromMaster_->ReadUInt( iExpt_ );
627 
629  UInt_t nEvent = inputFitData_->nEvents();
630  this->eventsPerExpt( nEvent );
631 
632  if ( nEvent < 1 ) {
633  std::cerr << "ERROR in LauAbsFitModel::fitSlave : Zero events in experiment " << firstExpt_ << ", the master should skip this experiment..." << std::endl;
634  }
635 
636  messageToMaster.Reset( kMESS_ANY );
637  messageToMaster.WriteUInt( slaveId_ );
638  messageToMaster.WriteUInt( nEvent );
639  sMaster_->Send( messageToMaster );
640 
641  } else if ( msgStr == "Cache" ) {
642 
643  // Perform the caching
644 
645  this->cacheInputFitVars();
646 
647  messageToMaster.Reset( kMESS_ANY );
648  messageToMaster.WriteUInt( slaveId_ );
649  messageToMaster.WriteBool( kTRUE );
650  sMaster_->Send( messageToMaster );
651 
652  } else if ( msgStr == "Write Results" ) {
653 
654  this->writeOutAllFitResults();
655 
656  messageToMaster.Reset( kMESS_ANY );
657  messageToMaster.WriteUInt( slaveId_ );
658  messageToMaster.WriteBool( kTRUE );
659  sMaster_->Send( messageToMaster );
660 
661  } else if ( msgStr == "Finish" ) {
662 
663  std::cout << "INFO in LauAbsFitModel::fitSlave : Message from master to finish" << std::endl;
664  break;
665  } else {
666  std::cerr << "ERROR in LauAbsFitModel::fitSlave : Unexpected message from master" << std::endl;
667  gSystem->Exit( EXIT_FAILURE );
668  }
669 
670  } else if ( messageFromMaster_->What() == kMESS_OBJECT ) {
671 
672  std::cout << "INFO in LauAbsFitModel::fitSlave : Received message from master: Finalise" << std::endl;
673 
674  messageFromMaster_->ReadInt( fitStatus_ );
675  messageFromMaster_->ReadDouble( NLL_ );
676 
677  TObjArray * objarray = dynamic_cast<TObjArray*>( messageFromMaster_->ReadObject( messageFromMaster_->GetClass() ) );
678  if ( ! objarray ) {
679  std::cerr << "ERROR in LauAbsFitModel::fitSlave : Error reading parameters from master" << std::endl;
680  gSystem->Exit( EXIT_FAILURE );
681  }
682 
683  TMatrixD * covMat = dynamic_cast<TMatrixD*>( messageFromMaster_->ReadObject( messageFromMaster_->GetClass() ) );
684  if ( ! covMat ) {
685  std::cerr << "ERROR in LauAbsFitModel::fitSlave : Error reading covariance matrix from master" << std::endl;
686  gSystem->Exit( EXIT_FAILURE );
687  }
688  covMatrix_.Clear();
689  covMatrix_.ResizeTo( covMat->GetNrows(), covMat->GetNcols() );
690  covMatrix_.SetMatrixArray( covMat->GetMatrixArray() );
691  delete covMat; covMat = 0;
692 
693 
694  UInt_t nPars = objarray->GetEntries();
695 
696  if ( nPars != nParams_ ) {
697  std::cerr << "ERROR in LauAbsFitModel::fitSlave : Unexpected number of parameters received from master" << std::endl;
698  std::cerr << " ::fitSlave : Received " << nPars << " when expecting " << nParams_ << std::endl;
699  gSystem->Exit( EXIT_FAILURE );
700  }
701 
702  for ( UInt_t iPar(0); iPar < nPars; ++iPar ) {
703  LauParameter* parameter = dynamic_cast<LauParameter*>( (*objarray)[iPar] );
704  if ( ! parameter ) {
705  std::cerr << "ERROR in LauAbsFitModel::fitSlave : Error reading parameter from master" << std::endl;
706  gSystem->Exit( EXIT_FAILURE );
707  }
708 
709  if ( parameter->name() != fitVars_[iPar]->name() ) {
710  std::cerr << "ERROR in LauAbsFitModel::fitSlave : Error reading parameter from master" << std::endl;
711  gSystem->Exit( EXIT_FAILURE );
712  }
713 
714  *(fitVars_[iPar]) = *parameter;
715  }
716 
717  // Write the results into the ntuple
718  this->finaliseFitResults( tableFileNameBase );
719 
720  // Store the per-event likelihood values
721  if ( this->writeSPlotData() ) {
722  this->storePerEvtLlhds();
723  }
724 
725  // Create a toy MC sample using the fitted parameters so that
726  // the user can compare the fit to the data.
727  if (compareFitData_ == kTRUE && fitStatus_ == 3) {
729  }
730 
731  // Send the finalised parameters
732  TObjArray array;
733  for ( LauParameterPList::iterator iter = fitVars_.begin(); iter != fitVars_.end(); ++iter ) {
734  array.Add( *iter );
735  }
736 
737  messageToMaster.Reset( kMESS_ANY );
738  messageToMaster.WriteUInt( slaveId_ );
739  messageToMaster.WriteBool( kTRUE );
740  messageToMaster.WriteObject( &array );
741  sMaster_->Send( messageToMaster );
742 
743  } else if ( messageFromMaster_->What() == kMESS_ANY ) {
744 
745  UInt_t nPars(0);
746  messageFromMaster_->ReadUInt( nPars );
747 
748  if ( nPars != nParams_ ) {
749  std::cerr << "ERROR in LauAbsFitModel::fitSlave : Unexpected number of parameters received from master" << std::endl;
750  std::cerr << " ::fitSlave : Received " << nPars << " when expecting " << nParams_ << std::endl;
751  gSystem->Exit( EXIT_FAILURE );
752  }
753 
754  messageFromMaster_->ReadFastArray( parValues_, nPars );
755 
756  for ( UInt_t iPar(0); iPar < nPars; ++iPar ) {
757  if ( ! fitVars_[iPar]->fixed() ) {
758  fitVars_[iPar]->value( parValues_[iPar] );
759  }
760  }
761  this->propagateParUpdates();
762 
763  Double_t negLogLike = this->getTotNegLogLikelihood();
764 
765  messageToMaster.Reset( kMESS_ANY );
766  messageToMaster.WriteDouble( negLogLike );
767  sMaster_->Send( messageToMaster );
768 
769  } else {
770  std::cerr << "ERROR in LauAbsFitModel::fitSlave : Unexpected message type" << std::endl;
771  gSystem->Exit( EXIT_FAILURE );
772  }
773 
774  delete messageFromMaster_;
775  messageFromMaster_ = 0;
776  }
777 }
778 
779 Bool_t LauAbsFitModel::cacheFitData(const TString& dataFileName, const TString& dataTreeName)
780 {
781  // From the input data stream, store the variables into the
782  // internal tree inputFitData_ that can be used by the sub-classes
783  // in calculating their likelihood functions for the fit
784  delete inputFitData_;
785  inputFitData_ = new LauFitDataTree(dataFileName,dataTreeName);
786  Bool_t dataOK = inputFitData_->findBranches();
787 
788  if (!dataOK) {
789  delete inputFitData_; inputFitData_ = 0;
790  }
791 
792  return dataOK;
793 }
794 
796 {
797  Bool_t hasBranch = inputFitData_->haveBranch( sWeightBranchName_ );
798  if ( ! hasBranch ) {
799  std::cerr << "ERROR in LauAbsFitModel::cacheInputSWeights : Input data does not contain variable \"" << sWeightBranchName_ << "\".\n";
800  std::cerr << " : Turning off sFit!" << std::endl;
801  doSFit_ = kFALSE;
802  return;
803  }
804 
805  UInt_t nEvents = this->eventsPerExpt();
806  sWeights_.clear();
807  sWeights_.reserve( nEvents );
808 
809  for (UInt_t iEvt = 0; iEvt < nEvents; ++iEvt) {
810 
811  const LauFitData& dataValues = inputFitData_->getData(iEvt);
812 
813  LauFitData::const_iterator iter = dataValues.find( sWeightBranchName_ );
814 
815  sWeights_.push_back( iter->second * sWeightScaleFactor_ );
816  }
817 }
818 
820 {
821  // Routine to perform the actual fit for the given experiment
822 
823  // Reset the worst likelihood found to its catch-all value
824  worstLogLike_ = std::numeric_limits<Double_t>::max();
825 
826  // Update initial fit parameters if required (e.g. if using random numbers).
827  this->checkInitFitParams();
828 
829  // Initialise the fitter
833 
836 
837  // Now ready for minimisation step
838  std::cout << "\nINFO in LauAbsFitModel::fitExpt : Start minimisation...\n";
839  std::pair<Int_t,Double_t> fitResult = LauFitter::fitter()->minimise();
840 
841  fitStatus_ = fitResult.first;
842  NLL_ = fitResult.second;
843 
844  // If we're doing a two stage fit we can now release (i.e. float)
845  // the 2nd stage parameters and re-fit
846  if (this->twoStageFit()) {
847 
848  if ( fitStatus_ != 3 ) {
849  std::cerr << "ERROR in LauAbsFitModel:fitExpt : Not running second stage fit since first stage failed." << std::endl;
851  } else {
856  fitResult = LauFitter::fitter()->minimise();
857  }
858  }
859 
860  fitStatus_ = fitResult.first;
861  NLL_ = fitResult.second;
862  const TMatrixD& covMat = LauFitter::fitter()->covarianceMatrix();
863  covMatrix_.Clear();
864  covMatrix_.ResizeTo( covMat.GetNrows(), covMat.GetNcols() );
865  covMatrix_.SetMatrixArray( covMat.GetMatrixArray() );
866 
867  // Store the final fit results and errors into protected internal vectors that
868  // all sub-classes can use within their own finalFitResults implementation
869  // used below (e.g. putting them into an ntuple in a root file)
872 }
873 
875 {
876  // Write out histograms at end
877  if (fitNtuple_ != 0) {
879  }
880 }
881 
883 {
884  if (sPlotNtuple_ != 0) {
887  LauSPlot splot(sPlotNtuple_->fileName(), sPlotNtuple_->treeName(), this->firstExpt(), this->nExpt(),
888  this->variableNames(), this->freeSpeciesNames(), this->fixdSpeciesNames(), this->twodimPDFs(),
889  this->splitSignal(), this->scfDPSmear());
891  splot.writeOutResults();
892  }
893 }
894 
895 void LauAbsFitModel::compareFitData(UInt_t toyMCScale, const TString& mcFileName, const TString& tableFileName, Bool_t poissonSmearing)
896 {
897  compareFitData_ = kTRUE;
898  fitToyMCScale_ = toyMCScale;
899  fitToyMCFileName_ = mcFileName;
900  fitToyMCTableName_ = tableFileName;
901  fitToyMCPoissonSmear_ = poissonSmearing;
902 }
903 
904 void LauAbsFitModel::createFitToyMC(const TString& mcFileName, const TString& tableFileName)
905 {
906  // Create a toy MC sample so that the user can compare the fitted
907  // result with the data.
908  // Generate more toy MC to reduce statistical fluctuations:
909  // - use the rescaling value fitToyMCScale_
910 
911  // Store the info on the number of experiments, first expt and current expt
912  UInt_t oldNExpt(this->nExpt());
913  UInt_t oldFirstExpt(this->firstExpt());
914  UInt_t oldIExpt(iExpt_);
915 
916  // Turn off Poisson smearing if required
917  Bool_t poissonSmearing(this->doPoissonSmearing());
919 
920  // Turn off embedding, since we need toy MC, not reco'ed events
921  Bool_t enableEmbeddingOrig(this->enableEmbedding());
922  this->enableEmbedding(kFALSE);
923 
924  // Need to make sure that the generation of the DP co-ordinates is
925  // switched on if any of our PDFs depend on it
926  Bool_t origUseDP = this->useDP();
927  if ( this->pdfsDependOnDP() && !origUseDP ) {
928  this->useDP( kTRUE );
929  this->initialiseDPModels();
930  }
931 
932  // Construct a unique filename for this experiment
933  TString exptString("_expt");
934  exptString += oldIExpt;
935  TString fileName( mcFileName );
936  fileName.Insert( fileName.Last('.'), exptString );
937 
938  // Generate the toy MC
939  std::cout << "INFO in LauAbsFitModel::createFitToyMC : Generating toy MC in " << fileName << " to compare fit with data..." << std::endl;
940  std::cout << " : Number of experiments to generate = " << fitToyMCScale_ << "." << std::endl;
941  std::cout << " : This is to allow the toy MC to be made with reduced statistical fluctuations." << std::endl;
942 
943  // Set the genValue of each parameter to its current (fitted) value
944  // but first store the original genValues for restoring later
945  std::vector<Double_t> origGenValues; origGenValues.reserve(nParams_);
946  for (LauParameterPList::iterator iter = fitVars_.begin(); iter != fitVars_.end(); ++iter) {
947  origGenValues.push_back((*iter)->genValue());
948  (*iter)->genValue((*iter)->value());
949  }
950 
951  // If we're asked to generate more than 100 experiments then split it
952  // up into multiple files since otherwise can run into memory issues
953  // when building the index
954 
955  UInt_t totalExpts = fitToyMCScale_;
956  if ( totalExpts > 100 ) {
957  UInt_t nFiles = totalExpts/100;
958  if ( totalExpts%100 ) {
959  nFiles += 1;
960  }
961  for ( UInt_t iFile(0); iFile < nFiles; ++iFile ) {
962 
963  UInt_t firstExp( iFile*100 );
964 
965  // Set number of experiments and first experiment to generate
966  UInt_t nExp = ((firstExp + 100)>totalExpts) ? totalExpts-firstExp : 100;
967  this->setNExpts(nExp, firstExp);
968 
969  // Create a unique filename and generate the events
970  TString extraname = "_file";
971  extraname += iFile;
972  fileName.Insert( fileName.Last('.'), extraname );
973  this->generate(fileName, "genResults", "dummy.root", tableFileName);
974  }
975  } else {
976  // Set number of experiments to new value
977  this->setNExpts(fitToyMCScale_, 0);
978  // Generate the toy
979  this->generate(fileName, "genResults", "dummy.root", tableFileName);
980  }
981 
982  // Reset number of experiments to original value
983  iExpt_ = oldIExpt;
984  this->setNExpts(oldNExpt, oldFirstExpt);
985 
986  // Restore the Poisson smearing to its former value
987  this->doPoissonSmearing(poissonSmearing);
988 
989  // Restore the embedding status
990  this->enableEmbedding(enableEmbeddingOrig);
991 
992  // Restore "useDP" to its former status
993  this->useDP( origUseDP );
994 
995  // Restore the original genValue to each parameter
996  for (UInt_t i(0); i<nParams_; ++i) {
997  fitVars_[i]->genValue(origGenValues[i]);
998  }
999 
1000  std::cout << "INFO in LauAbsFitModel::createFitToyMC : Finished in createFitToyMC." << std::endl;
1001 }
1002 
1004 {
1005  // Calculate the total negative log-likelihood over all events.
1006  // This function assumes that the fit parameters and data tree have
1007  // already been set-up correctly.
1008 
1009  // Loop over the data points to calculate the log likelihood
1010  Double_t logLike = this->getLogLikelihood( 0, this->eventsPerExpt() );
1011 
1012  // Include the Poisson term in the extended likelihood if required
1013  if (this->doEMLFit()) {
1014  logLike -= this->getEventSum();
1015  }
1016 
1017  // Calculate any penalty terms from Gaussian constrained variables
1018  if ( ! conVars_.empty() ){
1019  logLike -= this->getLogLikelihoodPenalty();
1020  }
1021 
1022  Double_t totNegLogLike = -logLike;
1023  return totNegLogLike;
1024 }
1025 
1027 {
1028  Double_t penalty(0.0);
1029 
1030  for ( LauAbsRValuePList::const_iterator iter = conVars_.begin(); iter != conVars_.end(); ++iter ) {
1031  Double_t val = (*iter)->value();
1032  Double_t mean = (*iter)->constraintMean();
1033  Double_t width = (*iter)->constraintWidth();
1034 
1035  Double_t term = ( val - mean )*( val - mean );
1036  penalty += term/( 2*width*width );
1037  }
1038 
1039  return penalty;
1040 }
1041 
1042 Double_t LauAbsFitModel::getLogLikelihood( UInt_t iStart, UInt_t iEnd )
1043 {
1044  // Calculate the total negative log-likelihood over all events.
1045  // This function assumes that the fit parameters and data tree have
1046  // already been set-up correctly.
1047 
1048  // Loop over the data points to calculate the log likelihood
1049  Double_t logLike(0.0);
1050 
1051  // Loop over the number of events in this experiment
1052  Bool_t ok(kTRUE);
1053  for (UInt_t iEvt = iStart; iEvt < iEnd; ++iEvt) {
1054 
1055  Double_t likelihood = this->getTotEvtLikelihood(iEvt);
1056 
1057  if (likelihood > DBL_MIN) { // Is the likelihood zero?
1058  Double_t evtLogLike = TMath::Log(likelihood);
1059  if ( doSFit_ ) {
1060  evtLogLike *= sWeights_[iEvt];
1061  }
1062  logLike += evtLogLike;
1063  } else {
1064  ok = kFALSE;
1065  std::cerr << "WARNING in LauAbsFitModel::getLogLikelihood : Strange likelihood value for event " << iEvt << ": " << likelihood << "\n";
1066  this->printEventInfo(iEvt);
1067  this->printVarsInfo(); //Write the values of the floated variables for which the likelihood is zero
1068  break;
1069  }
1070  }
1071 
1072  if (!ok) {
1073  std::cerr << " : Returning worst NLL found so far to force MINUIT out of this region." << std::endl;
1074  logLike = worstLogLike_;
1075  } else if (logLike < worstLogLike_) {
1076  worstLogLike_ = logLike;
1077  }
1078 
1079  return logLike;
1080 }
1081 
1082 void LauAbsFitModel::setParsFromMinuit(Double_t* par, Int_t npar)
1083 {
1084  // This function sets the internal parameters based on the values
1085  // that Minuit is using when trying to minimise the total likelihood function.
1086 
1087  // MINOS reports different numbers of free parameters depending on the
1088  // situation, so disable this check
1089  if ( ! withinAsymErrorCalc_ ) {
1090  if (static_cast<UInt_t>(npar) != nFreeParams_) {
1091  std::cerr << "ERROR in LauAbsFitModel::setParsFromMinuit : Unexpected number of free parameters: " << npar << ".\n";
1092  std::cerr << " Expected: " << nFreeParams_ << ".\n" << std::endl;
1093  gSystem->Exit(EXIT_FAILURE);
1094  }
1095  }
1096 
1097  // Despite npar being the number of free parameters
1098  // the par array actually contains all the parameters,
1099  // free and floating...
1100  // Update all the floating ones with their new values.
1101  for (UInt_t i(0); i<nParams_; ++i) {
1102  if (!fitVars_[i]->fixed()) {
1103  fitVars_[i]->value(par[i]);
1104  }
1105  }
1106 
1107  this->propagateParUpdates();
1108 }
1109 
1111 {
1112  UInt_t nParsAdded(0);
1113  for (LauPdfList::iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) {
1114  LauAbsPdf* pdf = (*pdf_iter);
1115  if ( pdf->isDPDependent() ) {
1116  this->pdfsDependOnDP( kTRUE );
1117  }
1118  LauAbsRValuePList& pars = pdf->getParameters();
1119  for (LauAbsRValuePList::iterator pars_iter = pars.begin(); pars_iter != pars.end(); ++pars_iter) {
1120  LauParameterPList params = (*pars_iter)->getPars();
1121  for (LauParameterPList::iterator params_iter = params.begin(); params_iter != params.end(); ++params_iter) {
1122  if ( !(*params_iter)->clone() && ( !(*params_iter)->fixed() ||
1123  (this->twoStageFit() && ( (*params_iter)->secondStage() || (*params_iter)->firstStage())) ) ) {
1124  fitVars_.push_back(*params_iter);
1125  ++nParsAdded;
1126  }
1127  }
1128  }
1129  }
1130  return nParsAdded;
1131 }
1132 
1133 void LauAbsFitModel::addConstraint(const TString& formula, const std::vector<TString>& pars, const Double_t mean, const Double_t width)
1134 {
1135  StoreConstraints newCon;
1136  newCon.formula_ = formula;
1137  newCon.conPars_ = pars;
1138  newCon.mean_ = mean;
1139  newCon.width_ = width;
1140  storeCon_.push_back(newCon);
1141 }
1142 
1144 {
1145  for ( LauParameterPList::const_iterator iter = fitVars_.begin(); iter != fitVars_.end(); ++iter ) {
1146  if ( (*iter)->gaussConstraint() ) {
1147  conVars_.push_back( *iter );
1148  std::cout << "INFO in LauAbsFitModel::addConParameters : Added Gaussian constraint to parameter "<< (*iter)->name() << std::endl;
1149  }
1150  }
1151 
1152  // Add penalties from the constraints to fit parameters
1153  for ( std::vector<StoreConstraints>::iterator iter = storeCon_.begin(); iter != storeCon_.end(); ++iter ) {
1154  std::vector<TString> names = (*iter).conPars_;
1155  std::vector<LauParameter*> params;
1156  for ( std::vector<TString>::iterator iternames = names.begin(); iternames != names.end(); ++iternames ) {
1157  for ( LauParameterPList::const_iterator iterfit = fitVars_.begin(); iterfit != fitVars_.end(); ++iterfit ) {
1158  if ( (*iternames) == (*iterfit)->name() ){
1159  params.push_back(*iterfit);
1160  }
1161  }
1162  }
1163 
1164  // If the parameters are not found, skip it
1165  if ( params.size() != (*iter).conPars_.size() ) {
1166  std::cerr << "WARNING in LauAbsFitModel::addConParameters: Could not find parameters to constrain in the formula... skipping" << std::endl;
1167  continue;
1168  }
1169 
1170  LauFormulaPar* formPar = new LauFormulaPar( (*iter).formula_, (*iter).formula_, params );
1171  formPar->addGaussianConstraint( (*iter).mean_, (*iter).width_ );
1172  conVars_.push_back(formPar);
1173 
1174  std::cout << "INFO in LauAbsFitModel::addConParameters : Added Gaussian constraint to formula\n";
1175  std::cout << " : Formula: " << (*iter).formula_ << std::endl;
1176  for ( std::vector<LauParameter*>::iterator iterparam = params.begin(); iterparam != params.end(); ++iterparam ) {
1177  std::cout << " : Parameter: " << (*iterparam)->name() << std::endl;
1178  }
1179  }
1180 
1181 }
1182 
1184 {
1185  for (LauPdfList::iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) {
1186  (*pdf_iter)->updatePulls();
1187  }
1188 }
1189 
1190 void LauAbsFitModel::printFitParameters(const LauPdfList& pdfList, std::ostream& fout) const
1191 {
1192  LauPrint print;
1193  for (LauPdfList::const_iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) {
1194  const LauAbsRValuePList& pars = (*pdf_iter)->getParameters();
1195  for (LauAbsRValuePList::const_iterator pars_iter = pars.begin(); pars_iter != pars.end(); ++pars_iter) {
1196  LauParameterPList params = (*pars_iter)->getPars();
1197  for (LauParameterPList::iterator params_iter = params.begin(); params_iter != params.end(); ++params_iter) {
1198  if (!(*params_iter)->clone()) {
1199  fout << (*params_iter)->name() << " & $";
1200  print.printFormat(fout, (*params_iter)->value());
1201  if ((*params_iter)->fixed() == kTRUE) {
1202  fout << "$ (fixed) \\\\";
1203  } else {
1204  fout << " \\pm ";
1205  print.printFormat(fout, (*params_iter)->error());
1206  fout << "$ \\\\" << std::endl;
1207  }
1208  }
1209  }
1210  }
1211  }
1212 }
1213 
1215 {
1216  for (LauPdfList::iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) {
1217  (*pdf_iter)->cacheInfo(theData);
1218  }
1219 }
1220 
1221 Double_t LauAbsFitModel::prodPdfValue(LauPdfList& pdfList, UInt_t iEvt)
1222 {
1223  Double_t pdfVal = 1.0;
1224  for (LauPdfList::iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) {
1225  (*pdf_iter)->calcLikelihoodInfo(iEvt);
1226  pdfVal *= (*pdf_iter)->getLikelihood();
1227  }
1228  return pdfVal;
1229 }
1230 
1231 void LauAbsFitModel::printEventInfo(UInt_t iEvt) const
1232 {
1233  const LauFitData& data = inputFitData_->getData(iEvt);
1234  std::cerr << " : Input data values for this event:" << std::endl;
1235  for (LauFitData::const_iterator iter = data.begin(); iter != data.end(); ++iter) {
1236  std::cerr << " " << iter->first << " = " << iter->second << std::endl;
1237  }
1238 }
1239 
1241 {
1242  std::cerr << " : Current values of fit parameters:" << std::endl;
1243  for (UInt_t i(0); i<nParams_; ++i) {
1244  std::cerr << " " << (fitVars_[i]->name()).Data() << " = " << fitVars_[i]->value() << std::endl;
1245  }
1246 }
1247 
Class for defining the SPlot technique.
Definition: LauSPlot.hh:55
TString fitToyMCTableName_
The output table name for Toy MC.
Double_t mean_
The mean value of the Gaussian constraint to be applied.
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 twoStageFit() const
Determine whether the two-stage fit is enabled.
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.
virtual void updateParameters()=0
Update the values and errors of the parameters based on the fit minimum.
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.
Double_t * parValues_
Parameter values array (for reading from the master)
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()
Calculates the total negative log-likelihood.
LauAbsRValuePList conVars_
Internal vectors of Gaussian parameters.
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.
void runCalculations(const TString &option="q")
Method to calculate the sWeights and cN coeffs.
Definition: LauSPlot.cc:626
virtual void setGenNtupleIntegerBranchValue(const TString &name, Int_t value)
Set the value of an integer branch in the gen tree.
ClassImp(LauAbsCoeffSet)
virtual void setupBkgndVectors()=0
Method to set up the storage for background-related quantities called by setBkgndClassNames.
LauGenNtuple * genNtuple_
The generated ntuple.
File containing declaration of LauAbsFitModel class.
Bool_t withinAsymErrorCalc_
Flag to indicate if the asymmetric error calculation (e.g. MINOS) is currently running.
TString formula_
The formula to be used in the LauFormulaPar.
Bool_t useAsymmFitErrors() const
Determine whether calculation of asymmetric errors is enabled.
const TString & name() const
The parameter name.
Class for defining combinations of fit parameter objects.
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.
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
std::vector< TString > conPars_
The list of LauParameter names to be used in the LauFormulaPar.
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 Bool_t twoStageFit() const =0
Determine whether the two-stage fit is enabled.
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.
File containing declaration of LauAbsFitter class.
TSocket * sMaster_
A socket to enable parallel setup.
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.
virtual void addConstraint(const TString &formula, const std::vector< TString > &pars, const Double_t mean, const Double_t width)
Store constraint information for fit parameters.
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 total 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.
virtual void fillGenNtupleBranches()
Fill the gen tuple branches.
void fitSlave(const TString &dataFileName, const TString &dataTreeName, const TString &histFileName, const TString &tableFileNameBase)
Slaves required when performing a simultaneous fit.
virtual void setSPlotNtupleDoubleBranchValue(const TString &name, Double_t value)
Set the value of a double branch in the sPlot tree.
void printFitParameters(const LauPdfList &pdfList, std::ostream &fout) const
Print the fit parameters for all PDFs in the list.
void addGaussianConstraint(Double_t newGaussMean, Double_t newGaussWidth)
Add a Gaussian constraint (or modify an existing one)
virtual void initialise(LauFitObject *fitObj, const std::vector< LauParameter * > &parameters)=0
Initialise the fitter, setting the information on the parameters.
void readExperimentData(UInt_t iExpt)
Read events only for the given experiment.
virtual std::pair< Int_t, Double_t > minimise()=0
Perform the minimisation of the fit function.
std::vector< LauAbsRValue * > LauAbsRValuePList
List of parameter pointers.
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 class.
virtual void propagateParUpdates()=0
This function (specific to each model) calculates anything that depends on the fit parameter values...
void addConParameters()
Add parameters to the list of Gaussian constrained parameters.
LauFitDataTree * inputFitData_
The input data.
Bool_t doSFit() const
Return the flag to store the status of using an sFit or not.
void run(const TString &applicationCode, const TString &dataFileName, const TString &dataTreeName, const TString &histFileName, const TString &tableFileName="")
Start the toy generation / fitting.
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.
std::vector< StoreConstraints > storeCon_
Store the constraints for fit parameters until initialisation is complete.
TMatrixD covMatrix_
The fit covariance matrix.
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.
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.
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.
virtual void fixFirstStageParameters()=0
Fix parameters marked as &quot;first stage&quot;.
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:111
virtual UInt_t nParameters() const =0
Get the total number of fit parameters.
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 worstLogLike_
The worst LL value found so far.
virtual void releaseSecondStageParameters()=0
Release parameters marked as &quot;second stage&quot;.
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.
Class for defining the fit parameter objects.
Definition: LauParameter.hh:33
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.
Struct to store constraint information until the fit is run.
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.
virtual const std::vector< LauAbsRValue * > & getParameters() const
Retrieve the parameters of the PDF, e.g. so that they can be loaded into a fit.
Definition: LauAbsPdf.hh:239
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.
UInt_t firstExpt_
The number of the first experiment to consider.
static LauAbsFitter * fitter()
Method that provides access to the singleton fitter.
Definition: LauFitter.cc:34
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.
virtual UInt_t nFreeParameters() const =0
Get the number of floating fit parameters.
UInt_t iExpt_
The number of the current experiment.
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 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.
UInt_t fitToyMCScale_
The scaling factor (toy vs data statistics)
virtual void fillSPlotNtupleBranches()
Fill the sPlot tuple.
TString fitToyMCFileName_
The output file name for Toy MC.
Double_t width_
The width of the Gaussian constraint to be applied.
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
void runSlave(const TString &dataFileName, const TString &dataTreeName, const TString &histFileName, const TString &tableFileName="", const TString &addressMaster="localhost", const UInt_t portMaster=9090)
Start the slave process for simultaneous fitting.
virtual LauSPlot::TwoDMap twodimPDFs() const =0
Returns the species and variables for all 2D PDFs in the fit.
Double_t getLogLikelihoodPenalty()
Calculate the penalty terms to the log likelihood from Gaussian constraints.
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.
virtual Bool_t useAsymmFitErrors() const =0
Determine whether calculation of asymmetric errors is enabled.
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:41
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.
virtual void releaseFirstStageParameters()=0
Release parameters marked as &quot;first stage&quot;.
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 const TMatrixD & covarianceMatrix() const =0
Retrieve the fit covariance matrix.
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.