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