laura is hosted by Hepforge, IPPP Durham
Laura++  v3r0p1
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauSimpleFitModel.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 <iomanip>
17 #include <fstream>
18 
19 #include "TFile.h"
20 #include "TH2.h"
21 #include "TMinuit.h"
22 #include "TRandom.h"
23 #include "TSystem.h"
24 #include "TVirtualFitter.h"
25 
26 #include "LauAbsBkgndDPModel.hh"
27 #include "LauAbsCoeffSet.hh"
28 #include "LauIsobarDynamics.hh"
29 #include "LauAbsPdf.hh"
30 #include "LauComplex.hh"
31 #include "LauConstants.hh"
32 #include "LauDaughters.hh"
33 #include "LauEmbeddedData.hh"
34 #include "LauEffModel.hh"
35 #include "LauFitNtuple.hh"
36 #include "LauGenNtuple.hh"
37 #include "LauKinematics.hh"
38 #include "LauPrint.hh"
39 #include "LauRandom.hh"
40 #include "LauScfMap.hh"
41 #include "LauSimpleFitModel.hh"
42 
44 
45 
47  sigDPModel_(sigModel),
48  kinematics_(sigModel ? sigModel->getKinematics() : 0),
49  usingBkgnd_(kFALSE),
50  nSigComp_(0),
51  nSigDPPar_(0),
52  nExtraPdfPar_(0),
53  nNormPar_(0),
54  meanEff_("meanEff",0.0,0.0,1.0),
55  dpRate_("dpRate",0.0,0.0,100.0),
56  //signalEvents_("signalEvents",1.0,0.0,1.0),
57  signalEvents_(0),
58  useSCF_(kFALSE),
59  useSCFHist_(kFALSE),
60  scfFrac_("scfFrac",0.0,0.0,1.0),
61  scfFracHist_(0),
62  scfMap_(0),
63  compareFitData_(kFALSE),
64  signalTree_(0),
65  reuseSignal_(kFALSE),
66  useReweighting_(kFALSE),
67  sigDPLike_(0.0),
68  scfDPLike_(0.0),
69  sigExtraLike_(0.0),
70  scfExtraLike_(0.0),
71  sigTotalLike_(0.0),
72  scfTotalLike_(0.0)
73 {
74 }
75 
77 {
78  delete signalTree_;
79  for (LauBkgndEmbDataList::iterator iter = bkgndTree_.begin(); iter != bkgndTree_.end(); ++iter) {
80  delete (*iter);
81  }
82  delete scfFracHist_;
83 }
84 
86 {
87  UInt_t nBkgnds = this->nBkgndClasses();
88  bkgndDPModels_.resize( nBkgnds );
89  bkgndPdfs_.resize( nBkgnds );
90  bkgndEvents_.resize( nBkgnds );
91  bkgndTree_.resize( nBkgnds );
92  reuseBkgnd_.resize( nBkgnds );
93  bkgndDPLike_.resize( nBkgnds );
94  bkgndExtraLike_.resize( nBkgnds );
95  bkgndTotalLike_.resize( nBkgnds );
96 }
97 
99 {
100  if ( nSigEvents == 0 ) {
101  std::cerr << "ERROR in LauSimpleFitModel::setNSigEvents : The signal yield LauParameter pointer is null." << std::endl;
102  gSystem->Exit(EXIT_FAILURE);
103  }
104  if ( signalEvents_ != 0 ) {
105  std::cerr << "ERROR in LauSimpleFitModel::setNSigEvents : You are trying to overwrite the signal yield." << std::endl;
106  return;
107  }
108 
109  signalEvents_ = nSigEvents;
110  TString name = signalEvents_->name();
111  if ( ! name.Contains("signalEvents") && !( name.BeginsWith("signal") && name.EndsWith("Events") ) ) {
112  signalEvents_->name("signalEvents");
113  }
114  Double_t value = nSigEvents->value();
115  signalEvents_->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0));
116 }
117 
119 {
120  if ( nBkgndEvents == 0 ) {
121  std::cerr << "ERROR in LauSimpleFitModel::setNBkgndEvents : The background yield LauParameter pointer is null." << std::endl;
122  gSystem->Exit(EXIT_FAILURE);
123  }
124 
125  if ( ! this->validBkgndClass( nBkgndEvents->name() ) ) {
126  std::cerr << "ERROR in LauSimpleFitModel::setNBkgndEvents : Invalid background class \"" << nBkgndEvents->name() << "\"." << std::endl;
127  std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
128  gSystem->Exit(EXIT_FAILURE);
129  }
130 
131  UInt_t bkgndID = this->bkgndClassID( nBkgndEvents->name() );
132 
133  if ( bkgndEvents_[bkgndID] != 0 ) {
134  std::cerr << "ERROR in LauSimpleFitModel::setNBkgndEvents : You are trying to overwrite the background yield." << std::endl;
135  return;
136  }
137 
138  bkgndEvents_[bkgndID] = nBkgndEvents;
139  bkgndEvents_[bkgndID]->name( nBkgndEvents->name()+"Events" );
140  Double_t value = nBkgndEvents->value();
141  bkgndEvents_[bkgndID]->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0));
142 }
143 
144 void LauSimpleFitModel::splitSignalComponent( const TH2* dpHisto, const Bool_t upperHalf, const Bool_t fluctuateBins, LauScfMap* scfMap )
145 {
146  if ( useSCF_ == kTRUE ) {
147  std::cerr << "ERROR in LauSimpleFitModel::splitSignalComponent : Have already setup SCF." << std::endl;
148  return;
149  }
150 
151  if ( dpHisto == 0 ) {
152  std::cerr << "ERROR in LauSimpleFitModel::splitSignalComponent : The histogram pointer is null." << std::endl;
153  return;
154  }
155 
156  LauDaughters* daughters = sigDPModel_->getDaughters();
157  scfFracHist_ = new LauEffModel( daughters, 0 );
158  scfFracHist_->setEffHisto( dpHisto, kTRUE, fluctuateBins, 0.0, 0.0, upperHalf, daughters->squareDP() );
159 
160  scfMap_ = scfMap;
161 
162  useSCF_ = kTRUE;
163  useSCFHist_ = kTRUE;
164 }
165 
166 void LauSimpleFitModel::splitSignalComponent( const Double_t scfFrac, const Bool_t fixed )
167 {
168  if ( useSCF_ == kTRUE ) {
169  std::cerr << "ERROR in LauSimpleFitModel::splitSignalComponent : Have already setup SCF." << std::endl;
170  return;
171  }
172 
173  scfFrac_.range( 0.0, 1.0 );
174  scfFrac_.value( scfFrac ); scfFrac_.initValue( scfFrac ); scfFrac_.genValue( scfFrac );
175  scfFrac_.fixed( fixed );
176 
177  useSCF_ = kTRUE;
178  useSCFHist_ = kFALSE;
179 }
180 
181 void LauSimpleFitModel::setBkgndDPModel(const TString& bkgndClass, LauAbsBkgndDPModel* bkgndDPModel)
182 {
183  if (bkgndDPModel == 0) {
184  std::cerr << "ERROR in LauSimpleFitModel::setBkgndDPModel : The model pointer is null." << std::endl;
185  return;
186  }
187 
188  // check that this background name is valid
189  if ( ! this->validBkgndClass( bkgndClass ) ) {
190  std::cerr << "ERROR in LauSimpleFitModel::setBkgndDPModel : Invalid background class \"" << bkgndClass << "\"." << std::endl;
191  std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
192  return;
193  }
194 
195  UInt_t bkgndID = this->bkgndClassID( bkgndClass );
196  bkgndDPModels_[bkgndID] = bkgndDPModel;
197 
198  usingBkgnd_ = kTRUE;
199 }
200 
202 {
203  if (pdf==0) {
204  std::cerr << "ERROR in LauSimpleFitModel::setSignalPdf : The PDF pointer is null." << std::endl;
205  return;
206  }
207  signalPdfs_.push_back(pdf);
208 }
209 
211 {
212  if (pdf==0) {
213  std::cerr << "ERROR in LauSimpleFitModel::setSCFPdf : The PDF pointer is null." << std::endl;
214  return;
215  }
216  scfPdfs_.push_back(pdf);
217 }
218 
219 void LauSimpleFitModel::setBkgndPdf(const TString& bkgndClass, LauAbsPdf* pdf)
220 {
221  if (pdf == 0) {
222  std::cerr << "ERROR in LauSimpleFitModel::setBkgndPdf : The PDF pointer is null." << std::endl;
223  return;
224  }
225 
226  // check that this background name is valid
227  if ( ! this->validBkgndClass( bkgndClass ) ) {
228  std::cerr << "ERROR in LauSimpleFitModel::setBkgndPdf : Invalid background class \"" << bkgndClass << "\"." << std::endl;
229  std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
230  return;
231  }
232 
233  UInt_t bkgndID = this->bkgndClassID( bkgndClass );
234  bkgndPdfs_[bkgndID].push_back(pdf);
235 
236  usingBkgnd_ = kTRUE;
237 }
238 
240 {
241  // Is there a component called compName in the signal model?
242  TString compName(coeffSet->name());
243  Bool_t ok = sigDPModel_->hasResonance(compName);
244  if (!ok) {
245  std::cerr << "ERROR in LauSimpleFitModel::setAmpCoeffSet : Signal DP model doesn't contain component \"" << compName << "\"." << std::endl;
246  return;
247  }
248 
249  // Do we already have it in our list of names?
250  for (std::vector<LauAbsCoeffSet*>::const_iterator iter = coeffPars_.begin(); iter != coeffPars_.end(); ++iter) {
251  if ((*iter)->name() == compName) {
252  std::cerr << "ERROR in LauSimpleFitModel::setAmpCoeffSet : Have already set coefficients for \"" << compName << "\"." << std::endl;
253  return;
254  }
255  }
256 
257  coeffSet->index(nSigComp_);
258  coeffPars_.push_back(coeffSet);
259 
260  ++nSigComp_;
261 
262  std::cout << "INFO in LauSimpleFitModel::setAmpCoeffSet : Added coefficients for component \"" << compName << "\" to the fit model." << std::endl;
263  coeffSet->printParValues();
264 }
265 
266 
268 {
269  // From the initial parameter values calculate the coefficients
270  // so they can be passed to the signal model
271  this->updateCoeffs();
272 
273  // Initialisation
274  if (this->useDP() == kTRUE) {
275  this->initialiseDPModels();
276  }
277 
278  if (!this->useDP() && signalPdfs_.empty()) {
279  std::cerr << "ERROR in LauSimpleFitModel::initialise : Signal model doesn't exist for any variable." << std::endl;
280  gSystem->Exit(EXIT_FAILURE);
281  }
282 
283  if (this->useDP()) {
284  // Check that we have all the Dalitz-plot models
285  if ( sigDPModel_ == 0 ) {
286  std::cerr << "ERROR in LauSimpleFitModel::initialise : The pointer to the signal DP model is null.\n";
287  std::cerr << " : Removing the Dalitz Plot from the model." << std::endl;
288  this->useDP(kFALSE);
289  }
290  if ( usingBkgnd_ ) {
291  for (LauBkgndDPModelList::const_iterator dpmodel_iter = bkgndDPModels_.begin(); dpmodel_iter != bkgndDPModels_.end(); ++dpmodel_iter ) {
292  if ( (*dpmodel_iter) == 0 ) {
293  std::cerr << "ERROR in LauSimpleFitModel::initialise : The pointer to one of the background DP models is null.\n";
294  std::cerr << " : Removing the Dalitz Plot from the model." << std::endl;
295  this->useDP(kFALSE);
296  break;
297  }
298  }
299  }
300  }
301 
302  // Next check that, if a given component is being used we've got the
303  // right number of PDFs for all the variables involved
304  // TODO - should probably check variable names and so on as well
305 
306  UInt_t nsigpdfvars(0);
307  for ( LauPdfList::const_iterator pdf_iter = signalPdfs_.begin(); pdf_iter != signalPdfs_.end(); ++pdf_iter ) {
308  std::vector<TString> varNames = (*pdf_iter)->varNames();
309  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
310  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
311  ++nsigpdfvars;
312  }
313  }
314  }
315  if (useSCF_) {
316  UInt_t nscfpdfvars(0);
317  for ( LauPdfList::const_iterator pdf_iter = scfPdfs_.begin(); pdf_iter != scfPdfs_.end(); ++pdf_iter ) {
318  std::vector<TString> varNames = (*pdf_iter)->varNames();
319  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
320  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
321  ++nscfpdfvars;
322  }
323  }
324  }
325  if (nscfpdfvars != nsigpdfvars) {
326  std::cerr << "ERROR in LauSimpleFitModel::initialise : There are " << nsigpdfvars << " TM signal PDF variables but " << nscfpdfvars << " SCF signal PDF variables." << std::endl;
327  gSystem->Exit(EXIT_FAILURE);
328  }
329  }
330  if (usingBkgnd_) {
331  for (LauBkgndPdfsList::const_iterator bgclass_iter = bkgndPdfs_.begin(); bgclass_iter != bkgndPdfs_.end(); ++bgclass_iter) {
332  UInt_t nbkgndpdfvars(0);
333  const LauPdfList& pdfList = (*bgclass_iter);
334  for ( LauPdfList::const_iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter ) {
335  std::vector<TString> varNames = (*pdf_iter)->varNames();
336  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
337  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
338  ++nbkgndpdfvars;
339  }
340  }
341  }
342  if (nbkgndpdfvars != nsigpdfvars) {
343  std::cerr << "ERROR in LauSimpleFitModel::initialise : There are " << nsigpdfvars << " signal PDF variables but " << nbkgndpdfvars << " bkgnd PDF variables." << std::endl;
344  gSystem->Exit(EXIT_FAILURE);
345  }
346  }
347  }
348 
349  // Clear the vectors of parameter information so we can start from scratch
350  this->clearFitParVectors();
351 
352  // Set the fit parameters for signal and background models
353  this->setSignalDPParameters();
354 
355  // Set the fit parameters for the various extra PDFs
356  this->setExtraPdfParameters();
357 
358  // Set the initial bg and signal events
359  this->setFitNEvents();
360 
361  // Check that we have the expected number of fit variables
362  const LauParameterPList& fitVars = this->fitPars();
363  if (fitVars.size() != (nSigDPPar_ + nExtraPdfPar_ + nNormPar_)) {
364  std::cerr << "ERROR in LauSimpleFitModel::initialise : Number of fit parameters not of expected size." << std::endl;
365  gSystem->Exit(EXIT_FAILURE);
366  }
367 
368  this->setExtraNtupleVars();
369 }
370 
372 {
373  // Need to check that the number of components we have and that the dynamics has matches up
374  UInt_t nAmp = sigDPModel_->getnTotAmp();
375  if (nAmp != nSigComp_) {
376  std::cerr << "ERROR in LauSimpleFitModel::initialiseDPModels : Number of signal DP components with magnitude and phase set not right." << std::endl;
377  gSystem->Exit(EXIT_FAILURE);
378  }
379 
380  std::cout << "INFO in LauSimpleFitModel::initialiseDPModels : Initialising DP models" << std::endl;
381 
383 
384  if (usingBkgnd_ == kTRUE) {
385  for (LauBkgndDPModelList::iterator iter = bkgndDPModels_.begin(); iter != bkgndDPModels_.end(); ++iter) {
386  (*iter)->initialise();
387  }
388  }
389 }
390 
392 {
393  //std::cout << "INFO in LauSimpleFitModel::recalculateNormalizationInDPModels : Recalc Norm in DP model" << std::endl;
396 }
397 
399 {
400  // Set the fit parameters for the signal model.
401 
402  nSigDPPar_ = 0;
403  if ( ! this->useDP() ) {
404  return;
405  }
406 
407  std::cout << "INFO in LauSimpleFitModel::setSignalDPParameters : Setting the initial fit parameters for the signal DP model." << std::endl;
408 
409  // Place isobar coefficient parameters in vector of fit variables
410  LauParameterPList& fitVars = this->fitPars();
411  for (UInt_t i = 0; i < nSigComp_; i++) {
412  LauParameterPList pars = coeffPars_[i]->getParameters();
413  for (LauParameterPList::iterator iter = pars.begin(); iter != pars.end(); ++iter) {
414  if ( !(*iter)->clone() ) {
415  fitVars.push_back(*iter);
416  ++nSigDPPar_;
417  }
418  }
419  }
420 
421  // Obtain the resonance parameters and place them in the vector of fit variables and in a separate vector
422  LauParameterPSet& resVars = this->resPars();
423  resVars.clear();
424 
426 
427  for ( LauParameterPList::iterator iter = sigDPPars.begin(); iter != sigDPPars.end(); ++iter ) {
428  if ( resVars.insert(*iter).second ) {
429  fitVars.push_back(*iter);
430  ++nSigDPPar_;
431  }
432  }
433 }
434 
436 {
437  // Include all the parameters of the various PDFs in the fit
438  // NB all of them are passed to the fit, even though some have been fixed through parameter.fixed(kTRUE)
439  // With the new "cloned parameter" scheme only "original" parameters are passed to the fit.
440  // Their clones are updated automatically when the originals are updated.
441 
442  nExtraPdfPar_ = 0;
443 
445 
446  if (useSCF_ == kTRUE) {
448  }
449 
450  if (usingBkgnd_ == kTRUE) {
451  for (LauBkgndPdfsList::iterator iter = bkgndPdfs_.begin(); iter != bkgndPdfs_.end(); ++iter) {
452  nExtraPdfPar_ += this->addFitParameters(*iter);
453  }
454  }
455 }
456 
458 {
459  if ( signalEvents_ == 0 ) {
460  std::cerr << "ERROR in LauSimpleFitModel::setFitNEvents : Signal yield not defined." << std::endl;
461  return;
462  }
463  nNormPar_ = 0;
464 
465  // initialise the total number of events to be the sum of all the hypotheses
466  Double_t nTotEvts = signalEvents_->value();
467  for (LauBkgndYieldList::const_iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
468  nTotEvts += (*iter)->value();
469  if ( (*iter) == 0 ) {
470  std::cerr << "ERROR in LauSimpleFitModel::setFitNEvents : Background yield not defined." << std::endl;
471  return;
472  }
473  }
474  this->eventsPerExpt(TMath::FloorNint(nTotEvts));
475 
476  LauParameterPList& fitVars = this->fitPars();
477 
478  // if doing an extended ML fit add the number of signal events into the fit parameters
479  if (this->doEMLFit()) {
480  std::cout << "INFO in LauSimpleFitModel::setFitNEvents : Initialising number of events for signal and background components..." << std::endl;
481  // add the signal events to the list of fit parameters
482  fitVars.push_back(signalEvents_);
483  ++nNormPar_;
484  } else {
485  std::cout << "INFO in LauSimpleFitModel::setFitNEvents : Initialising number of events for background components (and hence signal)..." << std::endl;
486  }
487 
488  if (useSCF_ && !useSCFHist_) {
489  fitVars.push_back(&scfFrac_);
490  ++nNormPar_;
491  }
492 
493  if (usingBkgnd_ == kTRUE) {
494  for (LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
495  LauParameter* parameter = (*iter);
496  fitVars.push_back(parameter);
497  ++nNormPar_;
498  }
499  }
500 }
501 
503 {
504  // Set-up other parameters derived from the fit results, e.g. fit fractions.
505  if (this->useDP() != kTRUE) {
506  return;
507  }
508 
509  // First clear the vectors so we start from scratch
510  this->clearExtraVarVectors();
511 
512  LauParameterList& extraVars = this->extraPars();
513 
514  // Add a fit fraction for each signal component
516  if (fitFrac_.size() != nSigComp_) {
517  std::cerr << "ERROR in LauSimpleFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << fitFrac_.size() << std::endl;
518  gSystem->Exit(EXIT_FAILURE);
519  }
520  for (UInt_t i(0); i<nSigComp_; ++i) {
521  if (fitFrac_[i].size() != nSigComp_) {
522  std::cerr << "ERROR in LauSimpleFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << fitFrac_[i].size() << std::endl;
523  gSystem->Exit(EXIT_FAILURE);
524  }
525  }
526 
527  // Add the fit fraction that has not been corrected for the efficiency for each signal component
529  if (fitFracEffUnCorr_.size() != nSigComp_) {
530  std::cerr << "ERROR in LauSimpleFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << fitFracEffUnCorr_.size() << std::endl;
531  gSystem->Exit(EXIT_FAILURE);
532  }
533  for (UInt_t i(0); i<nSigComp_; ++i) {
534  if (fitFracEffUnCorr_[i].size() != nSigComp_) {
535  std::cerr << "ERROR in LauSimpleFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << fitFracEffUnCorr_[i].size() << std::endl;
536  gSystem->Exit(EXIT_FAILURE);
537  }
538  }
539 
540  for (UInt_t i = 0; i < nSigComp_; i++) {
541  for (UInt_t j = i; j < nSigComp_; j++) {
542  extraVars.push_back(fitFrac_[i][j]);
543  extraVars.push_back(fitFracEffUnCorr_[i][j]);
544  }
545  }
546 
547  // Store any extra parameters/quantities from the DP model (e.g. K-matrix total fit fractions)
548  std::vector<LauParameter> extraParams = sigDPModel_->getExtraParameters();
549  std::vector<LauParameter>::iterator extraIter;
550  for (extraIter = extraParams.begin(); extraIter != extraParams.end(); ++extraIter) {
551  LauParameter extraParameter = (*extraIter);
552  extraVars.push_back(extraParameter);
553  }
554 
555  // Now add in the DP efficiency value
556  Double_t initMeanEff = sigDPModel_->getMeanEff().initValue();
557  meanEff_.value(initMeanEff); meanEff_.initValue(initMeanEff); meanEff_.genValue(initMeanEff);
558  extraVars.push_back(meanEff_);
559 
560  // Also add in the DP rate
561  Double_t initDPRate = sigDPModel_->getDPRate().initValue();
562  dpRate_.value(initDPRate); dpRate_.initValue(initDPRate); dpRate_.genValue(initDPRate);
563  extraVars.push_back(dpRate_);
564 }
565 
566 void LauSimpleFitModel::finaliseFitResults(const TString& tablePrefixName)
567 {
568  // Retrieve parameters from the fit results for calculations and toy generation
569  // and eventually store these in output root ntuples/text files
570 
571  // Now take the fit parameters and update them as necessary
572  // e.g. to make mag > 0.0 and phase in the right range.
573  // This function will also calculate any other values, such as the
574  // fit fractions, using any errors provided by fitParErrors as appropriate.
575  // Also obtain the pull values: (measured - generated)/(average error)
576 
577  if (this->useDP() == kTRUE) {
578  for (UInt_t i = 0; i < nSigComp_; i++) {
579  // Check whether we have mag > 0.0, and phase in the right range
580  coeffPars_[i]->finaliseValues();
581  }
582  }
583 
584  // update the pulls on the events
585  if (this->doEMLFit()) {
587  }
588  if (useSCF_ && !useSCFHist_) {
590  }
591  if (usingBkgnd_ == kTRUE) {
592  for (LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
593  (*iter)->updatePull();
594  }
595  }
596 
597  // Update the pulls on all the extra PDFs' parameters
599  if (useSCF_ == kTRUE) {
601  }
602  if (usingBkgnd_ == kTRUE) {
603  for (LauBkgndPdfsList::iterator iter = bkgndPdfs_.begin(); iter != bkgndPdfs_.end(); ++iter) {
604  this->updateFitParameters(*iter);
605  }
606  }
607 
608  // Fill the fit results to the ntuple for current experiment
609 
610  // update the coefficients and then calculate the fit fractions
611  if (this->useDP() == kTRUE) {
612  this->updateCoeffs();
615 
617  if (fitFrac.size() != nSigComp_) {
618  std::cerr << "ERROR in LauSimpleFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << fitFrac.size() << std::endl;
619  gSystem->Exit(EXIT_FAILURE);
620  }
621  for (UInt_t i(0); i<nSigComp_; ++i) {
622  if (fitFrac[i].size() != nSigComp_) {
623  std::cerr << "ERROR in LauSimpleFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << fitFrac[i].size() << std::endl;
624  gSystem->Exit(EXIT_FAILURE);
625  }
626  }
627 
629  if (fitFracEffUnCorr.size() != nSigComp_) {
630  std::cerr << "ERROR in LauSimpleFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << fitFracEffUnCorr.size() << std::endl;
631  gSystem->Exit(EXIT_FAILURE);
632  }
633  for (UInt_t i(0); i<nSigComp_; ++i) {
634  if (fitFracEffUnCorr[i].size() != nSigComp_) {
635  std::cerr << "ERROR in LauSimpleFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << fitFracEffUnCorr[i].size() << std::endl;
636  gSystem->Exit(EXIT_FAILURE);
637  }
638  }
639 
640  for (UInt_t i(0); i<nSigComp_; ++i) {
641  for (UInt_t j(i); j<nSigComp_; ++j) {
642  fitFrac_[i][j].value(fitFrac[i][j].value());
643  fitFracEffUnCorr_[i][j].value(fitFracEffUnCorr[i][j].value());
644  }
645  }
646 
649 
650  this->clearExtraVarVectors();
651  LauParameterList& extraVars = this->extraPars();
652 
653  // Then store the final fit parameters, and any extra parameters for
654  // the signal model (e.g. fit fractions)
655  for (UInt_t i(0); i<nSigComp_; ++i) {
656  for (UInt_t j(i); j<nSigComp_; ++j) {
657  extraVars.push_back(fitFrac_[i][j]);
658  extraVars.push_back(fitFracEffUnCorr_[i][j]);
659  }
660  }
661 
662  // Store any extra parameters/quantities from the DP model (e.g. K-matrix total fit fractions)
663  std::vector<LauParameter> extraParams = sigDPModel_->getExtraParameters();
664  std::vector<LauParameter>::iterator extraIter;
665  for (extraIter = extraParams.begin(); extraIter != extraParams.end(); ++extraIter) {
666  LauParameter extraParameter = (*extraIter);
667  extraVars.push_back(extraParameter);
668  }
669 
670  // Now add in the DP efficiency value
671  extraVars.push_back(meanEff_);
672 
673  // Also add in the DP rate
674  extraVars.push_back(dpRate_);
675 
676  this->printFitFractions(std::cout);
677  }
678 
679  const LauParameterPList& fitVars = this->fitPars();
680  const LauParameterList& extraVars = this->extraPars();
681  LauFitNtuple* ntuple = this->fitNtuple();
682  ntuple->storeParsAndErrors(fitVars, extraVars);
683 
684  // find out the correlation matrix for the parameters
685  ntuple->storeCorrMatrix(this->iExpt(), this->nll(), this->fitStatus(), this->covarianceMatrix());
686 
687  // Fill the data into ntuple
688  ntuple->updateFitNtuple();
689 
690  // Print out the partial fit fractions, phases and the
691  // averaged efficiency, reweighted by the dynamics (and anything else)
692  if (this->writeLatexTable()) {
693  TString sigOutFileName(tablePrefixName);
694  sigOutFileName += "_"; sigOutFileName += this->iExpt(); sigOutFileName += "Expt.tex";
695  this->writeOutTable(sigOutFileName);
696  }
697 }
698 
699 void LauSimpleFitModel::printFitFractions(std::ostream& output)
700 {
701  // Print out Fit Fractions, total DP rate and mean efficiency
702  for (UInt_t i = 0; i < nSigComp_; i++) {
703  output << "FitFraction for component " << i << " (" << coeffPars_[i]->name() << ") = " << fitFrac_[i][i] << std::endl;
704  }
705  output << "Overall DP rate (integral of matrix element squared) = " << dpRate_ << std::endl;
706  output << "Average efficiency weighted by whole DP dynamics = " << meanEff_ << std::endl;
707 }
708 
709 void LauSimpleFitModel::writeOutTable(const TString& outputFile)
710 {
711  // Write out the results of the fit to a tex-readable table
712  // TODO - need to include the yields in this table
713  std::ofstream fout(outputFile);
714  LauPrint print;
715 
716  std::cout << "INFO in LauSimpleFitModel::writeOutTable : Writing out results of the fit to the tex file " << outputFile << std::endl;
717 
718  if (this->useDP() == kTRUE) {
719  // print the fit coefficients in one table
720  coeffPars_.front()->printTableHeading(fout);
721  for (UInt_t i = 0; i < nSigComp_; i++) {
722  coeffPars_[i]->printTableRow(fout);
723  }
724  fout << "\\hline" << std::endl;
725  fout << "\\end{tabular}" << std::endl << std::endl;
726 
727  // print the fit fractions in another
728  fout << "\\begin{tabular}{|l|c|}" << std::endl;
729  fout << "\\hline" << std::endl;
730  fout << "Component & FitFraction \\\\" << std::endl;
731  fout << "\\hline" << std::endl;
732  Double_t fitFracSum(0.0);
733  for (UInt_t i = 0; i < nSigComp_; i++) {
734  TString resName = coeffPars_[i]->name();
735  resName = resName.ReplaceAll("_", "\\_");
736  fout << resName << " & $";
737  Double_t fitFrac = fitFrac_[i][i].value();
738  fitFracSum += fitFrac;
739  print.printFormat(fout, fitFrac);
740  fout << "$ \\\\" << std::endl;
741  }
742  fout << "\\hline" << std::endl;
743 
744  // Also print out sum of fit fractions
745  fout << "Fit Fraction Sum & $";
746  print.printFormat(fout, fitFracSum);
747  fout << "$ \\\\" << std::endl;
748  fout << "\\hline" << std::endl;
749 
750  fout << "DP rate & $";
751  print.printFormat(fout, dpRate_.value());
752  fout << "$ \\\\" << std::endl;
753  fout << "\\hline" << std::endl;
754 
755  fout << "$< \\varepsilon >$ & $";
756  print.printFormat(fout, meanEff_.value());
757  fout << "$ \\\\" << std::endl;
758  fout << "\\hline" << std::endl;
759  fout << "\\end{tabular}" << std::endl << std::endl;
760  }
761 
762  if (!signalPdfs_.empty()) {
763  fout << "\\begin{tabular}{|l|c|}" << std::endl;
764  fout << "\\hline" << std::endl;
765  if (useSCF_ == kTRUE) {
766  fout << "\\Extra TM Signal PDFs' Parameters: & \\\\" << std::endl;
767  } else {
768  fout << "\\Extra Signal PDFs' Parameters: & \\\\" << std::endl;
769  }
770  this->printFitParameters(signalPdfs_, fout);
771  if (useSCF_ == kTRUE && !scfPdfs_.empty()) {
772  fout << "\\hline" << std::endl;
773  fout << "\\Extra SCF Signal PDFs' Parameters: & \\\\" << std::endl;
774  this->printFitParameters(scfPdfs_, fout);
775  }
776  if (usingBkgnd_ == kTRUE && !bkgndPdfs_.empty()) {
777  fout << "\\hline" << std::endl;
778  fout << "\\Extra Background PDFs' Parameters: & \\\\" << std::endl;
779  for (LauBkgndPdfsList::const_iterator iter = bkgndPdfs_.begin(); iter != bkgndPdfs_.end(); ++iter) {
780  this->printFitParameters(*iter, fout);
781  }
782  }
783  fout << "\\hline \n\\end{tabular}" << std::endl << std::endl;
784  }
785 }
786 
788 {
789  // Update the number of signal events to be total-sum(background events)
790  this->updateSigEvents();
791 
792  // Check whether we want to have randomised initial fit parameters for the signal model
793  if (this->useRandomInitFitPars() == kTRUE) {
794  std::cout << "INFO in LauSimpleFitModel::checkInitFitParams : Setting random parameters for the signal DP model" << std::endl;
795  this->randomiseInitFitPars();
796  }
797 }
798 
800 {
801  // Only randomise those parameters that are not fixed!
802  std::cout << "INFO in LauSimpleFitModel::randomiseInitFitPars : Randomising the initial fit magnitudes and phases of the resonances..." << std::endl;
803 
804  for (UInt_t i = 0; i < nSigComp_; i++) {
805  coeffPars_[i]->randomiseInitValues();
806  }
807 }
808 
810 {
811  // Determine the number of events to generate for each hypothesis
812  // If we're smearing then smear each one individually
813 
814  LauGenInfo nEvtsGen;
815 
816  // Signal
817  Double_t evtWeight(1.0);
818  Int_t nEvts = TMath::FloorNint(signalEvents_->genValue());
819  if ( nEvts < 0 ) {
820  evtWeight = -1.0;
821  nEvts = TMath::Abs( nEvts );
822  }
823  if (this->doPoissonSmearing()) {
824  nEvts = LauRandom::randomFun()->Poisson(nEvts);
825  }
826  nEvtsGen["signal"] = std::make_pair( nEvts, evtWeight );
827 
828  // Backgrounds
829  const UInt_t nBkgnds = this->nBkgndClasses();
830  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
831  const TString& bkgndClass = this->bkgndClassName(bkgndID);
832  const LauParameter* evtsPar = bkgndEvents_[bkgndID];
833  evtWeight = 1.0;
834  nEvts = TMath::FloorNint( evtsPar->genValue() );
835  if ( nEvts < 0 ) {
836  evtWeight = -1.0;
837  nEvts = TMath::Abs( nEvts );
838  }
839  if (this->doPoissonSmearing()) {
840  nEvts = LauRandom::randomFun()->Poisson(nEvts);
841  }
842  nEvtsGen[bkgndClass] = std::make_pair( nEvts, evtWeight );
843  }
844 
845  return nEvtsGen;
846 }
847 
849 {
850  // Routine to generate toy Monte Carlo events according to the various models we have defined.
851 
852  // Determine the number of events to generate for each hypothesis
853  LauGenInfo nEvts = this->eventsToGenerate();
854 
855  Bool_t genOK(kTRUE);
856  Int_t evtNum(0);
857 
858  const UInt_t nBkgnds = this->nBkgndClasses();
859  std::vector<TString> bkgndClassNames(nBkgnds);
860  std::vector<TString> bkgndClassNamesGen(nBkgnds);
861  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
862  TString name( this->bkgndClassName(iBkgnd) );
863  bkgndClassNames[iBkgnd] = name;
864  bkgndClassNamesGen[iBkgnd] = "gen"+name;
865  }
866 
867  const Bool_t storeSCFTruthInfo = ( useSCF_ || ( this->enableEmbedding() && signalTree_ != 0 && signalTree_->haveBranch("mcMatch") ) );
868 
869  // Loop over the hypotheses and generate the requested number of events for each one
870  for (LauGenInfo::const_iterator iter = nEvts.begin(); iter != nEvts.end(); ++iter) {
871 
872  const TString& type(iter->first);
873  Int_t nEvtsGen( iter->second.first );
874  Double_t evtWeight( iter->second.second );
875 
876  for (Int_t iEvt(0); iEvt<nEvtsGen; ++iEvt) {
877 
878  this->setGenNtupleDoubleBranchValue( "evtWeight", evtWeight );
879  // Add efficiency information
880  this->setGenNtupleDoubleBranchValue( "efficiency", 1 );
881 
882  if (type == "signal") {
883  this->setGenNtupleIntegerBranchValue("genSig",1);
884  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
885  this->setGenNtupleIntegerBranchValue( bkgndClassNamesGen[iBkgnd], 0 );
886  }
887  genOK = this->generateSignalEvent();
888  this->setGenNtupleDoubleBranchValue( "efficiency", sigDPModel_->getEvtEff() );
889  } else {
890  this->setGenNtupleIntegerBranchValue("genSig",0);
891  if ( storeSCFTruthInfo ) {
892  this->setGenNtupleIntegerBranchValue("genTMSig",0);
893  this->setGenNtupleIntegerBranchValue("genSCFSig",0);
894  }
895  UInt_t bkgndID(0);
896  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
897  Int_t gen(0);
898  if ( bkgndClassNames[iBkgnd] == type ) {
899  gen = 1;
900  bkgndID = iBkgnd;
901  }
902  this->setGenNtupleIntegerBranchValue( bkgndClassNamesGen[iBkgnd], gen );
903  }
904  genOK = this->generateBkgndEvent(bkgndID);
905  }
906 
907  if (!genOK) {
908  // If there was a problem with the generation then break out and return.
909  // The problem model will have adjusted itself so that all should be OK next time.
910  break;
911  }
912 
913  if (this->useDP() == kTRUE) {
914  this->setDPBranchValues();
915  }
916 
917  // Store the event number (within this experiment)
918  this->setGenNtupleIntegerBranchValue("iEvtWithinExpt",evtNum);
919  // and then increment it
920  ++evtNum;
921 
922  this->fillGenNtupleBranches();
923  if (iEvt%500 == 0) {
924  std::cout << "INFO in LauSimpleFitModel::genExpt : Generated event number " << iEvt << " out of " << nEvtsGen << " " << type << " events." << std::endl;
925  }
926  }
927 
928  if (!genOK) {
929  break;
930  }
931  }
932 
933  if (this->useDP() && genOK) {
934 
935  sigDPModel_->checkToyMC(kTRUE,kTRUE);
936 
937  // Get the fit fractions if they're to be written into the latex table
938  if (this->writeLatexTable()) {
940  if (fitFrac.size() != nSigComp_) {
941  std::cerr << "ERROR in LauSimpleFitModel::genExpt : Fit Fraction array of unexpected dimension: " << fitFrac.size() << std::endl;
942  gSystem->Exit(EXIT_FAILURE);
943  }
944  for (UInt_t i(0); i<nSigComp_; ++i) {
945  if (fitFrac[i].size() != nSigComp_) {
946  std::cerr << "ERROR in LauSimpleFitModel::genExpt : Fit Fraction array of unexpected dimension: " << fitFrac[i].size() << std::endl;
947  gSystem->Exit(EXIT_FAILURE);
948  }
949  }
950 
951  for (UInt_t i(0); i<nSigComp_; ++i) {
952  for (UInt_t j = i; j < nSigComp_; j++) {
953  fitFrac_[i][j].value(fitFrac[i][j].value());
954  }
955  }
958  }
959  }
960 
961  // If we're reusing embedded events or if the generation is being
962  // reset then clear the lists of used events
963  if (reuseSignal_ || !genOK) {
964  if (signalTree_) {
966  }
967  }
968  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
969  LauEmbeddedData* data = bkgndTree_[bkgndID];
970  if (reuseBkgnd_[bkgndID] || !genOK) {
971  if (data) {
972  data->clearUsedList();
973  }
974  }
975  }
976 
977  return genOK;
978 }
979 
981 {
982  // Generate signal event
983  Bool_t genOK(kTRUE);
984  Bool_t genSCF(kFALSE);
985 
986  if (this->useDP()) {
987  if (this->enableEmbedding() && signalTree_) {
988  if (useReweighting_) {
989  // Select a (random) event from the generated data. Then store the
990  // reconstructed DP co-ords, together with other pdf information,
991  // as the embedded data.
993  } else {
994 
996  }
997  if (signalTree_->haveBranch("mcMatch")) {
998  Int_t match = static_cast<Int_t>(signalTree_->getValue("mcMatch"));
999  if (match) {
1000  this->setGenNtupleIntegerBranchValue("genTMSig",1);
1001  this->setGenNtupleIntegerBranchValue("genSCFSig",0);
1002  genSCF = kFALSE;
1003  } else {
1004  this->setGenNtupleIntegerBranchValue("genTMSig",0);
1005  this->setGenNtupleIntegerBranchValue("genSCFSig",1);
1006  genSCF = kTRUE;
1007  }
1008  }
1009  } else {
1010  genOK = sigDPModel_->generate();
1011  if ( genOK && useSCF_ ) {
1012  Double_t frac(0.0);
1013  if ( useSCFHist_ ) {
1015  } else {
1016  frac = scfFrac_.genValue();
1017  }
1018  if ( frac < LauRandom::randomFun()->Rndm() ) {
1019  this->setGenNtupleIntegerBranchValue("genTMSig",1);
1020  this->setGenNtupleIntegerBranchValue("genSCFSig",0);
1021  genSCF = kFALSE;
1022  } else {
1023  this->setGenNtupleIntegerBranchValue("genTMSig",0);
1024  this->setGenNtupleIntegerBranchValue("genSCFSig",1);
1025  genSCF = kTRUE;
1026 
1027  // Optionally smear the DP position
1028  // of the SCF event
1029  if ( scfMap_ != 0 ) {
1030  Double_t xCoord(0.0), yCoord(0.0);
1031  if ( kinematics_->squareDP() ) {
1032  xCoord = kinematics_->getmPrime();
1033  yCoord = kinematics_->getThetaPrime();
1034  } else {
1035  xCoord = kinematics_->getm13Sq();
1036  yCoord = kinematics_->getm23Sq();
1037  }
1038 
1039  // Find the bin number where this event is generated
1040  Int_t binNo = scfMap_->binNumber( xCoord, yCoord );
1041 
1042  // Retrieve the migration histogram
1043  TH2* histo = scfMap_->trueHist( binNo );
1044 
1045  LauAbsEffModel * effModel = sigDPModel_->getEffModel();
1046  do {
1047  // Get a random point from the histogram
1048  histo->GetRandom2( xCoord, yCoord );
1049 
1050  // Update the kinematics
1051  if ( kinematics_->squareDP() ) {
1052  kinematics_->updateSqDPKinematics( xCoord, yCoord );
1053  } else {
1054  kinematics_->updateKinematics( xCoord, yCoord );
1055  }
1056  } while ( ! effModel->passVeto( kinematics_ ) );
1057  }
1058  }
1059  }
1060  }
1061  } else {
1062  if (this->enableEmbedding() && signalTree_) {
1064  if (signalTree_->haveBranch("mcMatch")) {
1065  Int_t match = static_cast<Int_t>(signalTree_->getValue("mcMatch"));
1066  if (match) {
1067  this->setGenNtupleIntegerBranchValue("genTMSig",1);
1068  this->setGenNtupleIntegerBranchValue("genSCFSig",0);
1069  genSCF = kFALSE;
1070  } else {
1071  this->setGenNtupleIntegerBranchValue("genTMSig",0);
1072  this->setGenNtupleIntegerBranchValue("genSCFSig",1);
1073  genSCF = kTRUE;
1074  }
1075  }
1076  } else if ( useSCF_ ) {
1077  Double_t frac = scfFrac_.genValue();
1078  if ( frac < LauRandom::randomFun()->Rndm() ) {
1079  this->setGenNtupleIntegerBranchValue("genTMSig",1);
1080  this->setGenNtupleIntegerBranchValue("genSCFSig",0);
1081  genSCF = kFALSE;
1082  } else {
1083  this->setGenNtupleIntegerBranchValue("genTMSig",0);
1084  this->setGenNtupleIntegerBranchValue("genSCFSig",1);
1085  genSCF = kTRUE;
1086  }
1087  }
1088  }
1089  if (genOK) {
1090  if ( useSCF_ ) {
1091  if ( genSCF ) {
1093  } else {
1095  }
1096  } else {
1098  }
1099  }
1100  // Check for problems with the embedding
1101  if (this->enableEmbedding() && signalTree_ && (signalTree_->nEvents() == signalTree_->nUsedEvents())) {
1102  std::cerr << "WARNING in LauSimpleFitModel::generateSignalEvent : Source of embedded signal events used up, clearing the list of used events." << std::endl;
1104  }
1105  return genOK;
1106 }
1107 
1109 {
1110  // Generate background event
1111  Bool_t genOK(kTRUE);
1112 
1113  LauAbsBkgndDPModel* model = bkgndDPModels_[bkgndID];
1114  LauEmbeddedData* embeddedData(0);
1115  if (this->enableEmbedding()) {
1116  embeddedData = bkgndTree_[bkgndID];
1117  }
1118  LauPdfList* extraPdfs = &bkgndPdfs_[bkgndID];
1119  if (this->useDP()) {
1120  if (embeddedData) {
1121  embeddedData->getEmbeddedEvent(kinematics_);
1122  } else {
1123  if (model == 0) {
1124  std::cerr << "ERROR in LauSimpleFitModel::generateBkgndEvent : Can't find the DP model for background class \"" << this->bkgndClassName(bkgndID) << "\"." << std::endl;
1125  gSystem->Exit(EXIT_FAILURE);
1126  }
1127  genOK = model->generate();
1128  }
1129  } else {
1130  if (embeddedData) {
1131  embeddedData->getEmbeddedEvent(0);
1132  }
1133  }
1134  if (genOK) {
1135  this->generateExtraPdfValues(extraPdfs, embeddedData);
1136  }
1137  // Check for problems with the embedding
1138  if (embeddedData && (embeddedData->nEvents() == embeddedData->nUsedEvents())) {
1139  std::cerr << "WARNING in LauSimpleFitModel::generateBkgndEvent : Source of embedded " << this->bkgndClassName(bkgndID) << " events used up, clearing the list of used events." << std::endl;
1140  embeddedData->clearUsedList();
1141  }
1142  return genOK;
1143 }
1144 
1146 {
1147  // Setup the required ntuple branches
1148  this->addGenNtupleDoubleBranch("evtWeight");
1149  this->addGenNtupleIntegerBranch("genSig");
1150  this->addGenNtupleDoubleBranch("efficiency");
1151  if ( useSCF_ || ( this->enableEmbedding() && signalTree_ != 0 && signalTree_->haveBranch("mcMatch") ) ) {
1152  this->addGenNtupleIntegerBranch("genTMSig");
1153  this->addGenNtupleIntegerBranch("genSCFSig");
1154  }
1155  const UInt_t nBkgnds = this->nBkgndClasses();
1156  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
1157  TString name( this->bkgndClassName(iBkgnd) );
1158  name.Prepend("gen");
1159  this->addGenNtupleIntegerBranch(name);
1160  }
1161  if (this->useDP() == kTRUE) {
1162  this->addGenNtupleDoubleBranch("m12");
1163  this->addGenNtupleDoubleBranch("m23");
1164  this->addGenNtupleDoubleBranch("m13");
1165  this->addGenNtupleDoubleBranch("m12Sq");
1166  this->addGenNtupleDoubleBranch("m23Sq");
1167  this->addGenNtupleDoubleBranch("m13Sq");
1168  this->addGenNtupleDoubleBranch("cosHel12");
1169  this->addGenNtupleDoubleBranch("cosHel23");
1170  this->addGenNtupleDoubleBranch("cosHel13");
1171  if (kinematics_->squareDP()) {
1172  this->addGenNtupleDoubleBranch("mPrime");
1173  this->addGenNtupleDoubleBranch("thPrime");
1174  }
1175  }
1176  for (LauPdfList::const_iterator pdf_iter = signalPdfs_.begin(); pdf_iter != signalPdfs_.end(); ++pdf_iter) {
1177  std::vector<TString> varNames = (*pdf_iter)->varNames();
1178  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
1179  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
1180  this->addGenNtupleDoubleBranch( (*var_iter) );
1181  }
1182  }
1183  }
1184 }
1185 
1187 {
1188  // Store all the DP information
1195  this->setGenNtupleDoubleBranchValue("cosHel12", kinematics_->getc12());
1196  this->setGenNtupleDoubleBranchValue("cosHel23", kinematics_->getc23());
1197  this->setGenNtupleDoubleBranchValue("cosHel13", kinematics_->getc13());
1198  if (kinematics_->squareDP()) {
1201  }
1202 }
1203 
1205 {
1206  // Generate from the extra PDFs
1207  if (extraPdfs) {
1208  for (LauPdfList::iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end(); ++pdf_iter) {
1209  LauFitData genValues;
1210  if (embeddedData) {
1211  genValues = embeddedData->getValues( (*pdf_iter)->varNames() );
1212  } else {
1213  genValues = (*pdf_iter)->generate(kinematics_);
1214  }
1215  for ( LauFitData::const_iterator var_iter = genValues.begin(); var_iter != genValues.end(); ++var_iter ) {
1216  TString varName = var_iter->first;
1217  if ( varName != "m13Sq" && varName != "m23Sq" ) {
1218  Double_t value = var_iter->second;
1219  this->setGenNtupleDoubleBranchValue(varName,value);
1220  }
1221  }
1222  }
1223  }
1224 }
1225 
1227 {
1228  // Update the total normalisation for the signal likelihood
1229  if (this->useDP() == kTRUE) {
1230  this->updateCoeffs();
1232  }
1233 
1234  // Update the signal events from the background events if not doing an extended fit
1235  if ( !this->doEMLFit() && !signalEvents_->fixed() ) {
1236  this->updateSigEvents();
1237  }
1238 }
1239 
1241 {
1242  // The background parameters will have been set from Minuit.
1243  // We need to update the signal events using these.
1244  Double_t nTotEvts = this->eventsPerExpt();
1245 
1246  signalEvents_->range(-2.0*nTotEvts,2.0*nTotEvts);
1247  for (LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
1248  (*iter)->range(-2.0*nTotEvts,2.0*nTotEvts);
1249  }
1250 
1251  if ( signalEvents_->fixed() ) {
1252  return;
1253  }
1254 
1255  // Subtract background events (if any) from signal.
1256  Double_t signalEvents = nTotEvts;
1257  if ( usingBkgnd_ == kTRUE ) {
1258  for (LauBkgndYieldList::const_iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
1259  signalEvents -= (*iter)->value();
1260  }
1261  }
1262  signalEvents_->value(signalEvents);
1263 }
1264 
1266 {
1267  // Fill the internal data trees of the signal and backgrond models.
1268 
1269  LauFitDataTree* inputFitData = this->fitData();
1270 
1271  // First the Dalitz plot variables (m_ij^2)
1272  if (this->useDP() == kTRUE) {
1273 
1274  // need to append SCF smearing bins before caching DP amplitudes
1275  if ( scfMap_ != 0 ) {
1276  this->appendBinCentres( inputFitData );
1277  }
1278  sigDPModel_->fillDataTree(*inputFitData);
1279 
1280  if (usingBkgnd_ == kTRUE) {
1281  for (LauBkgndDPModelList::iterator iter = bkgndDPModels_.begin(); iter != bkgndDPModels_.end(); ++iter) {
1282  (*iter)->fillDataTree(*inputFitData);
1283  }
1284  }
1285  }
1286 
1287  // ...and then the extra PDFs
1288  this->cacheInfo(signalPdfs_, *inputFitData);
1289  this->cacheInfo(scfPdfs_, *inputFitData);
1290  for (LauBkgndPdfsList::iterator iter = bkgndPdfs_.begin(); iter != bkgndPdfs_.end(); ++iter) {
1291  this->cacheInfo((*iter), *inputFitData);
1292  }
1293 
1294  // and finally the SCF fractions and jacobians
1295  if ( useSCF_ && useSCFHist_ ) {
1296  if ( !inputFitData->haveBranch( "m13Sq" ) || !inputFitData->haveBranch( "m23Sq" ) ) {
1297  std::cerr << "ERROR in LauSimpleFitModel::cacheInputFitVars : Input data does not contain DP branches and so can't cache the SCF fraction." << std::endl;
1298  gSystem->Exit(EXIT_FAILURE);
1299  }
1300 
1301  UInt_t nEvents = inputFitData->nEvents();
1302  recoSCFFracs_.clear();
1303  recoSCFFracs_.reserve( nEvents );
1304  if ( kinematics_->squareDP() ) {
1305  recoJacobians_.clear();
1306  recoJacobians_.reserve( nEvents );
1307  }
1308  for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) {
1309  const LauFitData& dataValues = inputFitData->getData(iEvt);
1310  LauFitData::const_iterator m13_iter = dataValues.find("m13Sq");
1311  LauFitData::const_iterator m23_iter = dataValues.find("m23Sq");
1312  kinematics_->updateKinematics( m13_iter->second, m23_iter->second );
1313  Double_t scfFrac = scfFracHist_->calcEfficiency( kinematics_ );
1314  recoSCFFracs_.push_back( scfFrac );
1315  if ( kinematics_->squareDP() ) {
1317  }
1318  }
1319  }
1320 }
1321 
1323 {
1324  // We'll be caching the DP amplitudes and efficiencies of the centres of the true bins.
1325  // To do so, we attach some fake points at the end of inputData, the number of the entry
1326  // minus the total number of events corresponding to the number of the histogram for that
1327  // given true bin in the LauScfMap object. (What this means is that when Laura is provided with
1328  // the LauScfMap object by the user, it's the latter who has to make sure that it contains the
1329  // right number of histograms and in exactly the right order!)
1330 
1331  // Get the x and y co-ordinates of the bin centres
1332  std::vector<Double_t> binCentresXCoords;
1333  std::vector<Double_t> binCentresYCoords;
1334  scfMap_->listBinCentres(binCentresXCoords, binCentresYCoords);
1335 
1336  // The SCF histograms could be in square Dalitz plot histograms.
1337  // The dynamics takes normal Dalitz plot coords, so we might have to convert them back.
1338  Bool_t sqDP = kinematics_->squareDP();
1339  UInt_t nBins = binCentresXCoords.size();
1340  fakeSCFFracs_.clear();
1341  fakeSCFFracs_.reserve( nBins );
1342  if ( sqDP ) {
1343  fakeJacobians_.clear();
1344  fakeJacobians_.reserve( nBins );
1345  }
1346 
1347  for (UInt_t iBin = 0; iBin < nBins; ++iBin) {
1348 
1349  if ( sqDP ) {
1350  kinematics_->updateSqDPKinematics(binCentresXCoords[iBin],binCentresYCoords[iBin]);
1351  binCentresXCoords[iBin] = kinematics_->getm13Sq();
1352  binCentresYCoords[iBin] = kinematics_->getm23Sq();
1354  } else {
1355  kinematics_->updateKinematics(binCentresXCoords[iBin],binCentresYCoords[iBin]);
1356  }
1357 
1359  }
1360 
1361  // Set up inputFitVars_ object to hold the fake events
1362  inputData->appendFakePoints(binCentresXCoords,binCentresYCoords);
1363 }
1364 
1366 {
1367  // Get the DP likelihood for signal and backgrounds
1368  this->getEvtDPLikelihood(iEvt);
1369 
1370  // Get the combined extra PDFs likelihood for signal and backgrounds
1371  this->getEvtExtraLikelihoods(iEvt);
1372 
1373  // If appropriate, combine the TM and SCF likelihoods
1374  Double_t sigLike = sigDPLike_ * sigExtraLike_;
1375  if ( useSCF_ ) {
1376  Double_t scfFrac(0.0);
1377  if (useSCFHist_) {
1378  scfFrac = recoSCFFracs_[iEvt];
1379  } else {
1380  scfFrac = scfFrac_.value();
1381  }
1382  sigLike *= (1.0 - scfFrac);
1383  if ( (scfMap_ != 0) && (this->useDP() == kTRUE) ) {
1384  // if we're smearing the SCF DP PDF then the SCF frac
1385  // is already included in the SCF DP likelihood
1386  sigLike += (scfDPLike_ * scfExtraLike_);
1387  } else {
1388  sigLike += (scfFrac * scfDPLike_ * scfExtraLike_);
1389  }
1390  }
1391 
1392  // Construct the total event likelihood
1393  Double_t likelihood = signalEvents_->value() * sigLike;
1394  const UInt_t nBkgnds = this->nBkgndClasses();
1395  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
1396  likelihood += (bkgndEvents_[bkgndID]->value() * bkgndDPLike_[bkgndID] * bkgndExtraLike_[bkgndID]);
1397  }
1398 
1399  return likelihood;
1400 }
1401 
1403 {
1404  Double_t eventSum(0.0);
1405  eventSum += signalEvents_->value();
1406  for (LauBkgndYieldList::const_iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
1407  eventSum += (*iter)->value();
1408  }
1409  return eventSum;
1410 }
1411 
1413 {
1414  // Function to return the signal and background likelihoods for the
1415  // Dalitz plot for the given event evtNo.
1416 
1417  if (this->useDP() == kTRUE) {
1420 
1421  if ( useSCF_ == kTRUE ) {
1422  if ( scfMap_ == 0 ) {
1423  // we're not smearing the SCF DP position
1424  // so the likelihood is the same as the TM
1426  } else {
1427  // calculate the smeared SCF DP likelihood
1428  scfDPLike_ = this->getEvtSCFDPLikelihood(iEvt);
1429  }
1430  }
1431 
1432  const UInt_t nBkgnds = this->nBkgndClasses();
1433  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
1434  if (usingBkgnd_ == kTRUE) {
1435  bkgndDPLike_[bkgndID] = bkgndDPModels_[bkgndID]->getLikelihood(iEvt);
1436  } else {
1437  bkgndDPLike_[bkgndID] = 0.0;
1438  }
1439  }
1440  } else {
1441  // There's always going to be a term in the likelihood for the
1442  // signal, so we'd better not zero it.
1443  sigDPLike_ = 1.0;
1444  scfDPLike_ = 1.0;
1445 
1446  const UInt_t nBkgnds = this->nBkgndClasses();
1447  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
1448  if (usingBkgnd_ == kTRUE) {
1449  bkgndDPLike_[bkgndID] = 1.0;
1450  } else {
1451  bkgndDPLike_[bkgndID] = 0.0;
1452  }
1453  }
1454  }
1455 }
1456 
1458 {
1459  Double_t scfDPLike(0.0);
1460 
1461  Double_t recoJacobian(1.0);
1462  Double_t xCoord(0.0);
1463  Double_t yCoord(0.0);
1464  Bool_t squareDP = kinematics_->squareDP();
1465  if ( squareDP ) {
1466  xCoord = sigDPModel_->getEvtmPrime();
1467  yCoord = sigDPModel_->getEvtthPrime();
1468  recoJacobian = recoJacobians_[iEvt];
1469  } else {
1470  xCoord = sigDPModel_->getEvtm13Sq();
1471  yCoord = sigDPModel_->getEvtm23Sq();
1472  }
1473 
1474  // Find the bin that our reco event falls in
1475  Int_t recoBin = scfMap_->binNumber( xCoord, yCoord );
1476 
1477  // Find out which true Bins contribute to the given reco bin
1478  const std::vector<Int_t>* trueBins = scfMap_->trueBins(recoBin);
1479 
1480  Int_t nDataEvents = this->eventsPerExpt();
1481 
1482  // Loop over the true bins
1483  for (std::vector<Int_t>::const_iterator iter = trueBins->begin(); iter != trueBins->end(); ++iter)
1484  {
1485  Int_t trueBin = (*iter);
1486 
1487  // prob of a true event in the given true bin migrating to the reco bin
1488  Double_t pRecoGivenTrue = scfMap_->prob( recoBin, trueBin );
1489 
1490  // We've cached the DP amplitudes and the efficiency for the
1491  // true bin centres, just after the data points
1492  sigDPModel_->calcLikelihoodInfo( nDataEvents + trueBin );
1493  Double_t pTrue = sigDPModel_->getEvtLikelihood();
1494 
1495  // Get the cached SCF fraction (and jacobian if we're using the square DP)
1496  Double_t scfFraction = fakeSCFFracs_[ trueBin ];
1497  Double_t jacobian(1.0);
1498  if ( squareDP ) {
1499  jacobian = fakeJacobians_[ trueBin ];
1500  }
1501 
1502  scfDPLike += pTrue * jacobian * scfFraction * pRecoGivenTrue;
1503 
1504  }
1505 
1506  // Divide by the reco jacobian
1507  scfDPLike /= recoJacobian;
1508 
1509  return scfDPLike;
1510 }
1511 
1513 {
1514  // Function to return the signal and background likelihoods for the
1515  // extra variables for the given event evtNo.
1516 
1517  sigExtraLike_ = 1.0;
1518 
1519  for (LauPdfList::iterator iter = signalPdfs_.begin(); iter != signalPdfs_.end(); ++iter) {
1520  (*iter)->calcLikelihoodInfo(iEvt);
1521  sigExtraLike_ *= (*iter)->getLikelihood();
1522  }
1523  if (useSCF_) {
1524  scfExtraLike_ = 1.0;
1525  for (LauPdfList::iterator iter = scfPdfs_.begin(); iter != scfPdfs_.end(); ++iter) {
1526  (*iter)->calcLikelihoodInfo(iEvt);
1527  scfExtraLike_ *= (*iter)->getLikelihood();
1528  }
1529  }
1530 
1531  const UInt_t nBkgnds = this->nBkgndClasses();
1532  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
1533  if (usingBkgnd_) {
1534  bkgndExtraLike_[bkgndID] = 1.0;
1535  LauPdfList& pdfList = bkgndPdfs_[bkgndID];
1536  for (LauPdfList::iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) {
1537  (*pdf_iter)->calcLikelihoodInfo(iEvt);
1538  bkgndExtraLike_[bkgndID] *= (*pdf_iter)->getLikelihood();
1539  }
1540  } else {
1541  bkgndExtraLike_[bkgndID] = 0.0;
1542  }
1543  }
1544 }
1545 
1547 {
1548  coeffs_.clear();
1549  coeffs_.reserve(nSigComp_);
1550  for (UInt_t i = 0; i < nSigComp_; i++) {
1551  coeffs_.push_back(coeffPars_[i]->particleCoeff());
1552  }
1553 }
1554 
1556 {
1557  // add branches for storing the experiment number and the number of
1558  // the event within the current experiment
1559  this->addSPlotNtupleIntegerBranch("iExpt");
1560  this->addSPlotNtupleIntegerBranch("iEvtWithinExpt");
1561 
1562  // Store the efficiency of the event (for inclusive BF calculations).
1563  if (this->storeDPEff()) {
1564  this->addSPlotNtupleDoubleBranch("efficiency");
1565  if ( sigDPModel_->usingScfModel() ) {
1566  this->addSPlotNtupleDoubleBranch("scffraction");
1567  }
1568  }
1569 
1570  // Store the total event likelihood for each species.
1571  if (useSCF_) {
1572  this->addSPlotNtupleDoubleBranch("sigTMTotalLike");
1573  this->addSPlotNtupleDoubleBranch("sigSCFTotalLike");
1574  this->addSPlotNtupleDoubleBranch("sigSCFFrac");
1575  } else {
1576  this->addSPlotNtupleDoubleBranch("sigTotalLike");
1577  }
1578  if (usingBkgnd_) {
1579  const UInt_t nBkgnds = this->nBkgndClasses();
1580  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
1581  TString name( this->bkgndClassName(iBkgnd) );
1582  name += "TotalLike";
1583  this->addSPlotNtupleDoubleBranch(name);
1584  }
1585  }
1586 
1587  // Store the DP likelihoods
1588  if (this->useDP()) {
1589  if (useSCF_) {
1590  this->addSPlotNtupleDoubleBranch("sigTMDPLike");
1591  this->addSPlotNtupleDoubleBranch("sigSCFDPLike");
1592  } else {
1593  this->addSPlotNtupleDoubleBranch("sigDPLike");
1594  }
1595  if (usingBkgnd_) {
1596  const UInt_t nBkgnds = this->nBkgndClasses();
1597  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
1598  TString name( this->bkgndClassName(iBkgnd) );
1599  name += "DPLike";
1600  this->addSPlotNtupleDoubleBranch(name);
1601  }
1602  }
1603  }
1604 
1605  // Store the likelihoods for each extra PDF
1606  if (useSCF_) {
1607  this->addSPlotNtupleBranches(&signalPdfs_, "sigTM");
1608  this->addSPlotNtupleBranches(&scfPdfs_, "sigSCF");
1609  } else {
1610  this->addSPlotNtupleBranches(&signalPdfs_, "sig");
1611  }
1612  if (usingBkgnd_) {
1613  const UInt_t nBkgnds = this->nBkgndClasses();
1614  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
1615  const TString& bkgndClass = this->bkgndClassName(iBkgnd);
1616  const LauPdfList* pdfList = &(bkgndPdfs_[iBkgnd]);
1617  this->addSPlotNtupleBranches(pdfList, bkgndClass);
1618  }
1619  }
1620 }
1621 
1622 void LauSimpleFitModel::addSPlotNtupleBranches(const LauPdfList* extraPdfs, const TString& prefix)
1623 {
1624  if (extraPdfs) {
1625  // Loop through each of the PDFs
1626  for (LauPdfList::const_iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end(); ++pdf_iter) {
1627 
1628  // Count the number of input variables that are not
1629  // DP variables (used in the case where there is DP
1630  // dependence for e.g. MVA)
1631  UInt_t nVars(0);
1632  std::vector<TString> varNames = (*pdf_iter)->varNames();
1633  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
1634  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
1635  ++nVars;
1636  }
1637  }
1638 
1639  if ( nVars == 1 ) {
1640  // If the PDF only has one variable then
1641  // simply add one branch for that variable
1642  TString varName = (*pdf_iter)->varName();
1643  TString name(prefix);
1644  name += varName;
1645  name += "Like";
1646  this->addSPlotNtupleDoubleBranch(name);
1647  } else if ( nVars == 2 ) {
1648  // If the PDF has two variables then we
1649  // need a branch for them both together and
1650  // branches for each
1651  TString allVars("");
1652  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
1653  allVars += (*var_iter);
1654  TString name(prefix);
1655  name += (*var_iter);
1656  name += "Like";
1657  this->addSPlotNtupleDoubleBranch(name);
1658  }
1659  TString name(prefix);
1660  name += allVars;
1661  name += "Like";
1662  this->addSPlotNtupleDoubleBranch(name);
1663  } else {
1664  std::cerr << "WARNING in LauSimpleFitModel::addSPlotNtupleBranches : Can't yet deal with 3D PDFs." << std::endl;
1665  }
1666  }
1667  }
1668 }
1669 
1670 Double_t LauSimpleFitModel::setSPlotNtupleBranchValues(LauPdfList* extraPdfs, const TString& prefix, UInt_t iEvt)
1671 {
1672  // Store the PDF value for each variable in the list
1673  Double_t totalLike(1.0);
1674  Double_t extraLike(0.0);
1675  if (extraPdfs) {
1676  for (LauPdfList::iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end(); ++pdf_iter) {
1677 
1678  // calculate the likelihood for this event
1679  (*pdf_iter)->calcLikelihoodInfo(iEvt);
1680  extraLike = (*pdf_iter)->getLikelihood();
1681  totalLike *= extraLike;
1682 
1683  // Count the number of input variables that are not
1684  // DP variables (used in the case where there is DP
1685  // dependence for e.g. MVA)
1686  UInt_t nVars(0);
1687  std::vector<TString> varNames = (*pdf_iter)->varNames();
1688  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
1689  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
1690  ++nVars;
1691  }
1692  }
1693 
1694  if ( nVars == 1 ) {
1695  // If the PDF only has one variable then
1696  // simply store the value for that variable
1697  TString varName = (*pdf_iter)->varName();
1698  TString name(prefix);
1699  name += varName;
1700  name += "Like";
1701  this->setSPlotNtupleDoubleBranchValue(name, extraLike);
1702  } else if ( nVars == 2 ) {
1703  // If the PDF has two variables then we
1704  // store the value for them both together
1705  // and for each on their own
1706  TString allVars("");
1707  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
1708  allVars += (*var_iter);
1709  TString name(prefix);
1710  name += (*var_iter);
1711  name += "Like";
1712  Double_t indivLike = (*pdf_iter)->getLikelihood( (*var_iter) );
1713  this->setSPlotNtupleDoubleBranchValue(name, indivLike);
1714  }
1715  TString name(prefix);
1716  name += allVars;
1717  name += "Like";
1718  this->setSPlotNtupleDoubleBranchValue(name, extraLike);
1719  } else {
1720  std::cerr << "WARNING in LauSimpleFitModel::setSPlotNtupleBranchValues : Can't yet deal with 3D PDFs." << std::endl;
1721  }
1722  }
1723  }
1724  return totalLike;
1725 }
1726 
1728 {
1729  LauSPlot::NameSet nameSet;
1730  if (this->useDP()) {
1731  nameSet.insert("DP");
1732  }
1733  // Loop through all the signal PDFs
1734  for (LauPdfList::const_iterator pdf_iter = signalPdfs_.begin(); pdf_iter != signalPdfs_.end(); ++pdf_iter) {
1735  // Loop over the variables involved in each PDF
1736  std::vector<TString> varNames = (*pdf_iter)->varNames();
1737  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
1738  // If they are not DP coordinates then add them
1739  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
1740  nameSet.insert( (*var_iter) );
1741  }
1742  }
1743  }
1744  return nameSet;
1745 }
1746 
1748 {
1749  LauSPlot::NumbMap numbMap;
1750  if (!signalEvents_->fixed() && this->doEMLFit()) {
1751  numbMap["sig"] = signalEvents_->genValue();
1752  }
1753  if ( usingBkgnd_ ) {
1754  const UInt_t nBkgnds = this->nBkgndClasses();
1755  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
1756  const TString& bkgndClass = this->bkgndClassName(iBkgnd);
1757  const LauParameter* par = bkgndEvents_[iBkgnd];
1758  if (!par->fixed()) {
1759  numbMap[bkgndClass] = par->genValue();
1760  }
1761  }
1762  }
1763  return numbMap;
1764 }
1765 
1767 {
1768  LauSPlot::NumbMap numbMap;
1769  if (signalEvents_->fixed() && this->doEMLFit()) {
1770  numbMap["sig"] = signalEvents_->genValue();
1771  }
1772  if ( usingBkgnd_ ) {
1773  const UInt_t nBkgnds = this->nBkgndClasses();
1774  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
1775  const TString& bkgndClass = this->bkgndClassName(iBkgnd);
1776  const LauParameter* par = bkgndEvents_[iBkgnd];
1777  if (par->fixed()) {
1778  numbMap[bkgndClass] = par->genValue();
1779  }
1780  }
1781  }
1782  return numbMap;
1783 }
1784 
1786 {
1787  LauSPlot::TwoDMap twodimMap;
1788 
1789  for (LauPdfList::const_iterator pdf_iter = signalPdfs_.begin(); pdf_iter != signalPdfs_.end(); ++pdf_iter) {
1790  // Count the number of input variables that are not DP variables
1791  UInt_t nVars(0);
1792  std::vector<TString> varNames = (*pdf_iter)->varNames();
1793  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
1794  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
1795  ++nVars;
1796  }
1797  }
1798  if ( nVars == 2 ) {
1799  if (useSCF_) {
1800  twodimMap.insert( std::make_pair( "sigTM", std::make_pair( varNames[0], varNames[1] ) ) );
1801  } else {
1802  twodimMap.insert( std::make_pair( "sig", std::make_pair( varNames[0], varNames[1] ) ) );
1803  }
1804  }
1805  }
1806 
1807  if ( useSCF_ ) {
1808  for (LauPdfList::const_iterator pdf_iter = scfPdfs_.begin(); pdf_iter != scfPdfs_.end(); ++pdf_iter) {
1809  // Count the number of input variables that are not DP variables
1810  UInt_t nVars(0);
1811  std::vector<TString> varNames = (*pdf_iter)->varNames();
1812  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
1813  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
1814  ++nVars;
1815  }
1816  }
1817  if ( nVars == 2 ) {
1818  twodimMap.insert( std::make_pair( "sigSCF", std::make_pair( varNames[0], varNames[1] ) ) );
1819  }
1820  }
1821  }
1822 
1823  if (usingBkgnd_) {
1824  const UInt_t nBkgnds = this->nBkgndClasses();
1825  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
1826  const TString& bkgndClass = this->bkgndClassName(iBkgnd);
1827  const LauPdfList& pdfList = bkgndPdfs_[iBkgnd];
1828  for (LauPdfList::const_iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) {
1829  // Count the number of input variables that are not DP variables
1830  UInt_t nVars(0);
1831  std::vector<TString> varNames = (*pdf_iter)->varNames();
1832  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
1833  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
1834  ++nVars;
1835  }
1836  }
1837  if ( nVars == 2 ) {
1838  twodimMap.insert( std::make_pair( bkgndClass, std::make_pair( varNames[0], varNames[1] ) ) );
1839  }
1840  }
1841  }
1842  }
1843 
1844  return twodimMap;
1845 }
1846 
1848 {
1849  std::cout << "INFO in LauSimpleFitModel::storePerEvtLlhds : Storing per-event likelihood values..." << std::endl;
1850 
1851  // if we've not been using the DP model then we need to cache all
1852  // the info here so that we can get the efficiency from it
1853 
1854  LauFitDataTree* inputFitData = this->fitData();
1855 
1856  if (!this->useDP() && this->storeDPEff()) {
1858  sigDPModel_->fillDataTree(*inputFitData);
1859  }
1860 
1861  UInt_t evtsPerExpt(this->eventsPerExpt());
1862 
1863  for (UInt_t iEvt = 0; iEvt < evtsPerExpt; ++iEvt) {
1864 
1865  this->setSPlotNtupleIntegerBranchValue("iExpt",this->iExpt());
1866  this->setSPlotNtupleIntegerBranchValue("iEvtWithinExpt",iEvt);
1867 
1868  // the DP information
1869  this->getEvtDPLikelihood(iEvt);
1870  if (this->storeDPEff()) {
1871  if (!this->useDP()) {
1873  }
1874  this->setSPlotNtupleDoubleBranchValue("efficiency",sigDPModel_->getEvtEff());
1875  if ( sigDPModel_->usingScfModel() ) {
1877  }
1878  }
1879  if (this->useDP()) {
1881  if (useSCF_) {
1883  this->setSPlotNtupleDoubleBranchValue("sigTMDPLike",sigDPLike_);
1884  this->setSPlotNtupleDoubleBranchValue("sigSCFDPLike",scfDPLike_);
1885  } else {
1886  this->setSPlotNtupleDoubleBranchValue("sigDPLike",sigDPLike_);
1887  }
1888  if (usingBkgnd_) {
1889  const UInt_t nBkgnds = this->nBkgndClasses();
1890  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
1891  TString name = this->bkgndClassName(iBkgnd);
1892  name += "DPLike";
1893  this->setSPlotNtupleDoubleBranchValue(name,bkgndDPLike_[iBkgnd]);
1894  }
1895  }
1896  } else {
1897  sigTotalLike_ = 1.0;
1898  if (useSCF_) {
1899  scfTotalLike_ = 1.0;
1900  }
1901  if (usingBkgnd_) {
1902  const UInt_t nBkgnds = this->nBkgndClasses();
1903  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
1904  bkgndTotalLike_[iBkgnd] = 1.0;
1905  }
1906  }
1907  }
1908 
1909  // the signal PDF values
1910  if ( useSCF_ ) {
1911  sigTotalLike_ *= this->setSPlotNtupleBranchValues(&signalPdfs_, "sigTM", iEvt);
1912  scfTotalLike_ *= this->setSPlotNtupleBranchValues(&scfPdfs_, "sigSCF", iEvt);
1913  } else {
1914  sigTotalLike_ *= this->setSPlotNtupleBranchValues(&signalPdfs_, "sig", iEvt);
1915  }
1916 
1917  // the background PDF values
1918  if (usingBkgnd_) {
1919  const UInt_t nBkgnds = this->nBkgndClasses();
1920  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
1921  const TString& bkgndClass = this->bkgndClassName(iBkgnd);
1922  LauPdfList& pdfs = bkgndPdfs_[iBkgnd];
1923  bkgndTotalLike_[iBkgnd] *= this->setSPlotNtupleBranchValues(&(pdfs), bkgndClass, iEvt);
1924  }
1925  }
1926 
1927  // the total likelihoods
1928  if (useSCF_) {
1929  Double_t scfFrac(0.0);
1930  if ( useSCFHist_ ) {
1931  scfFrac = recoSCFFracs_[iEvt];
1932  } else {
1933  scfFrac = scfFrac_.value();
1934  }
1935  this->setSPlotNtupleDoubleBranchValue("sigSCFFrac",scfFrac);
1936  sigTotalLike_ *= ( 1.0 - scfFrac );
1937  if ( scfMap_ == 0 ) {
1938  scfTotalLike_ *= scfFrac;
1939  }
1940  this->setSPlotNtupleDoubleBranchValue("sigTMTotalLike",sigTotalLike_);
1941  this->setSPlotNtupleDoubleBranchValue("sigSCFTotalLike",scfTotalLike_);
1942  } else {
1943  this->setSPlotNtupleDoubleBranchValue("sigTotalLike",sigTotalLike_);
1944  }
1945  if (usingBkgnd_) {
1946  const UInt_t nBkgnds = this->nBkgndClasses();
1947  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
1948  TString name = this->bkgndClassName(iBkgnd);
1949  name += "TotalLike";
1951  }
1952  }
1953  // fill the tree
1954  this->fillSPlotNtupleBranches();
1955  }
1956  std::cout << "INFO in LauSimpleFitModel::storePerEvtLlhds : Finished storing per-event likelihood values." << std::endl;
1957 }
1958 
1959 void LauSimpleFitModel::embedSignal(const TString& fileName, const TString& treeName,
1960  Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment,
1961  Bool_t useReweighting)
1962 {
1963  if (signalTree_) {
1964  std::cerr << "ERROR in LauSimpleFitModel::embedSignal : Already embedding signal from a file." << std::endl;
1965  return;
1966  }
1967 
1968  if (!reuseEventsWithinEnsemble && reuseEventsWithinExperiment) {
1969  std::cerr << "WARNING in LauSimpleFitModel::embedSignal : Conflicting options provided, will not reuse events at all." << std::endl;
1970  reuseEventsWithinExperiment = kFALSE;
1971  }
1972 
1973  signalTree_ = new LauEmbeddedData(fileName,treeName,reuseEventsWithinExperiment);
1974 
1975  Bool_t dataOK = signalTree_->findBranches();
1976  if (!dataOK) {
1977  delete signalTree_; signalTree_ = 0;
1978  std::cerr << "ERROR in LauSimpleFitModel::embedSignal : Problem creating data tree for embedding." << std::endl;
1979  return;
1980  }
1981 
1982  reuseSignal_ = reuseEventsWithinEnsemble;
1983  useReweighting_ = useReweighting;
1984 
1985  if (this->enableEmbedding() == kFALSE) {
1986  this->enableEmbedding(kTRUE);
1987  }
1988 }
1989 
1990 void LauSimpleFitModel::embedBkgnd(const TString& bkgndClass, const TString& fileName, const TString& treeName,
1991  Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment)
1992 {
1993  if ( ! this->validBkgndClass( bkgndClass ) ) {
1994  std::cerr << "ERROR in LauSimpleFitModel::embedBkgnd : Invalid background class \"" << bkgndClass << "\"." << std::endl;
1995  std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
1996  return;
1997  }
1998 
1999  UInt_t bkgndID = this->bkgndClassID( bkgndClass );
2000 
2001  if (bkgndTree_[bkgndID]) {
2002  std::cerr << "ERROR in LauSimpleFitModel::embedBkgnd : Already embedding background from a file." << std::endl;
2003  return;
2004  }
2005 
2006  if (!reuseEventsWithinEnsemble && reuseEventsWithinExperiment) {
2007  std::cerr << "WARNING in LauSimpleFitModel::embedBkgnd : Conflicting options provided, will not reuse events at all." << std::endl;
2008  reuseEventsWithinExperiment = kFALSE;
2009  }
2010 
2011  bkgndTree_[bkgndID] = new LauEmbeddedData(fileName,treeName,reuseEventsWithinExperiment);
2012 
2013  Bool_t dataOK = bkgndTree_[bkgndID]->findBranches();
2014  if (!dataOK) {
2015  delete bkgndTree_[bkgndID]; bkgndTree_[bkgndID] = 0;
2016  std::cerr << "ERROR in LauSimpleFitModel::embedBkgnd : Problem creating data tree for embedding." << std::endl;
2017  return;
2018  }
2019 
2020  reuseBkgnd_[bkgndID] = reuseEventsWithinEnsemble;
2021 
2022  if (this->enableEmbedding() == kFALSE) {
2023  this->enableEmbedding(kTRUE);
2024  }
2025 }
2026 
2027 void LauSimpleFitModel::weightEvents( const TString& dataFileName, const TString& dataTreeName )
2028 {
2029  // Routine to provide weights for events that are uniformly distributed
2030  // in the DP (or square DP) so as to reproduce the given DP model
2031 
2032  if ( kinematics_->squareDP() ) {
2033  std::cout << "INFO in LauSimpleFitModel::weightEvents : will create weights assuming events were generated flat in the square DP" << std::endl;
2034  } else {
2035  std::cout << "INFO in LauSimpleFitModel::weightEvents : will create weights assuming events were generated flat in phase space" << std::endl;
2036  }
2037 
2038  // This reads in the given dataFile and creates an input
2039  // fit data tree that stores them for all events and experiments.
2040  Bool_t dataOK = this->cacheFitData(dataFileName,dataTreeName);
2041  if (!dataOK) {
2042  std::cerr << "ERROR in LauSimpleFitModel::weightEvents : Problem caching the data." << std::endl;
2043  return;
2044  }
2045 
2046  LauFitDataTree* inputFitData = this->fitData();
2047 
2048  if ( ! inputFitData->haveBranch( "m13Sq_MC" ) || ! inputFitData->haveBranch( "m23Sq_MC" ) ) {
2049  std::cerr << "WARNING in LauSimpleFitModel::weightEvents : Cannot find MC truth DP coordinate branches in supplied data, aborting." << std::endl;
2050  return;
2051  }
2052 
2053  // Create the ntuple to hold the DP weights
2054  TString weightsFileName( dataFileName );
2055  Ssiz_t index = weightsFileName.Last('.');
2056  weightsFileName.Insert( index, "_DPweights" );
2057  LauGenNtuple * weightsTuple = new LauGenNtuple( weightsFileName, dataTreeName );
2058  weightsTuple->addIntegerBranch("iExpt");
2059  weightsTuple->addIntegerBranch("iEvtWithinExpt");
2060  weightsTuple->addDoubleBranch("dpModelWeight");
2061 
2062  UInt_t iExpmt = this->iExpt();
2063  UInt_t nExpmt = this->nExpt();
2064  UInt_t firstExpmt = this->firstExpt();
2065  for (iExpmt = firstExpmt; iExpmt < (firstExpmt+nExpmt); ++iExpmt) {
2066 
2067  inputFitData->readExperimentData(iExpmt);
2068  UInt_t nEvents = inputFitData->nEvents();
2069 
2070  if (nEvents < 1) {
2071  std::cerr << "WARNING in LauSimpleFitModel::weightEvents : Zero events in experiment " << iExpmt << ", skipping..." << std::endl;
2072  continue;
2073  }
2074 
2075  weightsTuple->setIntegerBranchValue( "iExpt", iExpmt );
2076 
2077  // Calculate and store the weights for the events in this experiment
2078  for ( UInt_t iEvent(0); iEvent < nEvents; ++iEvent ) {
2079 
2080  weightsTuple->setIntegerBranchValue( "iEvtWithinExpt", iEvent );
2081 
2082  const LauFitData& evtData = inputFitData->getData( iEvent );
2083 
2084  Double_t m13Sq_MC = evtData.find("m13Sq_MC")->second;
2085  Double_t m23Sq_MC = evtData.find("m23Sq_MC")->second;
2086 
2087  Double_t dpModelWeight(0.0);
2088 
2089  if ( kinematics_->withinDPLimits( m13Sq_MC, m23Sq_MC ) ) {
2090 
2091  kinematics_->updateKinematics( m13Sq_MC, m23Sq_MC );
2092  dpModelWeight = sigDPModel_->getEventWeight();
2093 
2094  if ( kinematics_->squareDP() ) {
2095  dpModelWeight *= kinematics_->calcSqDPJacobian();
2096  }
2097 
2098  if (LauIsobarDynamics::GenOK != sigDPModel_->checkToyMC(kTRUE,kFALSE)) {
2099  std::cerr << "WARNING in LauSimpleFitModel::weightEvents : Problem in calculating the weight, aborting." << std::endl;
2100  delete weightsTuple;
2101  return;
2102  }
2103  }
2104 
2105  weightsTuple->setDoubleBranchValue( "dpModelWeight", dpModelWeight );
2106  weightsTuple->fillBranches();
2107  }
2108 
2109  }
2110 
2111  weightsTuple->buildIndex( "iExpt", "iEvtWithinExpt" );
2112  weightsTuple->addFriendTree(dataFileName, dataTreeName);
2113  weightsTuple->writeOutGenResults();
2114 
2115  delete weightsTuple;
2116 }
2117 
void setSignalPdf(LauAbsPdf *pdf)
Set the signal PDF for a given variable.
Double_t range() const
The range allowed for the parameter.
void setBkgndDPModel(const TString &bkgndClass, LauAbsBkgndDPModel *bkgndModel)
Set the background DP models.
Bool_t reuseSignal_
Boolean to reuse signal events.
void modifyDataTree()
Recache the amplitude values for those that have changed.
void setDPBranchValues()
Store all of the DP information.
virtual void addSPlotNtupleIntegerBranch(const TString &name)
Add a branch to the sPlot tree for storing an integer.
The abstract interface for a background Dalitz plot model.
Bool_t squareDP() const
Are the square Dalitz plot co-ordinates being calculated?
UInt_t firstExpt() const
Obtain the number of the first experiment.
Double_t getc23() const
Get the cosine of the helicity angle theta23.
LauParameter meanEff_
The mean efficiency.
TRandom * randomFun()
Access the singleton random number generator with a particular seed.
Definition: LauRandom.cc:20
Double_t sigExtraLike_
Signal likelihood from extra PDFs.
std::vector< LauParameter > LauParameterList
List of parameters.
Bool_t fixed() const
Check whether the parameter is fixed or floated.
UInt_t nEvents() const
Retrieve the number of events.
UInt_t eventsPerExpt() const
Obtain the total number of events in the current experiment.
Bool_t writeLatexTable() const
Determine whether writing out of the latex table is enabled.
Double_t sigDPLike_
Signal likelihood value.
void randomiseInitFitPars()
Randomise the initial fit parameters.
LauEffModel * scfFracHist_
The histogram giving the DP-dependence of the SCF fraction.
void updateCoeffs(const std::vector< LauComplex > &coeffs)
Update the complex coefficients for the resonances.
void updateSqDPKinematics(Double_t mPrime, Double_t thetaPrime)
Update all kinematic quantities based on the square DP co-ordinates m&#39; and theta&#39;.
virtual Double_t getTotEvtLikelihood(UInt_t iEvt)
Get the total likelihood for each event.
void setEffHisto(const TH2 *effHisto, Bool_t useInterpolation=kTRUE, Bool_t fluctuateBins=kFALSE, Double_t avEff=-1.0, Double_t absError=-1.0, Bool_t useUpperHalfOnly=kFALSE, Bool_t squareDP=kFALSE)
Set the efficiency variation across the phase space using a predetermined 2D histogram.
Definition: LauEffModel.cc:55
void cacheInfo(LauPdfList &pdfList, const LauFitDataTree &theData)
Have all PDFs in the list cache the data.
virtual void setGenNtupleIntegerBranchValue(const TString &name, Int_t value)
Set the value of an integer branch in the gen tree.
ClassImp(LauAbsCoeffSet)
Bool_t getReweightedEvent(LauIsobarDynamics *dynamics)
Retrieve an event from the data sample, applying an accept/reject based on the given DP model...
std::vector< Double_t > recoJacobians_
The cached values of the sqDP jacobians for each event.
Bool_t hasResonance(const TString &resName) const
Check whether this model includes a named resonance.
void generateExtraPdfValues(LauPdfList *extraPdfs, LauEmbeddedData *embeddedData)
Generate from the extra PDFs.
Double_t getm23() const
Get the m23 invariant mass.
virtual void initialise()
Initialise the fit.
const TString & name() const
The parameter name.
Double_t getmPrime() const
Get m&#39; value.
Class that defines the particular 3-body decay under study.
Definition: LauDaughters.hh:33
UInt_t addFitParameters(LauPdfList &pdfList)
Add parameters of the PDFs in the list to the list of all fit parameters.
void splitSignalComponent(const TH2 *dpHisto, const Bool_t upperHalf=kFALSE, const Bool_t fluctuateBins=kFALSE, LauScfMap *scfMap=0)
Split the signal component into well-reconstructed and mis-reconstructed parts.
std::map< TString, std::pair< Int_t, Double_t > > LauGenInfo
Define a map to be used to store a category name and numbers.
Bool_t useDP() const
Is the Dalitz plot term in the likelihood.
virtual void setNBkgndEvents(LauParameter *nBkgndEvents)
Set the background event yield(s)
Double_t getc13() const
Get the cosine of the helicity angle theta13.
File containing declaration of LauAbsCoeffSet class.
void clearUsedList()
Clear the list of used events.
virtual void weightEvents(const TString &dataFileName, const TString &dataTreeName)
Weight events based on the DP model.
void calcExtraInfo(const Bool_t init=kFALSE)
Calculate the fit fractions, mean efficiency and total DP rate.
TH2 * trueHist(Int_t trueBin)
Retrieve the migration histogram for trueBin.
Definition: LauScfMap.cc:176
virtual void setAmpCoeffSet(LauAbsCoeffSet *coeffSet)
Set the DP amplitude coefficients.
std::multimap< TString, std::pair< TString, TString > > TwoDMap
Type to associate the name of the species that have 2D PDFs with the names of the two variables invol...
Definition: LauSPlot.hh:68
UInt_t nExtraPdfPar_
Number of extra PDF parameters.
Int_t buildIndex(const TString &majorName, const TString &minorName="0")
Create an index table using leaves of the tree.
Double_t scfDPLike_
SCF likelihood value.
void listBinCentres(std::vector< Double_t > &xCoords, std::vector< Double_t > &yCoords) const
Create lists of the co-ordinates of the centres of true bins.
Definition: LauScfMap.cc:114
File containing declaration of LauDaughters class.
const std::vector< Int_t > * trueBins(Int_t recoBin) const
Find which true bins contribute to the given reco bin.
Definition: LauScfMap.cc:155
virtual void setGenNtupleDoubleBranchValue(const TString &name, Double_t value)
Set the value of a double branch in the gen tree.
virtual void finaliseFitResults(const TString &tablePrefixName)
Get the fit results and store them.
virtual void propagateParUpdates()
Calculate things that depend on the fit parameters after they have been updated by Minuit...
Bool_t useSCFHist_
Is the SCF fraction DP-dependent.
void addFriendTree(const TString &rootFileName, const TString &rootTreeName)
Add a friend tree.
File containing declaration of LauScfMap class.
Double_t getEvtthPrime() const
Retrieve the square Dalitz plot coordinate, theta&#39;, for the current event.
LauBkgndReuseEventsList reuseBkgnd_
Vector of booleans to reuse background events.
std::vector< std::vector< LauParameter > > LauParArray
Type to define an array of parameters.
File containing declaration of LauSimpleFitModel class.
LauAbsEffModel * getEffModel()
Retrieve the model for the efficiency across the Dalitz plot.
Int_t fitStatus() const
Access the fit status information.
virtual void printParValues() const =0
Print the current values of the parameters.
File containing declaration of LauPrint class.
Pure abstract base class for defining the efficiency description across the signal Dalitz plot...
void setBkgndPdf(const TString &bkgndClass, LauAbsPdf *pdf)
Set the background PDF.
std::vector< LauAbsPdf * > LauPdfList
List of Pdfs.
virtual void updateCoeffs()
Update the coefficients.
std::map< TString, Double_t > LauFitData
Type for holding event data.
void fillBranches()
Fill branches in the ntuple.
virtual void getEvtDPLikelihood(UInt_t iEvt)
Calculate the signal and background likelihoods for the DP for a given event.
LauEmbeddedData * signalTree_
The signal event tree.
void setSCFPdf(LauAbsPdf *pdf)
Set the SCF PDF for a given variable.
virtual void setNSigEvents(LauParameter *nSigEvents)
Set the signal event yield.
Double_t getm23Sq() const
Get the m23 invariant mass square.
Bool_t useReweighting_
Boolean to use reweighting.
virtual void printFitFractions(std::ostream &output)
Print the fit fractions, total DP rate and mean efficiency.
const LauParArray & getFitFractions() const
Retrieve the fit fractions for the amplitude components.
void recalculateNormalisation()
recalculate Normalization
Bool_t squareDP() const
Determine to use or not the square Dalitz plot.
Definition: LauDaughters.hh:78
void addSPlotNtupleBranches(const LauPdfList *extraPdfs, const TString &prefix)
Add sPlot branches for the extra PDFs.
virtual LauSPlot::TwoDMap twodimPDFs() const
Returns the species and variables for all 2D PDFs in the fit.
LauIsobarDynamics * sigDPModel_
The signal Dalitz plot model.
Double_t getEvtScfFraction() const
Retrieve the fraction of events that are poorly reconstructed (the self cross feed fraction) for the ...
Class to define various output print commands.
Definition: LauPrint.hh:29
Double_t getValue(const TString &name) const
Get the value of a specified branch.
LauDaughters * getDaughters()
Retrieve the daughters.
void calcLikelihoodInfo(const UInt_t iEvt)
Calculate the likelihood (and all associated information) for the given event number.
std::vector< LauComplex > coeffs_
The complex coefficients.
virtual void fillGenNtupleBranches()
Fill the gen tuple branches.
LauParameter getMeanEff() const
Retrieve the mean efficiency across the Dalitz plot.
std::vector< Double_t > bkgndExtraLike_
Background likelihood value(s) from extra PDFs.
void updateFitNtuple()
Update the fit ntuple.
std::set< LauParameter * > LauParameterPSet
Set of parameter pointers.
File containing declaration of LauKinematics class.
virtual void setSPlotNtupleDoubleBranchValue(const TString &name, Double_t value)
Set the value of a double branch in the sPlot tree.
File containing declaration of LauEmbeddedData class.
LauBkgndYieldList bkgndEvents_
Background yield(s)
void printFitParameters(const LauPdfList &pdfList, std::ostream &fout) const
Print the fit parameters for all PDFs in the list.
LauPdfList signalPdfs_
The signal PDFs.
Bool_t haveBranch(const TString &name) const
Boolean to determine whether branch exists.
void setExtraPdfParameters()
Set the fit parameters for the extra PDFs.
Bool_t generateBkgndEvent(UInt_t bkgndID)
Generate background event.
Double_t scfTotalLike_
Total SCF likelihood.
Double_t getEvtmPrime() const
Retrieve the square Dalitz plot coordinate, m&#39;, for the current event.
void readExperimentData(UInt_t iExpt)
Read events only for the given experiment.
void updateSigEvents()
Update the signal events after Minuit sets background parameters.
File containing declaration of LauIsobarDynamics class.
virtual UInt_t index() const
Retrieve the index number of the coefficient set.
virtual LauSPlot::NumbMap fixdSpeciesNames() const
Returns the names and yields of species that are fixed in the fit.
void updateKinematics(Double_t m13Sq, Double_t m23Sq)
Update all kinematic quantities based on the DP co-ordinates m13Sq and m23Sq.
void fillDataTree(const LauFitDataTree &fitDataTree)
Fill the internal data structure that caches the resonance dynamics.
LauParameter dpRate_
The average DP rate.
const LauParameterList & extraPars() const
Access the extra variables.
std::vector< Double_t > bkgndDPLike_
Background likelihood value(s)
std::vector< LauParameter * > LauParameterPList
List of parameter pointers.
UInt_t nNormPar_
Number of normalisation parameters (i.e. yields)
Double_t getEvtm23Sq() const
Retrieve the invariant mass squared of the second and third daughters in the current event...
Double_t calcEfficiency(const LauKinematics *kinematics) const
Determine the efficiency for a given point in the Dalitz plot.
Definition: LauEffModel.cc:261
void setupGenNtupleBranches()
Setup the required ntuple branches.
void addIntegerBranch(const TString &name)
Add integer branch to tree.
Definition: LauGenNtuple.cc:73
void appendBinCentres(LauFitDataTree *inputData)
Append fake data points to the inputData for each bin in the SCF smearing matrix. ...
virtual void cacheInputFitVars()
Read in the input fit data variables, e.g. m13Sq and m23Sq.
Abstract interface to the fitting and toy MC model.
std::vector< LauParameter * > & getFloatingParameters()
Retrieve the floating parameters of the resonance models.
void setSignalDPParameters()
Set the fit parameters for the DP model.
UInt_t iExpt() const
Obtain the number of the current experiment.
virtual void addGenNtupleIntegerBranch(const TString &name)
Add a branch to the gen tree for storing an integer.
const LauFitNtuple * fitNtuple() const
Access the fit ntuple.
std::vector< LauAbsCoeffSet * > coeffPars_
Magnitudes and Phases.
Int_t binNumber(Double_t xCoord, Double_t yCoord) const
Find the global bin number for the given co-ordinates.
Definition: LauScfMap.cc:144
virtual LauSPlot::NumbMap freeSpeciesNames() const
Returns the names and yields of species that are free in the fit.
UInt_t nExpt() const
Obtain the number of experiments.
virtual LauSPlot::NameSet variableNames() const
Returns the names of all variables in the fit.
Bool_t doPoissonSmearing() const
Determine whether Poisson smearing is enabled for the toy MC generation.
virtual void getEvtExtraLikelihoods(UInt_t iEvt)
Determine the signal and background likelihood for the extra variables for a given event...
LauBkgndPdfsList bkgndPdfs_
The background PDFs.
Double_t getc12() const
Get the cosine of the helicity angle theta12.
Double_t getm13() const
Get the m13 invariant mass.
std::map< TString, Double_t > NumbMap
Type to associate a category name with a double precision number, e.g. a yield or PDF value for a giv...
Definition: LauSPlot.hh:62
Bool_t doEMLFit() const
Determine whether an extended maximum likelihood fit it being performed.
LauParameter getDPRate() const
Retrieve the overall Dalitz plot rate.
Bool_t useSCF_
Is the signal split into TM and SCF.
File containing declaration of LauComplex class.
LauFitData getValues(const std::vector< TString > &names) const
Get values of specified branches.
virtual void setupBkgndVectors()
Define the length of the background vectors.
LauGenInfo eventsToGenerate()
Determine the number of events to generate for each hypothesis.
const LauParArray & getFitFractionsEfficiencyUncorrected() const
Retrieve the fit fractions for the amplitude components.
Class to store the results from the toy MC generation into an ntuple.
Definition: LauGenNtuple.hh:32
Class for defining the abstract interface for complex coefficient classes.
Double_t nll() const
Access the current NLL value.
std::vector< Double_t > fakeJacobians_
The cached values of the sqDP jacobians for each true bin.
void storeCorrMatrix(UInt_t iExpt, Double_t NLL, Int_t fitStatus, const TMatrixD &covMatrix)
Store the correlation matrix and other fit information.
Definition: LauFitNtuple.cc:64
void embedSignal(const TString &fileName, const TString &treeName, Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment=kFALSE, Bool_t useReweighting=kFALSE)
Embed full simulation events for the signal, rather than generating toy from the PDFs.
std::vector< Double_t > bkgndTotalLike_
Total background likelihood(s)
void clearFitParVectors()
Clear the vectors containing fit parameters.
Class for defining the fit parameter objects.
Definition: LauParameter.hh:34
LauBkgndDPModelList bkgndDPModels_
The background Dalitz Plot model.
Class to store the data for embedding in toy experiments.
UInt_t getnTotAmp() const
Retrieve the total number of amplitude components.
virtual void addSPlotNtupleDoubleBranch(const TString &name)
Add a branch to the sPlot tree for storing a double.
LauParArray fitFracEffUnCorr_
Fit fractions (uncorrected for the efficiency)
void clearExtraVarVectors()
Clear the vectors containing extra ntuple variables.
void getEmbeddedEvent(LauKinematics *kinematics)
Retrieve an event from the data sample.
Bool_t withinDPLimits(Double_t m13Sq, Double_t m23Sq) const
Check whether a given (m13Sq,m23Sq) point is within the kinematic limits of the Dalitz plot...
virtual void recalculateNormalisation()
Recalculate Normalization the signal DP models.
std::vector< LauParameter > getExtraParameters()
Retrieve any extra parameters/quantities (e.g. K-matrix total fit fractions)
LauScfMap * scfMap_
The smearing matrix for the SCF DP PDF.
UInt_t nSigComp_
Number of signal components.
void appendFakePoints(const std::vector< Double_t > &xCoords, const std::vector< Double_t > &yCoords)
Add fake events to the data.
File containing declaration of LauAbsPdf class.
Bool_t storeDPEff() const
Determine whether the efficiency information should be stored in the sPlot ntuple.
virtual Double_t getEvtSCFDPLikelihood(UInt_t iEvt)
Calculate the SCF likelihood for the DP for a given event.
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.
File containing LauRandom namespace.
Bool_t generate()
Generate a toy MC signal event.
Bool_t useRandomInitFitPars() const
Determine whether the initial values of the fit parameters, in particular the isobar coefficient para...
File containing declaration of LauEffModel class.
const LauFitDataTree * fitData() const
Access the data store.
void setIntegerBranchValue(const TString &name, Int_t value)
Set value of an integer branch.
Definition: LauGenNtuple.cc:91
Double_t setSPlotNtupleBranchValues(LauPdfList *extraPdfs, const TString &prefix, UInt_t iEvt)
Set the branches for the sPlot ntuple with extra PDFs.
Bool_t validBkgndClass(const TString &className) const
Check if the given background class is in the list.
virtual void writeOutTable(const TString &outputFile)
Write the fit results in latex table format.
Double_t getm12Sq() const
Get the m12 invariant mass square.
Define a Dalitz plot according to the isobar model.
LauPdfList scfPdfs_
The SCF PDFs.
Class that implements the efficiency description across the signal Dalitz plot.
Definition: LauEffModel.hh:37
Bool_t cacheFitData(const TString &dataFileName, const TString &dataTreeName)
Store variables from the input file into the internal data storage.
const LauParameterPSet & resPars() const
Access the fit variables which affect the DP normalisation.
UInt_t nUsedEvents() const
Retrieve the number of events that have already been sampled.
Bool_t usingBkgnd_
Background boolean.
LauParameter * signalEvents_
Signal yield.
Double_t getm13Sq() const
Get the m13 invariant mass square.
void storeParsAndErrors(const std::vector< LauParameter * > &fitVars, const std::vector< LauParameter > &extraVars)
Store parameters and their errors.
std::vector< Double_t > recoSCFFracs_
The cached values of the SCF fraction for each event.
virtual void setSPlotNtupleIntegerBranchValue(const TString &name, Int_t value)
Set the value of an integer branch in the sPlot tree.
void embedBkgnd(const TString &bkgndClass, const TString &fileName, const TString &treeName, Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment=kFALSE)
Embed full simulation events for the given background class, rather than generating toy from the PDFs...
Double_t initValue() const
The initial value of the parameter.
Double_t scfExtraLike_
SCF likelihood from extra PDFs.
Class for defining signal dynamics using the isobar model.
Double_t getm12() const
Get the m12 invariant mass.
UInt_t nSigDPPar_
Number of signal DP parameters.
Double_t getEventWeight()
Calculate the acceptance rate, for events with the current kinematics, when generating events accordi...
virtual void fillSPlotNtupleBranches()
Fill the sPlot tuple.
File containing LauConstants namespace.
std::vector< Double_t > fakeSCFFracs_
The cached values of the SCF fraction for each bin centre.
LauParArray fitFrac_
Fit fractions.
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
ToyMCStatus checkToyMC(Bool_t printErrorMessages=kTRUE, Bool_t printInfoMessages=kFALSE)
Check the status of the toy MC generation.
Class to store the results from the fit into an ntuple.
Definition: LauFitNtuple.hh:42
Double_t sigTotalLike_
Total signal likelihood.
Bool_t haveBranch(const TString &name) const
Check if the named branch is stored.
Class for representing the 4D smearing matrix for mis-reconstructed signal (self cross feed) ...
Definition: LauScfMap.hh:29
LauParameter scfFrac_
The (global) SCF fraction parameter.
virtual ~LauSimpleFitModel()
Destructor.
Double_t getThetaPrime() const
Get theta&#39; value.
void updatePull()
Call to update the bias and pull values.
Double_t getEvtEff() const
Retrieve the efficiency for the current event.
UInt_t nBkgndClasses() const
Returns the number of background classes.
const TString & bkgndClassName(UInt_t classID) const
Get the name of a background class from the number.
Double_t getEvtLikelihood() const
Retrieve the likelihood for the current event.
Bool_t enableEmbedding() const
Determine whether embedding of events is enabled in the generation.
virtual TString name() const
Retrieve the name of the coefficient set.
Class for defining the abstract interface for PDF classes.
Definition: LauAbsPdf.hh:41
const TMatrixD & covarianceMatrix() const
Access the fit covariance matrix.
const LauParameterPList & fitPars() const
Access the fit variables.
Bool_t usingScfModel()
Check whether a self cross feed fraction model is being used.
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.
Double_t calcSqDPJacobian()
Calculate the Jacobian for the transformation m23^2, m13^2 -&gt; m&#39;, theta&#39; (square DP) ...
virtual void setupSPlotNtupleBranches()
Add branches to store experiment number and the event number within the experiment.
void updateFitParameters(LauPdfList &pdfList)
Update the fit parameters for all PDFs in the list.
virtual void checkInitFitParams()
Check the initial fit parameters.
void setExtraNtupleVars()
Set-up other parameters that are derived from the fit results, e.g. fit fractions.
void addDoubleBranch(const TString &name)
Add double branch to tree.
Definition: LauGenNtuple.cc:82
virtual void storePerEvtLlhds()
Store the per event likelihood values.
virtual Bool_t genExpt()
Toy MC generation and fitting overloaded functions.
Bool_t findBranches()
Find and read the branches in data tree.
File containing declaration of LauAbsBkgndDPModel class.
void setFitNEvents()
Set the initial yields.
Double_t prob(Int_t recoBin, Int_t trueBin) const
Probability of a true event in the given true bin migrating to the reco bin.
Definition: LauScfMap.cc:165
void initialise(const std::vector< LauComplex > &coeffs)
Initialise the Dalitz plot dynamics.
Double_t getEvtm13Sq() const
Retrieve the invariant mass squared of the first and third daughters in the current event...
Double_t genValue() const
The value generated for the parameter.
virtual Double_t getEventSum() const
Get the total number of events.
LauKinematics * kinematics_
The Dalitz plot kinematics object.
Class to store the input fit variables.
std::set< TString > NameSet
Type to store names, e.g. of the discriminating/control variables.
Definition: LauSPlot.hh:59
File containing declaration of LauFitNtuple class.
virtual void initialiseDPModels()
Initialise the signal DP model.
File containing declaration of LauGenNtuple class.
void setDoubleBranchValue(const TString &name, Double_t value)
Set value of a double branch.
Definition: LauGenNtuple.cc:96
LauBkgndEmbDataList bkgndTree_
The background event tree.
virtual Bool_t generate()=0
Generate a toy MC event from the model.
virtual Bool_t passVeto(const LauKinematics *kinematics) const =0
Determine whether the given DP position is outside the vetoes.
Bool_t generateSignalEvent()
Generate signal event.