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