laura is hosted by Hepforge, IPPP Durham
Laura++  v3r0
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauCPFitModel.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2004 - 2014.
3 // Distributed under the Boost Software License, Version 1.0.
4 // (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5 
6 // Authors:
7 // Thomas Latham
8 // John Back
9 // Paul Harrison
10 
15 #include <iostream>
16 #include <iomanip>
17 #include <fstream>
18 #include <vector>
19 
20 #include "TVirtualFitter.h"
21 #include "TSystem.h"
22 #include "TMinuit.h"
23 #include "TRandom.h"
24 #include "TFile.h"
25 #include "TMath.h"
26 #include "TH2.h"
27 
28 #include "LauAbsBkgndDPModel.hh"
29 #include "LauAbsCoeffSet.hh"
30 #include "LauIsobarDynamics.hh"
31 #include "LauAbsPdf.hh"
32 #include "LauAsymmCalc.hh"
33 #include "LauComplex.hh"
34 #include "LauConstants.hh"
35 #include "LauCPFitModel.hh"
36 #include "LauDaughters.hh"
37 #include "LauEffModel.hh"
38 #include "LauEmbeddedData.hh"
39 #include "LauFitNtuple.hh"
40 #include "LauKinematics.hh"
41 #include "LauPrint.hh"
42 #include "LauRandom.hh"
43 #include "LauScfMap.hh"
44 
46 
47 
48 LauCPFitModel::LauCPFitModel(LauIsobarDynamics* negModel, LauIsobarDynamics* posModel, Bool_t tagged, const TString& tagVarName) : LauAbsFitModel(),
49  negSigModel_(negModel), posSigModel_(posModel),
50  negKinematics_(negModel ? negModel->getKinematics() : 0),
51  posKinematics_(posModel ? posModel->getKinematics() : 0),
52  usingBkgnd_(kFALSE),
53  nSigComp_(0),
54  nSigDPPar_(0),
55  nExtraPdfPar_(0),
56  nNormPar_(0),
57  negMeanEff_("negMeanEff",0.0,0.0,1.0), posMeanEff_("posMeanEff",0.0,0.0,1.0),
58  negDPRate_("negDPRate",0.0,0.0,100.0), posDPRate_("posDPRate",0.0,0.0,100.0),
59  signalEvents_(0),
60  signalAsym_(0),
61  forceAsym_(kFALSE),
62  tagged_(tagged),
63  tagVarName_(tagVarName),
64  curEvtCharge_(0),
65  useSCF_(kFALSE),
66  useSCFHist_(kFALSE),
67  scfFrac_("scfFrac",0.0,0.0,1.0),
68  scfFracHist_(0),
69  scfMap_(0),
70  compareFitData_(kFALSE),
71  negParent_("B-"), posParent_("B+"),
72  negSignalTree_(0), posSignalTree_(0),
73  reuseSignal_(kFALSE),
74  useNegReweighting_(kFALSE), usePosReweighting_(kFALSE),
75  sigDPLike_(0.0),
76  scfDPLike_(0.0),
77  sigExtraLike_(0.0),
78  scfExtraLike_(0.0),
79  sigTotalLike_(0.0),
80  scfTotalLike_(0.0)
81 {
82  LauDaughters* negDaug = negSigModel_->getDaughters();
83  if (negDaug != 0) {negParent_ = negDaug->getNameParent();}
84  LauDaughters* posDaug = posSigModel_->getDaughters();
85  if (posDaug != 0) {posParent_ = posDaug->getNameParent();}
86 }
87 
89 {
90  delete negSignalTree_;
91  delete posSignalTree_;
92  for (LauBkgndEmbDataList::iterator iter = negBkgndTree_.begin(); iter != negBkgndTree_.end(); ++iter) {
93  delete (*iter);
94  }
95  for (LauBkgndEmbDataList::iterator iter = posBkgndTree_.begin(); iter != posBkgndTree_.end(); ++iter) {
96  delete (*iter);
97  }
98  delete scfFracHist_;
99 }
100 
102 {
103  UInt_t nBkgnds = this->nBkgndClasses();
104  negBkgndDPModels_.resize( nBkgnds );
105  posBkgndDPModels_.resize( nBkgnds );
106  negBkgndPdfs_.resize( nBkgnds );
107  posBkgndPdfs_.resize( nBkgnds );
108  bkgndEvents_.resize( nBkgnds );
109  bkgndAsym_.resize( nBkgnds );
110  negBkgndTree_.resize( nBkgnds );
111  posBkgndTree_.resize( nBkgnds );
112  reuseBkgnd_.resize( nBkgnds );
113  bkgndDPLike_.resize( nBkgnds );
114  bkgndExtraLike_.resize( nBkgnds );
115  bkgndTotalLike_.resize( nBkgnds );
116 }
117 
119 {
120  if ( nSigEvents == 0 ) {
121  std::cerr << "ERROR in LauCPFitModel::setNSigEvents : The LauParameter pointer is null." << std::endl;
122  gSystem->Exit(EXIT_FAILURE);
123  }
124  if ( signalEvents_ != 0 ) {
125  std::cerr << "ERROR in LauCPFitModel::setNSigEvents : You are trying to overwrite the signal yield." << std::endl;
126  return;
127  }
128  if ( signalAsym_ != 0 ) {
129  std::cerr << "ERROR in LauCPFitModel::setNSigEvents : You are trying to overwrite the signal asymmetry." << std::endl;
130  return;
131  }
132 
133  signalEvents_ = nSigEvents;
134  TString name = signalEvents_->name();
135  if ( ! name.Contains("signalEvents") && !( name.BeginsWith("signal") && name.EndsWith("Events") ) ) {
136  signalEvents_->name("signalEvents");
137  }
138  Double_t value = nSigEvents->value();
139  signalEvents_->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0));
140 
141  signalAsym_ = new LauParameter("signalAsym",0.0,-1.0,1.0,kTRUE);
142 }
143 
144 void LauCPFitModel::setNSigEvents( LauParameter* nSigEvents, LauParameter* sigAsym, Bool_t forceAsym )
145 {
146  if ( nSigEvents == 0 ) {
147  std::cerr << "ERROR in LauCPFitModel::setNSigEvents : The event LauParameter pointer is null." << std::endl;
148  gSystem->Exit(EXIT_FAILURE);
149  }
150  if ( sigAsym == 0 ) {
151  std::cerr << "ERROR in LauCPFitModel::setNSigEvents : The asym LauParameter pointer is null." << std::endl;
152  gSystem->Exit(EXIT_FAILURE);
153  }
154 
155  if ( signalEvents_ != 0 ) {
156  std::cerr << "ERROR in LauCPFitModel::setNSigEvents : You are trying to overwrite the signal yield." << std::endl;
157  return;
158  }
159  if ( signalAsym_ != 0 ) {
160  std::cerr << "ERROR in LauCPFitModel::setNSigEvents : You are trying to overwrite the signal asymmetry." << std::endl;
161  return;
162  }
163 
164  signalEvents_ = nSigEvents;
165  signalEvents_->name("signalEvents");
166  Double_t value = nSigEvents->value();
167  signalEvents_->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0));
168 
169  signalAsym_ = sigAsym;
170  signalAsym_->name("signalAsym");
171  signalAsym_->range(-1.0,1.0);
172 
173  forceAsym_ = forceAsym;
174 }
175 
177 {
178  if ( nBkgndEvents == 0 ) {
179  std::cerr << "ERROR in LauCPFitModel::setNBgkndEvents : The background yield LauParameter pointer is null." << std::endl;
180  gSystem->Exit(EXIT_FAILURE);
181  }
182 
183  if ( ! this->validBkgndClass( nBkgndEvents->name() ) ) {
184  std::cerr << "ERROR in LauCPFitModel::setNBkgndEvents : Invalid background class \"" << nBkgndEvents->name() << "\"." << std::endl;
185  std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
186  gSystem->Exit(EXIT_FAILURE);
187  }
188 
189  UInt_t bkgndID = this->bkgndClassID( nBkgndEvents->name() );
190 
191  if ( bkgndEvents_[bkgndID] != 0 ) {
192  std::cerr << "ERROR in LauCPFitModel::setNBkgndEvents : You are trying to overwrite the background yield." << std::endl;
193  return;
194  }
195 
196  if ( bkgndAsym_[bkgndID] != 0 ) {
197  std::cerr << "ERROR in LauCPFitModel::setNBkgndEvents : You are trying to overwrite the background asymmetry." << std::endl;
198  return;
199  }
200 
201  bkgndEvents_[bkgndID] = nBkgndEvents;
202  bkgndEvents_[bkgndID]->name( nBkgndEvents->name()+"Events" );
203  Double_t value = nBkgndEvents->value();
204  bkgndEvents_[bkgndID]->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0));
205 
206  bkgndAsym_[bkgndID] = new LauParameter(nBkgndEvents->name()+"Asym",0.0,-1.0,1.0,kTRUE);
207 }
208 
210 {
211  if ( nBkgndEvents == 0 ) {
212  std::cerr << "ERROR in LauCPFitModel::setNBkgndEvents : The background yield LauParameter pointer is null." << std::endl;
213  gSystem->Exit(EXIT_FAILURE);
214  }
215 
216  if ( bkgndAsym == 0 ) {
217  std::cerr << "ERROR in LauCPFitModel::setNBkgndEvents : The background asym LauParameter pointer is null." << std::endl;
218  gSystem->Exit(EXIT_FAILURE);
219  }
220 
221  if ( ! this->validBkgndClass( nBkgndEvents->name() ) ) {
222  std::cerr << "ERROR in LauCPFitModel::setNBkgndEvents : Invalid background class \"" << nBkgndEvents->name() << "\"." << std::endl;
223  std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
224  gSystem->Exit(EXIT_FAILURE);
225  }
226 
227  UInt_t bkgndID = this->bkgndClassID( nBkgndEvents->name() );
228 
229  if ( bkgndEvents_[bkgndID] != 0 ) {
230  std::cerr << "ERROR in LauCPFitModel::setNBkgndEvents : You are trying to overwrite the background yield." << std::endl;
231  return;
232  }
233 
234  if ( bkgndAsym_[bkgndID] != 0 ) {
235  std::cerr << "ERROR in LauCPFitModel::setNBkgndEvents : You are trying to overwrite the background asymmetry." << std::endl;
236  return;
237  }
238 
239  bkgndEvents_[bkgndID] = nBkgndEvents;
240  bkgndEvents_[bkgndID]->name( nBkgndEvents->name()+"Events" );
241  Double_t value = nBkgndEvents->value();
242  bkgndEvents_[bkgndID]->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0));
243 
244  bkgndAsym_[bkgndID] = bkgndAsym;
245  bkgndAsym_[bkgndID]->name( nBkgndEvents->name()+"Asym" );
246  bkgndAsym_[bkgndID]->range(-1.0,1.0);
247 }
248 
249 void LauCPFitModel::splitSignalComponent( const TH2* dpHisto, const Bool_t upperHalf, const Bool_t fluctuateBins, LauScfMap* scfMap )
250 {
251  if ( useSCF_ == kTRUE ) {
252  std::cerr << "ERROR in LauCPFitModel::splitSignalComponent : Have already setup SCF." << std::endl;
253  return;
254  }
255 
256  if ( dpHisto == 0 ) {
257  std::cerr << "ERROR in LauCPFitModel::splitSignalComponent : The histogram pointer is null." << std::endl;
258  return;
259  }
260 
261  LauDaughters* daughters = negSigModel_->getDaughters();
262  scfFracHist_ = new LauEffModel( daughters, 0 );
263  scfFracHist_->setEffHisto( dpHisto, kTRUE, fluctuateBins, 0.0, 0.0, upperHalf, daughters->squareDP() );
264 
265  scfMap_ = scfMap;
266 
267  useSCF_ = kTRUE;
268  useSCFHist_ = kTRUE;
269 }
270 
271 void LauCPFitModel::splitSignalComponent( const Double_t scfFrac, const Bool_t fixed )
272 {
273  if ( useSCF_ == kTRUE ) {
274  std::cerr << "ERROR in LauCPFitModel::splitSignalComponent : Have already setup SCF." << std::endl;
275  return;
276  }
277 
278  scfFrac_.range( 0.0, 1.0 );
279  scfFrac_.value( scfFrac ); scfFrac_.initValue( scfFrac ); scfFrac_.genValue( scfFrac );
280  scfFrac_.fixed( fixed );
281 
282  useSCF_ = kTRUE;
283  useSCFHist_ = kFALSE;
284 }
285 
286 void LauCPFitModel::setBkgndDPModels(const TString& bkgndClass, LauAbsBkgndDPModel* negModel, LauAbsBkgndDPModel* posModel)
287 {
288  if ((negModel==0) || (posModel==0)) {
289  std::cerr << "ERROR in LauCPFitModel::setBkgndDPModels : One or both of the model pointers is null." << std::endl;
290  return;
291  }
292 
293  // check that this background name is valid
294  if ( ! this->validBkgndClass( bkgndClass) ) {
295  std::cerr << "ERROR in LauCPFitModel::setBkgndDPModel : Invalid background class \"" << bkgndClass << "\"." << std::endl;
296  std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
297  return;
298  }
299 
300  UInt_t bkgndID = this->bkgndClassID( bkgndClass );
301  negBkgndDPModels_[bkgndID] = negModel;
302  posBkgndDPModels_[bkgndID] = posModel;
303 
304  usingBkgnd_ = kTRUE;
305 }
306 
308 {
309  if ( tagged_ ) {
310  if (negPdf==0 || posPdf==0) {
311  std::cerr << "ERROR in LauCPFitModel::setSignalPdfs : One or both of the PDF pointers is null." << std::endl;
312  return;
313  }
314  } else {
315  // if we're doing an untagged analysis we will only use the negative PDFs
316  if ( negPdf==0 ) {
317  std::cerr << "ERROR in LauCPFitModel::setSignalPdfs : The negative PDF pointer is null." << std::endl;
318  return;
319  }
320  if ( posPdf!=0 ) {
321  std::cerr << "WARNING in LauCPFitModel::setSignalPdfs : Doing an untagged fit so will not use the positive PDF." << std::endl;
322  }
323  }
324  negSignalPdfs_.push_back(negPdf);
325  posSignalPdfs_.push_back(posPdf);
326 }
327 
329 {
330  if ( tagged_ ) {
331  if (negPdf==0 || posPdf==0) {
332  std::cerr << "ERROR in LauCPFitModel::setSCFPdfs : One or both of the PDF pointers is null." << std::endl;
333  return;
334  }
335  } else {
336  // if we're doing an untagged analysis we will only use the negative PDFs
337  if ( negPdf==0 ) {
338  std::cerr << "ERROR in LauCPFitModel::setSCFPdfs : The negative PDF pointer is null." << std::endl;
339  return;
340  }
341  if ( posPdf!=0 ) {
342  std::cerr << "WARNING in LauCPFitModel::setSCFPdfs : Doing an untagged fit so will not use the positive PDF." << std::endl;
343  }
344  }
345  negScfPdfs_.push_back(negPdf);
346  posScfPdfs_.push_back(posPdf);
347 }
348 
349 void LauCPFitModel::setBkgndPdfs(const TString& bkgndClass, LauAbsPdf* negPdf, LauAbsPdf* posPdf)
350 {
351  if ( tagged_ ) {
352  if (negPdf==0 || posPdf==0) {
353  std::cerr << "ERROR in LauCPFitModel::setBkgndPdfs : One or both of the PDF pointers is null." << std::endl;
354  return;
355  }
356  } else {
357  // if we're doing an untagged analysis we will only use the negative PDFs
358  if ( negPdf==0 ) {
359  std::cerr << "ERROR in LauCPFitModel::setBkgndPdfs : The negative PDF pointer is null." << std::endl;
360  return;
361  }
362  if ( posPdf!=0 ) {
363  std::cerr << "WARNING in LauCPFitModel::setBkgndPdfs : Doing an untagged fit so will not use the positive PDF." << std::endl;
364  }
365  }
366 
367  // check that this background name is valid
368  if ( ! this->validBkgndClass( bkgndClass ) ) {
369  std::cerr << "ERROR in LauCPFitModel::setBkgndPdfs : Invalid background class \"" << bkgndClass << "\"." << std::endl;
370  std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
371  return;
372  }
373 
374  UInt_t bkgndID = this->bkgndClassID( bkgndClass );
375  negBkgndPdfs_[bkgndID].push_back(negPdf);
376  posBkgndPdfs_[bkgndID].push_back(posPdf);
377 
378  usingBkgnd_ = kTRUE;
379 }
380 
382 {
383  // Is there a component called compName in the signal model?
384  TString compName(coeffSet->name());
385  Bool_t negOK = negSigModel_->hasResonance(compName);
386  TString conjName = negSigModel_->getConjResName(compName);
387  Bool_t posOK = posSigModel_->hasResonance(conjName);
388  if (!negOK) {
389  std::cerr << "ERROR in LauCPFitModel::setMagPhase : " << negParent_ << " signal DP model doesn't contain component \"" << compName << "\"." << std::endl;
390  return;
391  }
392  if (!posOK) {
393  std::cerr << "ERROR in LauCPFitModel::setMagPhase : " << posParent_ << " signal DP model doesn't contain component \"" << conjName << "\"." << std::endl;
394  return;
395  }
396 
397  // Do we already have it in our list of names?
398  for (std::vector<LauAbsCoeffSet*>::const_iterator iter=coeffPars_.begin(); iter!=coeffPars_.end(); ++iter) {
399  if ((*iter)->name() == compName) {
400  std::cerr << "ERROR in LauCPFitModel::setAmpCoeffSet : Have already set coefficients for \"" << compName << "\"." << std::endl;
401  return;
402  }
403  }
404 
405  coeffSet->index(nSigComp_);
406  coeffPars_.push_back(coeffSet);
407 
408  TString parName = coeffSet->baseName(); parName += "FitFracAsym";
409  fitFracAsymm_.push_back(LauParameter(parName, 0.0, -1.0, 1.0));
410 
411  acp_.push_back(coeffSet->acp());
412 
413  ++nSigComp_;
414 
415  std::cout << "INFO in LauCPFitModel::setAmpCoeffSet : Added coefficients for component \"" << compName << "\" to the fit model." << std::endl;
416  coeffSet->printParValues();
417 }
418 
420 {
421  // From the initial parameter values calculate the coefficients
422  // so they can be passed to the signal model
423  this->updateCoeffs();
424 
425  // Initialisation
426  if (this->useDP() == kTRUE) {
427  this->initialiseDPModels();
428  }
429 
430  if (!this->useDP() && negSignalPdfs_.empty()) {
431  std::cerr << "ERROR in LauCPFitModel::initialise : Signal model doesn't exist for any variable." << std::endl;
432  gSystem->Exit(EXIT_FAILURE);
433  }
434 
435  if ( this->useDP() ) {
436  // Check that we have all the Dalitz-plot models
437  if ((negSigModel_ == 0) || (posSigModel_ == 0)) {
438  std::cerr << "ERROR in LauCPFitModel::initialise : the pointer to one (neg or pos) of the signal DP models is null.\n";
439  std::cerr << " : Removing the Dalitz Plot from the model." << std::endl;
440  this->useDP(kFALSE);
441  }
442  if ( usingBkgnd_ ) {
443  if ( negBkgndDPModels_.empty() || posBkgndDPModels_.empty() ) {
444  std::cerr << "ERROR in LauCPFitModel::initialise : No background DP models found.\n";
445  std::cerr << " : Removing the Dalitz plot from the model." << std::endl;
446  this->useDP(kFALSE);
447  }
448  for (LauBkgndDPModelList::const_iterator dpmodel_iter = negBkgndDPModels_.begin(); dpmodel_iter != negBkgndDPModels_.end(); ++dpmodel_iter ) {
449  if ( (*dpmodel_iter) == 0 ) {
450  std::cerr << "ERROR in LauCPFitModel::initialise : The pointer to one of the background DP models is null.\n";
451  std::cerr << " : Removing the Dalitz Plot from the model." << std::endl;
452  this->useDP(kFALSE);
453  break;
454  }
455  }
456  for (LauBkgndDPModelList::const_iterator dpmodel_iter = posBkgndDPModels_.begin(); dpmodel_iter != posBkgndDPModels_.end(); ++dpmodel_iter ) {
457  if ( (*dpmodel_iter) == 0 ) {
458  std::cerr << "ERROR in LauCPFitModel::initialise : The pointer to one of the background DP models is null.\n";
459  std::cerr << " : Removing the Dalitz Plot from the model." << std::endl;
460  this->useDP(kFALSE);
461  break;
462  }
463  }
464  }
465  }
466 
467  // Next check that, if a given component is being used we've got the
468  // right number of PDFs for all the variables involved
469  // TODO - should probably check variable names and so on as well
470 
471  UInt_t nsigpdfvars(0);
472  for ( LauPdfList::const_iterator pdf_iter = negSignalPdfs_.begin(); pdf_iter != negSignalPdfs_.end(); ++pdf_iter ) {
473  std::vector<TString> varNames = (*pdf_iter)->varNames();
474  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
475  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
476  ++nsigpdfvars;
477  }
478  }
479  }
480  if (useSCF_) {
481  UInt_t nscfpdfvars(0);
482  for ( LauPdfList::const_iterator pdf_iter = negScfPdfs_.begin(); pdf_iter != negScfPdfs_.end(); ++pdf_iter ) {
483  std::vector<TString> varNames = (*pdf_iter)->varNames();
484  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
485  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
486  ++nscfpdfvars;
487  }
488  }
489  }
490  if (nscfpdfvars != nsigpdfvars) {
491  std::cerr << "ERROR in LauCPFitModel::initialise : There are " << nsigpdfvars << " TM signal PDF variables but " << nscfpdfvars << " SCF signal PDF variables." << std::endl;
492  gSystem->Exit(EXIT_FAILURE);
493  }
494  }
495  if (usingBkgnd_) {
496  for (LauBkgndPdfsList::const_iterator bgclass_iter = negBkgndPdfs_.begin(); bgclass_iter != negBkgndPdfs_.end(); ++bgclass_iter) {
497  UInt_t nbkgndpdfvars(0);
498  const LauPdfList& pdfList = (*bgclass_iter);
499  for ( LauPdfList::const_iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter ) {
500  std::vector<TString> varNames = (*pdf_iter)->varNames();
501  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
502  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
503  ++nbkgndpdfvars;
504  }
505  }
506  }
507  if (nbkgndpdfvars != nsigpdfvars) {
508  std::cerr << "ERROR in LauCPFitModel::initialise : There are " << nsigpdfvars << " signal PDF variables but " << nbkgndpdfvars << " bkgnd PDF variables." << std::endl;
509  gSystem->Exit(EXIT_FAILURE);
510  }
511  }
512  }
513 
514  // Clear the vectors of parameter information so we can start from scratch
515  this->clearFitParVectors();
516 
517  // Set the fit parameters for signal and background models
518  this->setSignalDPParameters();
519 
520  // Set the fit parameters for the various extra PDFs
521  this->setExtraPdfParameters();
522 
523  // Set the initial bg and signal events
524  this->setFitNEvents();
525 
526  // Check that we have the expected number of fit variables
527  const LauParameterPList& fitVars = this->fitPars();
528  if (fitVars.size() != (nSigDPPar_ + nExtraPdfPar_ + nNormPar_)) {
529  std::cerr << "ERROR in LauCPFitModel::initialise : Number of fit parameters not of expected size. Exiting" << std::endl;
530  gSystem->Exit(EXIT_FAILURE);
531  }
532 
533  this->setExtraNtupleVars();
534 }
535 
537 {
538  //std::cout << "INFO in LauCPFitModel::recalculateNormalizationInDPModels : Recalc Norm in DP model" << std::endl;
543 }
544 
546 {
547  // Need to check that the number of components we have and that the dynamics has matches up
548  UInt_t nNegAmp = negSigModel_->getnTotAmp();
549  UInt_t nPosAmp = posSigModel_->getnTotAmp();
550  if ( nNegAmp != nPosAmp ) {
551  std::cerr << "ERROR in LauCPFitModel::initialiseDPModels : Unequal number of signal DP components in the negative and positive models: " << nNegAmp << " != " << nPosAmp << std::endl;
552  gSystem->Exit(EXIT_FAILURE);
553  }
554  if ( nNegAmp != nSigComp_ ) {
555  std::cerr << "ERROR in LauCPFitModel::initialiseDPModels : Number of signal DP components in the model (" << nNegAmp << ") not equal to number of coefficients supplied (" << nSigComp_ << ")." << std::endl;
556  gSystem->Exit(EXIT_FAILURE);
557  }
558 
559 
560  std::cout << "INFO in LauCPFitModel::initialiseDPModels : Initialising signal DP model" << std::endl;
563 
564  if (usingBkgnd_ == kTRUE) {
565  for (LauBkgndDPModelList::iterator iter = negBkgndDPModels_.begin(); iter != negBkgndDPModels_.end(); ++iter) {
566  (*iter)->initialise();
567  }
568  for (LauBkgndDPModelList::iterator iter = posBkgndDPModels_.begin(); iter != posBkgndDPModels_.end(); ++iter) {
569  (*iter)->initialise();
570  }
571  }
572 }
573 
575 {
576  // Set the fit parameters for the signal model.
577 
578  nSigDPPar_ = 0;
579  if ( ! this->useDP() ) {
580  return;
581  }
582 
583  std::cout << "INFO in LauCPFitModel::setSignalDPParameters : Setting the initial fit parameters for the signal DP model." << std::endl;
584 
585  // Place isobar coefficient parameters in vector of fit variables
586  LauParameterPList& fitVars = this->fitPars();
587  for (UInt_t i = 0; i < nSigComp_; i++) {
588  LauParameterPList pars = coeffPars_[i]->getParameters();
589  for (LauParameterPList::iterator iter = pars.begin(); iter != pars.end(); ++iter) {
590  if ( !(*iter)->clone() ) {
591  fitVars.push_back(*iter);
592  ++nSigDPPar_;
593  }
594  }
595  }
596 
597  // Obtain the resonance parameters and place them in the vector of fit variables and in a separate vector
598  // Need to make sure that they are unique because some might appear in both DP models
599  LauParameterPSet& resVars = this->resPars();
600  resVars.clear();
601 
604 
605  for ( LauParameterPList::iterator iter = negSigDPPars.begin(); iter != negSigDPPars.end(); ++iter ) {
606  if ( resVars.insert(*iter).second ) {
607  fitVars.push_back(*iter);
608  ++nSigDPPar_;
609  }
610  }
611  for ( LauParameterPList::iterator iter = posSigDPPars.begin(); iter != posSigDPPars.end(); ++iter ) {
612  if ( resVars.insert(*iter).second ) {
613  fitVars.push_back(*iter);
614  ++nSigDPPar_;
615  }
616  }
617 }
618 
620 {
621  // Include all the parameters of the PDF in the fit
622  // NB all of them are passed to the fit, even though some have been fixed through parameter.fixed(kTRUE)
623  // With the new "cloned parameter" scheme only "original" parameters are passed to the fit.
624  // Their clones are updated automatically when the originals are updated.
625 
626  nExtraPdfPar_ = 0;
627 
629  if ( tagged_ ) {
631  }
632 
633  if (useSCF_ == kTRUE) {
635  if ( tagged_ ) {
637  }
638  }
639 
640  if (usingBkgnd_ == kTRUE) {
641  for (LauBkgndPdfsList::iterator iter = negBkgndPdfs_.begin(); iter != negBkgndPdfs_.end(); ++iter) {
642  nExtraPdfPar_ += this->addFitParameters(*iter);
643  }
644  if ( tagged_ ) {
645  for (LauBkgndPdfsList::iterator iter = posBkgndPdfs_.begin(); iter != posBkgndPdfs_.end(); ++iter) {
646  nExtraPdfPar_ += this->addFitParameters(*iter);
647  }
648  }
649  }
650 }
651 
653 {
654  if ( signalEvents_ == 0 ) {
655  std::cerr << "ERROR in LauCPFitModel::setFitNEvents : Signal yield not defined." << std::endl;
656  return;
657  }
658  nNormPar_ = 0;
659 
660  // initialise the total number of events to be the sum of all the hypotheses
661  Double_t nTotEvts = signalEvents_->value();
662  for (LauBkgndYieldList::const_iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
663  nTotEvts += (*iter)->value();
664  if ( (*iter) == 0 ) {
665  std::cerr << "ERROR in LauCPFitModel::setFitNEvents : Background yield not defined." << std::endl;
666  return;
667  }
668  }
669  this->eventsPerExpt(TMath::FloorNint(nTotEvts));
670 
671  LauParameterPList& fitVars = this->fitPars();
672 
673  // if doing an extended ML fit add the number of signal events into the fit parameters
674  if (this->doEMLFit()) {
675  std::cout << "INFO in LauCPFitModel::setFitNEvents : Initialising number of events for signal and background components..." << std::endl;
676  // add the signal fraction to the list of fit parameters
677  fitVars.push_back(signalEvents_);
678  ++nNormPar_;
679  } else {
680  std::cout << "INFO in LauCPFitModel::setFitNEvents : Initialising number of events for background components (and hence signal)..." << std::endl;
681  }
682 
683  // if not using the DP in the model we need an explicit signal asymmetry parameter
684  if (this->useDP() == kFALSE) {
685  fitVars.push_back(signalAsym_);
686  ++nNormPar_;
687  }
688 
689  if (useSCF_ && !useSCFHist_) {
690  fitVars.push_back(&scfFrac_);
691  ++nNormPar_;
692  }
693 
694  if (usingBkgnd_ == kTRUE) {
695  for (LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
696  LauParameter* parameter = (*iter);
697  fitVars.push_back(parameter);
698  ++nNormPar_;
699  }
700  for (LauBkgndYieldList::iterator iter = bkgndAsym_.begin(); iter != bkgndAsym_.end(); ++iter) {
701  LauParameter* parameter = (*iter);
702  fitVars.push_back(parameter);
703  ++nNormPar_;
704  }
705  }
706 }
707 
709 {
710  // Set-up other parameters derived from the fit results, e.g. fit fractions.
711 
712  if (this->useDP() != kTRUE) {
713  return;
714  }
715 
716  // First clear the vectors so we start from scratch
717  this->clearExtraVarVectors();
718 
719  LauParameterList& extraVars = this->extraPars();
720 
721  // Add the positive and negative fit fractions for each signal component
723  if (negFitFrac_.size() != nSigComp_) {
724  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << negFitFrac_.size() << std::endl;
725  gSystem->Exit(EXIT_FAILURE);
726  }
727  for (UInt_t i(0); i<nSigComp_; ++i) {
728  if (negFitFrac_[i].size() != nSigComp_) {
729  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << negFitFrac_[i].size() << std::endl;
730  gSystem->Exit(EXIT_FAILURE);
731  }
732  }
733 
735  if (posFitFrac_.size() != nSigComp_) {
736  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << posFitFrac_.size() << std::endl;
737  gSystem->Exit(EXIT_FAILURE);
738  }
739  for (UInt_t i(0); i<nSigComp_; ++i) {
740  if (posFitFrac_[i].size() != nSigComp_) {
741  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << posFitFrac_[i].size() << std::endl;
742  gSystem->Exit(EXIT_FAILURE);
743  }
744  }
745 
746  // Add the positive and negative fit fractions that have not been corrected for the efficiency for each signal component
748  if (negFitFracEffUnCorr_.size() != nSigComp_) {
749  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << negFitFracEffUnCorr_.size() << std::endl;
750  gSystem->Exit(EXIT_FAILURE);
751  }
752  for (UInt_t i(0); i<nSigComp_; ++i) {
753  if (negFitFracEffUnCorr_[i].size() != nSigComp_) {
754  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << negFitFracEffUnCorr_[i].size() << std::endl;
755  gSystem->Exit(EXIT_FAILURE);
756  }
757  }
758 
760  if (posFitFracEffUnCorr_.size() != nSigComp_) {
761  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << posFitFracEffUnCorr_.size() << std::endl;
762  gSystem->Exit(EXIT_FAILURE);
763  }
764  for (UInt_t i(0); i<nSigComp_; ++i) {
765  if (posFitFracEffUnCorr_[i].size() != nSigComp_) {
766  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << posFitFracEffUnCorr_[i].size() << std::endl;
767  gSystem->Exit(EXIT_FAILURE);
768  }
769  }
770 
771  for (UInt_t i(0); i<nSigComp_; ++i) {
772  for (UInt_t j(i); j<nSigComp_; ++j) {
773  TString name = negFitFrac_[i][j].name();
774  name.Insert( name.Index("FitFrac"), "Neg" );
775  negFitFrac_[i][j].name(name);
776  extraVars.push_back(negFitFrac_[i][j]);
777 
778  name = posFitFrac_[i][j].name();
779  name.Insert( name.Index("FitFrac"), "Pos" );
780  posFitFrac_[i][j].name(name);
781  extraVars.push_back(posFitFrac_[i][j]);
782 
783  name = negFitFracEffUnCorr_[i][j].name();
784  name.Insert( name.Index("FitFrac"), "Neg" );
785  negFitFracEffUnCorr_[i][j].name(name);
786  extraVars.push_back(negFitFracEffUnCorr_[i][j]);
787 
788  name = posFitFracEffUnCorr_[i][j].name();
789  name.Insert( name.Index("FitFrac"), "Pos" );
790  posFitFracEffUnCorr_[i][j].name(name);
791  extraVars.push_back(posFitFracEffUnCorr_[i][j]);
792  }
793  }
794 
795  // Store any extra parameters/quantities from the DP model (e.g. K-matrix total fit fractions)
796  std::vector<LauParameter> negExtraPars = negSigModel_->getExtraParameters();
797  std::vector<LauParameter>::iterator negExtraIter;
798  for (negExtraIter = negExtraPars.begin(); negExtraIter != negExtraPars.end(); ++negExtraIter) {
799  LauParameter negExtraParameter = (*negExtraIter);
800  extraVars.push_back(negExtraParameter);
801  }
802 
803  std::vector<LauParameter> posExtraPars = posSigModel_->getExtraParameters();
804  std::vector<LauParameter>::iterator posExtraIter;
805  for (posExtraIter = posExtraPars.begin(); posExtraIter != posExtraPars.end(); ++posExtraIter) {
806  LauParameter posExtraParameter = (*posExtraIter);
807  extraVars.push_back(posExtraParameter);
808  }
809 
810  // Now add in the DP efficiency value
811  Double_t initMeanEff = negSigModel_->getMeanEff().initValue();
812  negMeanEff_.value(initMeanEff);
813  negMeanEff_.genValue(initMeanEff);
814  negMeanEff_.initValue(initMeanEff);
815  extraVars.push_back(negMeanEff_);
816 
817  initMeanEff = posSigModel_->getMeanEff().initValue();
818  posMeanEff_.value(initMeanEff);
819  posMeanEff_.genValue(initMeanEff);
820  posMeanEff_.initValue(initMeanEff);
821  extraVars.push_back(posMeanEff_);
822 
823  // Also add in the DP rates
824  Double_t initDPRate = negSigModel_->getDPRate().initValue();
825  negDPRate_.value(initDPRate);
826  negDPRate_.genValue(initDPRate);
827  negDPRate_.initValue(initDPRate);
828  extraVars.push_back(negDPRate_);
829 
830  initDPRate = posSigModel_->getDPRate().initValue();
831  posDPRate_.value(initDPRate);
832  posDPRate_.genValue(initDPRate);
833  posDPRate_.initValue(initDPRate);
834  extraVars.push_back(posDPRate_);
835 
836  // Calculate the CPC and CPV Fit Fractions, ACPs and FitFrac asymmetries
837  this->calcExtraFractions(kTRUE);
838  this->calcAsymmetries(kTRUE);
839 
840  // Add the CP violating and CP conserving fit fractions for each signal component
841  for (UInt_t i = 0; i < nSigComp_; i++) {
842  for (UInt_t j = i; j < nSigComp_; j++) {
843  extraVars.push_back(CPVFitFrac_[i][j]);
844  }
845  }
846  for (UInt_t i = 0; i < nSigComp_; i++) {
847  for (UInt_t j = i; j < nSigComp_; j++) {
848  extraVars.push_back(CPCFitFrac_[i][j]);
849  }
850  }
851 
852  // Add the Fit Fraction asymmetry for each signal component
853  for (UInt_t i = 0; i < nSigComp_; i++) {
854  extraVars.push_back(fitFracAsymm_[i]);
855  }
856 
857  // Add the calculated CP asymmetry for each signal component
858  for (UInt_t i = 0; i < nSigComp_; i++) {
859  extraVars.push_back(acp_[i]);
860  }
861 }
862 
863 void LauCPFitModel::calcExtraFractions(Bool_t initValues)
864 {
865  // Calculate the CP-conserving and CP-violating fit fractions
866 
867  if (initValues) {
868  // create the structure
869  CPCFitFrac_.clear();
870  CPVFitFrac_.clear();
871  CPCFitFrac_.resize(nSigComp_);
872  CPVFitFrac_.resize(nSigComp_);
873  for (UInt_t i(0); i<nSigComp_; ++i) {
874  CPCFitFrac_[i].resize(nSigComp_);
875  CPVFitFrac_[i].resize(nSigComp_);
876 
877  for (UInt_t j(i); j<nSigComp_; ++j) {
878  TString name = negFitFrac_[i][j].name();
879  name.Replace( name.Index("Neg"), 3, "CPC" );
880  CPCFitFrac_[i][j].name( name );
881  CPCFitFrac_[i][j].valueAndRange( 0.0, -100.0, 100.0 );
882 
883  name = negFitFrac_[i][j].name();
884  name.Replace( name.Index("Neg"), 3, "CPV" );
885  CPVFitFrac_[i][j].name( name );
886  CPVFitFrac_[i][j].valueAndRange( 0.0, -100.0, 100.0 );
887  }
888  }
889  }
890 
891  Double_t denom = negDPRate_ + posDPRate_;
892 
893  for (UInt_t i(0); i<nSigComp_; ++i) {
894  for (UInt_t j(i); j<nSigComp_; ++j) {
895 
896  Double_t negTerm = negFitFrac_[i][j]*negDPRate_;
897  Double_t posTerm = posFitFrac_[i][j]*posDPRate_;
898 
899  Double_t cpcFitFrac = (negTerm + posTerm)/denom;
900  Double_t cpvFitFrac = (negTerm - posTerm)/denom;
901 
902  CPCFitFrac_[i][j] = cpcFitFrac;
903  CPVFitFrac_[i][j] = cpvFitFrac;
904 
905  if (initValues) {
906  CPCFitFrac_[i][j].genValue(cpcFitFrac);
907  CPCFitFrac_[i][j].initValue(cpcFitFrac);
908 
909  CPVFitFrac_[i][j].genValue(cpvFitFrac);
910  CPVFitFrac_[i][j].initValue(cpvFitFrac);
911  }
912  }
913  }
914 }
915 
916 void LauCPFitModel::calcAsymmetries(Bool_t initValues)
917 {
918  // Calculate the CP asymmetries
919  // Also calculate the fit fraction asymmetries
920 
921  for (UInt_t i = 0; i < nSigComp_; i++) {
922 
923  acp_[i] = coeffPars_[i]->acp();
924 
925  LauAsymmCalc asymmCalc(negFitFrac_[i][i].value(), posFitFrac_[i][i].value());
926  Double_t asym = asymmCalc.getAsymmetry();
927  fitFracAsymm_[i] = asym;
928  if (initValues) {
929  fitFracAsymm_[i].genValue(asym);
930  fitFracAsymm_[i].initValue(asym);
931  }
932  }
933 }
934 
935 void LauCPFitModel::finaliseFitResults(const TString& tablePrefixName)
936 {
937  // Retrieve parameters from the fit results for calculations and toy generation
938  // and eventually store these in output root ntuples/text files
939 
940  // Now take the fit parameters and update them as necessary
941  // i.e. to make mag > 0.0, phase in the right range.
942  // This function will also calculate any other values, such as the
943  // fit fractions, using any errors provided by fitParErrors as appropriate.
944  // Also obtain the pull values: (measured - generated)/(average error)
945 
946  if (this->useDP() == kTRUE) {
947  for (UInt_t i = 0; i < nSigComp_; ++i) {
948  // Check whether we have "a/b > 0.0", and phases in the right range
949  coeffPars_[i]->finaliseValues();
950  }
951  }
952 
953  // update the pulls on the event fractions and asymmetries
954  if (this->doEMLFit()) {
956  }
957  if (this->useDP() == kFALSE) {
959  }
960  if (useSCF_ && !useSCFHist_) {
962  }
963  if (usingBkgnd_ == kTRUE) {
964  for (LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
965  (*iter)->updatePull();
966  }
967  for (LauBkgndYieldList::iterator iter = bkgndAsym_.begin(); iter != bkgndAsym_.end(); ++iter) {
968  (*iter)->updatePull();
969  }
970  }
971 
972  // Update the pulls on all the extra PDFs' parameters
975  if (useSCF_ == kTRUE) {
978  }
979  if (usingBkgnd_ == kTRUE) {
980  for (LauBkgndPdfsList::iterator iter = negBkgndPdfs_.begin(); iter != negBkgndPdfs_.end(); ++iter) {
981  this->updateFitParameters(*iter);
982  }
983  for (LauBkgndPdfsList::iterator iter = posBkgndPdfs_.begin(); iter != posBkgndPdfs_.end(); ++iter) {
984  this->updateFitParameters(*iter);
985  }
986  }
987 
988  // Fill the fit results to the ntuple
989 
990  // update the coefficients and then calculate the fit fractions and ACP's
991  if (this->useDP() == kTRUE) {
992  this->updateCoeffs();
995 
996  LauParArray negFitFrac = negSigModel_->getFitFractions();
997  if (negFitFrac.size() != nSigComp_) {
998  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << negFitFrac.size() << std::endl;
999  gSystem->Exit(EXIT_FAILURE);
1000  }
1001  for (UInt_t i(0); i<nSigComp_; ++i) {
1002  if (negFitFrac[i].size() != nSigComp_) {
1003  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << negFitFrac[i].size() << std::endl;
1004  gSystem->Exit(EXIT_FAILURE);
1005  }
1006  }
1007 
1008  LauParArray posFitFrac = posSigModel_->getFitFractions();
1009  if (posFitFrac.size() != nSigComp_) {
1010  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << posFitFrac.size() << std::endl;
1011  gSystem->Exit(EXIT_FAILURE);
1012  }
1013  for (UInt_t i(0); i<nSigComp_; ++i) {
1014  if (posFitFrac[i].size() != nSigComp_) {
1015  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << posFitFrac[i].size() << std::endl;
1016  gSystem->Exit(EXIT_FAILURE);
1017  }
1018  }
1019 
1021  if (negFitFracEffUnCorr.size() != nSigComp_) {
1022  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << negFitFracEffUnCorr.size() << std::endl;
1023  gSystem->Exit(EXIT_FAILURE);
1024  }
1025  for (UInt_t i(0); i<nSigComp_; ++i) {
1026  if (negFitFracEffUnCorr[i].size() != nSigComp_) {
1027  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << negFitFracEffUnCorr[i].size() << std::endl;
1028  gSystem->Exit(EXIT_FAILURE);
1029  }
1030  }
1031 
1033  if (posFitFracEffUnCorr.size() != nSigComp_) {
1034  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << posFitFracEffUnCorr.size() << std::endl;
1035  gSystem->Exit(EXIT_FAILURE);
1036  }
1037  for (UInt_t i(0); i<nSigComp_; ++i) {
1038  if (posFitFracEffUnCorr[i].size() != nSigComp_) {
1039  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << posFitFracEffUnCorr[i].size() << std::endl;
1040  gSystem->Exit(EXIT_FAILURE);
1041  }
1042  }
1043 
1044  for (UInt_t i(0); i<nSigComp_; ++i) {
1045  for (UInt_t j(i); j<nSigComp_; ++j) {
1046  negFitFrac_[i][j].value(negFitFrac[i][j].value());
1047  posFitFrac_[i][j].value(posFitFrac[i][j].value());
1048  negFitFracEffUnCorr_[i][j].value(negFitFracEffUnCorr[i][j].value());
1049  posFitFracEffUnCorr_[i][j].value(posFitFracEffUnCorr[i][j].value());
1050  }
1051  }
1052 
1057 
1058  this->calcExtraFractions();
1059  this->calcAsymmetries();
1060 
1061  // Then store the final fit parameters, and any extra parameters for
1062  // the signal model (e.g. fit fractions, FF asymmetries, ACPs, mean efficiency and DP rate)
1063 
1064  this->clearExtraVarVectors();
1065  LauParameterList& extraVars = this->extraPars();
1066 
1067  // Add the positive and negative fit fractions for each signal component
1068  for (UInt_t i(0); i<nSigComp_; ++i) {
1069  for (UInt_t j(i); j<nSigComp_; ++j) {
1070  extraVars.push_back(negFitFrac_[i][j]);
1071  extraVars.push_back(posFitFrac_[i][j]);
1072  extraVars.push_back(negFitFracEffUnCorr_[i][j]);
1073  extraVars.push_back(posFitFracEffUnCorr_[i][j]);
1074  }
1075  }
1076 
1077  // Store any extra parameters/quantities from the DP model (e.g. K-matrix total fit fractions)
1078  std::vector<LauParameter> negExtraPars = negSigModel_->getExtraParameters();
1079  std::vector<LauParameter>::iterator negExtraIter;
1080  for (negExtraIter = negExtraPars.begin(); negExtraIter != negExtraPars.end(); ++negExtraIter) {
1081  LauParameter negExtraParameter = (*negExtraIter);
1082  extraVars.push_back(negExtraParameter);
1083  }
1084 
1085  std::vector<LauParameter> posExtraPars = posSigModel_->getExtraParameters();
1086  std::vector<LauParameter>::iterator posExtraIter;
1087  for (posExtraIter = posExtraPars.begin(); posExtraIter != posExtraPars.end(); ++posExtraIter) {
1088  LauParameter posExtraParameter = (*posExtraIter);
1089  extraVars.push_back(posExtraParameter);
1090  }
1091 
1092  extraVars.push_back(negMeanEff_);
1093  extraVars.push_back(posMeanEff_);
1094  extraVars.push_back(negDPRate_);
1095  extraVars.push_back(posDPRate_);
1096 
1097  for (UInt_t i = 0; i < nSigComp_; i++) {
1098  for (UInt_t j(i); j<nSigComp_; ++j) {
1099  extraVars.push_back(CPVFitFrac_[i][j]);
1100  }
1101  }
1102  for (UInt_t i = 0; i < nSigComp_; i++) {
1103  for (UInt_t j(i); j<nSigComp_; ++j) {
1104  extraVars.push_back(CPCFitFrac_[i][j]);
1105  }
1106  }
1107  for (UInt_t i(0); i<nSigComp_; ++i) {
1108  extraVars.push_back(fitFracAsymm_[i]);
1109  }
1110  for (UInt_t i(0); i<nSigComp_; ++i) {
1111  extraVars.push_back(acp_[i]);
1112  }
1113 
1114  this->printFitFractions(std::cout);
1115  this->printAsymmetries(std::cout);
1116  }
1117 
1118  const LauParameterPList& fitVars = this->fitPars();
1119  const LauParameterList& extraVars = this->extraPars();
1120  LauFitNtuple* ntuple = this->fitNtuple();
1121  ntuple->storeParsAndErrors(fitVars, extraVars);
1122 
1123  // find out the correlation matrix for the parameters
1124  ntuple->storeCorrMatrix(this->iExpt(), this->nll(), this->fitStatus(), this->covarianceMatrix());
1125 
1126  // Fill the data into ntuple
1127  ntuple->updateFitNtuple();
1128 
1129  // Print out the partial fit fractions, phases and the
1130  // averaged efficiency, reweighted by the dynamics (and anything else)
1131  if (this->writeLatexTable()) {
1132  TString sigOutFileName(tablePrefixName);
1133  sigOutFileName += "_"; sigOutFileName += this->iExpt(); sigOutFileName += "Expt.tex";
1134  this->writeOutTable(sigOutFileName);
1135  }
1136 }
1137 
1138 void LauCPFitModel::printFitFractions(std::ostream& output)
1139 {
1140  // Print out Fit Fractions, total DP rate and mean efficiency
1141  // First for the B- events
1142  for (UInt_t i = 0; i < nSigComp_; i++) {
1143  const TString compName(coeffPars_[i]->name());
1144  output << negParent_ << " FitFraction for component " << i << " (" << compName << ") = " << negFitFrac_[i][i] << std::endl;
1145  }
1146  output << negParent_ << " overall DP rate (integral of matrix element squared) = " << negDPRate_ << std::endl;
1147  output << negParent_ << " average efficiency weighted by whole DP dynamics = " << negMeanEff_ << std::endl;
1148 
1149  // Then for the positive sample
1150  for (UInt_t i = 0; i < nSigComp_; i++) {
1151  const TString compName(coeffPars_[i]->name());
1152  const TString conjName(negSigModel_->getConjResName(compName));
1153  output << posParent_ << " FitFraction for component " << i << " (" << conjName << ") = " << posFitFrac_[i][i] << std::endl;
1154  }
1155  output << posParent_ << " overall DP rate (integral of matrix element squared) = " << posDPRate_ << std::endl;
1156  output << posParent_ << " average efficiency weighted by whole DP dynamics = " << posMeanEff_ << std::endl;
1157 }
1158 
1159 void LauCPFitModel::printAsymmetries(std::ostream& output)
1160 {
1161  for (UInt_t i = 0; i < nSigComp_; i++) {
1162  const TString compName(coeffPars_[i]->name());
1163  output << "Fit Fraction asymmetry for component " << i << " (" << compName << ") = " << fitFracAsymm_[i] << std::endl;
1164  }
1165  for (UInt_t i = 0; i < nSigComp_; i++) {
1166  const TString compName(coeffPars_[i]->name());
1167  output << "ACP for component " << i << " (" << compName << ") = " << acp_[i] << std::endl;
1168  }
1169 }
1170 
1171 void LauCPFitModel::writeOutTable(const TString& outputFile)
1172 {
1173  // Write out the results of the fit to a tex-readable table
1174  // TODO - need to include the yields in this table
1175  std::ofstream fout(outputFile);
1176  LauPrint print;
1177 
1178  std::cout << "INFO in LauCPFitModel::writeOutTable : Writing out results of the fit to the tex file " << outputFile << std::endl;
1179 
1180  if (this->useDP() == kTRUE) {
1181  // print the fit coefficients in one table
1182  coeffPars_.front()->printTableHeading(fout);
1183  for (UInt_t i = 0; i < nSigComp_; i++) {
1184  coeffPars_[i]->printTableRow(fout);
1185  }
1186  fout << "\\hline" << std::endl;
1187  fout << "\\end{tabular}" << std::endl << std::endl;
1188 
1189  // print the fit fractions and asymmetries in another
1190  fout << "\\begin{tabular}{|l|c|c|c|c|}" << std::endl;
1191  fout << "\\hline" << std::endl;
1192  fout << "Component & " << negParent_ << " Fit Fraction & " << posParent_ << " Fit Fraction & Fit Fraction Asymmetry & ACP \\\\" << std::endl;
1193  fout << "\\hline" << std::endl;
1194 
1195  Double_t negFitFracSum(0.0);
1196  Double_t posFitFracSum(0.0);
1197 
1198  for (UInt_t i = 0; i < nSigComp_; i++) {
1199  TString resName = coeffPars_[i]->name();
1200  resName = resName.ReplaceAll("_", "\\_");
1201 
1202  Double_t negFitFrac = negFitFrac_[i][i].value();
1203  Double_t posFitFrac = posFitFrac_[i][i].value();
1204  negFitFracSum += negFitFrac;
1205  posFitFracSum += posFitFrac;
1206 
1207  Double_t fitFracAsymm = fitFracAsymm_[i].value();
1208 
1209  Double_t acp = acp_[i].value();
1210  Double_t acpErr = acp_[i].error();
1211 
1212  fout << resName << " & $";
1213 
1214  print.printFormat(fout, negFitFrac);
1215  fout << "$ & $";
1216  print.printFormat(fout, posFitFrac);
1217  fout << "$ & $";
1218 
1219  print.printFormat(fout, fitFracAsymm);
1220  fout << "$ & $";
1221 
1222  print.printFormat(fout, acp);
1223  fout << " \\pm ";
1224  print.printFormat(fout, acpErr);
1225  fout << "$ \\\\" << std::endl;
1226  }
1227  fout << "\\hline" << std::endl;
1228 
1229  // Also print out sum of fit fractions
1230  fout << "Fit Fraction Sum & $";
1231  print.printFormat(fout, negFitFracSum);
1232  fout << "$ & $";
1233  print.printFormat(fout, posFitFracSum);
1234  fout << "$ & & \\\\" << std::endl;
1235  fout << "\\hline" << std::endl;
1236 
1237  fout << "DP rate & $";
1238  print.printFormat(fout, negDPRate_.value());
1239  fout << "$ & $";
1240  print.printFormat(fout, posDPRate_.value());
1241  fout << "$ & & \\\\" << std::endl;
1242 
1243  fout << "$< \\varepsilon > $ & $";
1244  print.printFormat(fout, negMeanEff_.value());
1245  fout << "$ & $";
1246  print.printFormat(fout, posMeanEff_.value());
1247  fout << "$ & & \\\\" << std::endl;
1248  fout << "\\hline" << std::endl;
1249  fout << "\\end{tabular}" << std::endl << std::endl;
1250  }
1251 
1252  if (!negSignalPdfs_.empty()) {
1253  fout << "\\begin{tabular}{|l|c|}" << std::endl;
1254  fout << "\\hline" << std::endl;
1255  if (useSCF_ == kTRUE) {
1256  fout << "\\Extra TM Signal PDFs' Parameters: & \\\\" << std::endl;
1257  } else {
1258  fout << "\\Extra Signal PDFs' Parameters: & \\\\" << std::endl;
1259  }
1260  this->printFitParameters(negSignalPdfs_, fout);
1261  if ( tagged_ ) {
1262  this->printFitParameters(posSignalPdfs_, fout);
1263  }
1264  if (useSCF_ == kTRUE && !negScfPdfs_.empty()) {
1265  fout << "\\hline" << std::endl;
1266  fout << "\\Extra SCF Signal PDFs' Parameters: & \\\\" << std::endl;
1267  this->printFitParameters(negScfPdfs_, fout);
1268  if ( tagged_ ) {
1269  this->printFitParameters(posScfPdfs_, fout);
1270  }
1271  }
1272  if (usingBkgnd_ == kTRUE && !negBkgndPdfs_.empty()) {
1273  fout << "\\hline" << std::endl;
1274  fout << "\\Extra Background PDFs' Parameters: & \\\\" << std::endl;
1275  for (LauBkgndPdfsList::const_iterator iter = negBkgndPdfs_.begin(); iter != negBkgndPdfs_.end(); ++iter) {
1276  this->printFitParameters(*iter, fout);
1277  }
1278  if ( tagged_ ) {
1279  for (LauBkgndPdfsList::const_iterator iter = posBkgndPdfs_.begin(); iter != posBkgndPdfs_.end(); ++iter) {
1280  this->printFitParameters(*iter, fout);
1281  }
1282  }
1283  }
1284  fout << "\\hline \n\\end{tabular}" << std::endl << std::endl;
1285  }
1286 }
1287 
1289 {
1290  // Update the number of signal events to be total-sum(background events)
1291  this->updateSigEvents();
1292 
1293  // Check whether we want to have randomised initial fit parameters for the signal model
1294  if (this->useRandomInitFitPars() == kTRUE) {
1295  std::cout << "INFO in LauCPFitModel::checkInitFitParams : Setting random parameters for the signal model" << std::endl;
1296  this->randomiseInitFitPars();
1297  }
1298 }
1299 
1301 {
1302  // Only randomise those parameters that are not fixed!
1303  std::cout << "INFO in LauCPFitModel::randomiseInitFitPars : Randomising the initial fit magnitudes and phases of the components..." << std::endl;
1304 
1305  for (UInt_t i = 0; i < nSigComp_; i++) {
1306  coeffPars_[i]->randomiseInitValues();
1307  }
1308 }
1309 
1311 {
1312  // Determine the number of events to generate for each hypothesis
1313  // If we're smearing then smear each one individually
1314 
1315  LauGenInfo nEvtsGen;
1316 
1317  // Signal
1318  Double_t evtWeight(1.0);
1319  Double_t nEvts = signalEvents_->genValue();
1320  if ( nEvts < 0.0 ) {
1321  evtWeight = -1.0;
1322  nEvts = TMath::Abs( nEvts );
1323  }
1324  Double_t asym(0.0);
1325  Double_t sigAsym(0.0);
1326  // need to include this as an alternative in case the DP isn't in the model
1327  if ( !this->useDP() || forceAsym_ ) {
1328  sigAsym = signalAsym_->genValue();
1329  } else {
1330  Double_t negRate = negSigModel_->getDPNorm();
1331  Double_t posRate = posSigModel_->getDPNorm();
1332  if (negRate+posRate>1e-30) {
1333  sigAsym = (negRate-posRate)/(negRate+posRate);
1334  }
1335  }
1336  asym = sigAsym;
1337  Int_t nPosEvts = static_cast<Int_t>((nEvts/2.0 * (1.0 - asym)) + 0.5);
1338  Int_t nNegEvts = static_cast<Int_t>((nEvts/2.0 * (1.0 + asym)) + 0.5);
1339  if (this->doPoissonSmearing()) {
1340  nNegEvts = LauRandom::randomFun()->Poisson(nNegEvts);
1341  nPosEvts = LauRandom::randomFun()->Poisson(nPosEvts);
1342  }
1343  nEvtsGen[std::make_pair("signal",-1)] = std::make_pair(nNegEvts,evtWeight);
1344  nEvtsGen[std::make_pair("signal",+1)] = std::make_pair(nPosEvts,evtWeight);
1345 
1346  // backgrounds
1347  const UInt_t nBkgnds = this->nBkgndClasses();
1348  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
1349  const TString& bkgndClass = this->bkgndClassName(bkgndID);
1350  const LauParameter* evtsPar = bkgndEvents_[bkgndID];
1351  const LauParameter* asymPar = bkgndAsym_[bkgndID];
1352  evtWeight = 1.0;
1353  nEvts = TMath::FloorNint( evtsPar->genValue() );
1354  if ( nEvts < 0 ) {
1355  evtWeight = -1.0;
1356  nEvts = TMath::Abs( nEvts );
1357  }
1358  asym = asymPar->genValue();
1359  nPosEvts = static_cast<Int_t>((nEvts/2.0 * (1.0 - asym)) + 0.5);
1360  nNegEvts = static_cast<Int_t>((nEvts/2.0 * (1.0 + asym)) + 0.5);
1361  if (this->doPoissonSmearing()) {
1362  nNegEvts = LauRandom::randomFun()->Poisson(nNegEvts);
1363  nPosEvts = LauRandom::randomFun()->Poisson(nPosEvts);
1364  }
1365  nEvtsGen[std::make_pair(bkgndClass,-1)] = std::make_pair(nNegEvts,evtWeight);
1366  nEvtsGen[std::make_pair(bkgndClass,+1)] = std::make_pair(nPosEvts,evtWeight);
1367  }
1368 
1369  std::cout << "INFO in LauCPFitModel::eventsToGenerate : Generating toy MC with:" << std::endl;
1370  std::cout << " : Signal asymmetry = " << sigAsym << " and number of signal events = " << signalEvents_->genValue() << std::endl;
1371  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
1372  const TString& bkgndClass = this->bkgndClassName(bkgndID);
1373  const LauParameter* evtsPar = bkgndEvents_[bkgndID];
1374  const LauParameter* asymPar = bkgndAsym_[bkgndID];
1375  std::cout << " : " << bkgndClass << " asymmetry = " << asymPar->genValue() << " and number of " << bkgndClass << " events = " << evtsPar->genValue() << std::endl;
1376  }
1377 
1378  return nEvtsGen;
1379 }
1380 
1382 {
1383  // Routine to generate toy Monte Carlo events according to the various models we have defined.
1384 
1385  // Determine the number of events to generate for each hypothesis
1386  LauGenInfo nEvts = this->eventsToGenerate();
1387 
1388  Bool_t genOK(kTRUE);
1389  Int_t evtNum(0);
1390 
1391  const UInt_t nBkgnds = this->nBkgndClasses();
1392  std::vector<TString> bkgndClassNames(nBkgnds);
1393  std::vector<TString> bkgndClassNamesGen(nBkgnds);
1394  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
1395  TString name( this->bkgndClassName(iBkgnd) );
1396  bkgndClassNames[iBkgnd] = name;
1397  bkgndClassNamesGen[iBkgnd] = "gen"+name;
1398  }
1399 
1400  const Bool_t storeSCFTruthInfo = ( useSCF_ || ( this->enableEmbedding() &&
1401  negSignalTree_ != 0 && negSignalTree_->haveBranch("mcMatch") &&
1402  posSignalTree_ != 0 && posSignalTree_->haveBranch("mcMatch") ) );
1403 
1404  // Loop over the hypotheses and generate the requested number of events for each one
1405  for (LauGenInfo::const_iterator iter = nEvts.begin(); iter != nEvts.end(); ++iter) {
1406 
1407  const TString& type(iter->first.first);
1408  curEvtCharge_ = iter->first.second;
1409  Double_t evtWeight( iter->second.second );
1410  Int_t nEvtsGen( iter->second.first );
1411 
1412  for (Int_t iEvt(0); iEvt<nEvtsGen; ++iEvt) {
1413 
1414  this->setGenNtupleDoubleBranchValue( "evtWeight", evtWeight );
1415  this->setGenNtupleDoubleBranchValue( "efficiency", 1.0 );
1416 
1417  if (type == "signal") {
1418  this->setGenNtupleIntegerBranchValue("genSig",1);
1419  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
1420  this->setGenNtupleIntegerBranchValue( bkgndClassNamesGen[iBkgnd], 0 );
1421  }
1422  genOK = this->generateSignalEvent();
1423  if ( curEvtCharge_ > 0 ){
1424  this->setGenNtupleDoubleBranchValue( "efficiency", posSigModel_->getEvtEff() );
1425  } else {
1426  this->setGenNtupleDoubleBranchValue( "efficiency", negSigModel_->getEvtEff() );
1427  }
1428 
1429  } else {
1430  this->setGenNtupleIntegerBranchValue("genSig",0);
1431  if ( storeSCFTruthInfo ) {
1432  this->setGenNtupleIntegerBranchValue("genTMSig",0);
1433  this->setGenNtupleIntegerBranchValue("genSCFSig",0);
1434  }
1435  UInt_t bkgndID(0);
1436  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
1437  Int_t gen(0);
1438  if ( bkgndClassNames[iBkgnd] == type ) {
1439  gen = 1;
1440  bkgndID = iBkgnd;
1441  }
1442  this->setGenNtupleIntegerBranchValue( bkgndClassNamesGen[iBkgnd], gen );
1443  }
1444  genOK = this->generateBkgndEvent(bkgndID);
1445  }
1446  if (!genOK) {
1447  // If there was a problem with the generation then break out and return.
1448  // The problem model will have adjusted itself so that all should be OK next time.
1449  break;
1450  }
1451  if (this->useDP() == kTRUE) {
1452  this->setDPBranchValues();
1453  }
1454 
1455  // Store the event charge
1457 
1458  // Store the event number (within this experiment)
1459  // and then increment it
1460  this->setGenNtupleIntegerBranchValue("iEvtWithinExpt",evtNum);
1461  ++evtNum;
1462 
1463  this->fillGenNtupleBranches();
1464  if (iEvt%500 == 0) {std::cout << "INFO in LauCPFitModel::genExpt : Generated event number " << iEvt << " out of " << nEvtsGen << " " << type << " events." << std::endl;}
1465  }
1466 
1467  if (!genOK) {
1468  break;
1469  }
1470  }
1471 
1472  if (this->useDP() && genOK) {
1473 
1474  negSigModel_->checkToyMC(kTRUE,kTRUE);
1475  posSigModel_->checkToyMC(kTRUE,kTRUE);
1476 
1477  // Get the fit fractions if they're to be written into the latex table
1478  if (this->writeLatexTable()) {
1479  LauParArray negFitFrac = negSigModel_->getFitFractions();
1480  if (negFitFrac.size() != nSigComp_) {
1481  std::cerr << "ERROR in LauCPFitModel::genExpt : Fit Fraction array of unexpected dimension: " << negFitFrac.size() << std::endl;
1482  gSystem->Exit(EXIT_FAILURE);
1483  }
1484  for (UInt_t i(0); i<nSigComp_; ++i) {
1485  if (negFitFrac[i].size() != nSigComp_) {
1486  std::cerr << "ERROR in LauCPFitModel::genExpt : Fit Fraction array of unexpected dimension: " << negFitFrac[i].size() << std::endl;
1487  gSystem->Exit(EXIT_FAILURE);
1488  }
1489  }
1490  LauParArray posFitFrac = posSigModel_->getFitFractions();
1491  if (posFitFrac.size() != nSigComp_) {
1492  std::cerr << "ERROR in LauCPFitModel::genExpt : Fit Fraction array of unexpected dimension: " << posFitFrac.size() << std::endl;
1493  gSystem->Exit(EXIT_FAILURE);
1494  }
1495  for (UInt_t i(0); i<nSigComp_; ++i) {
1496  if (posFitFrac[i].size() != nSigComp_) {
1497  std::cerr << "ERROR in LauCPFitModel::genExpt : Fit Fraction array of unexpected dimension: " << posFitFrac[i].size() << std::endl;
1498  gSystem->Exit(EXIT_FAILURE);
1499  }
1500  }
1501 
1502  for (UInt_t i(0); i<nSigComp_; ++i) {
1503  for (UInt_t j(i); j<nSigComp_; ++j) {
1504  negFitFrac_[i][j].value(negFitFrac[i][j].value());
1505  posFitFrac_[i][j].value(posFitFrac[i][j].value());
1506  }
1507  }
1512  }
1513  }
1514 
1515  // If we're reusing embedded events or if the generation is being
1516  // reset then clear the lists of used events
1517  if (reuseSignal_ || !genOK) {
1518  if (negSignalTree_) {
1520  }
1521  if (posSignalTree_) {
1523  }
1524  }
1525  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
1526  LauEmbeddedData* data = negBkgndTree_[bkgndID];
1527  if (reuseBkgnd_[bkgndID] || !genOK) {
1528  if (data) {
1529  data->clearUsedList();
1530  }
1531  }
1532  }
1533  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
1534  LauEmbeddedData* data = posBkgndTree_[bkgndID];
1535  if (reuseBkgnd_[bkgndID] || !genOK) {
1536  if (data) {
1537  data->clearUsedList();
1538  }
1539  }
1540  }
1541 
1542  return genOK;
1543 }
1544 
1546 {
1547  // Generate signal event
1548  Bool_t genOK(kTRUE);
1549  Bool_t genSCF(kFALSE);
1550 
1551  LauIsobarDynamics* model(0);
1552  LauKinematics* kinematics(0);
1553  LauEmbeddedData* embeddedData(0);
1554  LauPdfList* sigPdfs(0);
1555  LauPdfList* scfPdfs(0);
1556 
1557  Bool_t doReweighting(kFALSE);
1558 
1559  if (curEvtCharge_<0) {
1560  model = negSigModel_;
1561  kinematics = negKinematics_;
1562  sigPdfs = &negSignalPdfs_;
1563  scfPdfs = &negScfPdfs_;
1564  if (this->enableEmbedding()) {
1565  embeddedData = negSignalTree_;
1566  doReweighting = useNegReweighting_;
1567  }
1568  } else {
1569  model = posSigModel_;
1570  kinematics = posKinematics_;
1571  if ( tagged_ ) {
1572  sigPdfs = &posSignalPdfs_;
1573  scfPdfs = &posScfPdfs_;
1574  } else {
1575  sigPdfs = &negSignalPdfs_;
1576  scfPdfs = &negScfPdfs_;
1577  }
1578  if (this->enableEmbedding()) {
1579  embeddedData = posSignalTree_;
1580  doReweighting = usePosReweighting_;
1581  }
1582  }
1583 
1584  if (this->useDP()) {
1585  if (embeddedData) {
1586  if (doReweighting) {
1587  // Select a (random) event from the generated data. Then store the
1588  // reconstructed DP co-ords, together with other pdf information,
1589  // as the embedded data.
1590  genOK = embeddedData->getReweightedEvent(model);
1591  } else {
1592  // Just get the information of a (randomly) selected event in the
1593  // embedded data
1594  embeddedData->getEmbeddedEvent(kinematics);
1595  }
1596  genSCF = this->storeSignalMCMatch( embeddedData );
1597  } else {
1598  genOK = model->generate();
1599  if ( genOK && useSCF_ ) {
1600  Double_t frac(0.0);
1601  if ( useSCFHist_ ) {
1602  frac = scfFracHist_->calcEfficiency( kinematics );
1603  } else {
1604  frac = scfFrac_.genValue();
1605  }
1606  if ( frac < LauRandom::randomFun()->Rndm() ) {
1607  this->setGenNtupleIntegerBranchValue("genTMSig",1);
1608  this->setGenNtupleIntegerBranchValue("genSCFSig",0);
1609  genSCF = kFALSE;
1610  } else {
1611  this->setGenNtupleIntegerBranchValue("genTMSig",0);
1612  this->setGenNtupleIntegerBranchValue("genSCFSig",1);
1613  genSCF = kTRUE;
1614 
1615  // Optionally smear the DP position
1616  // of the SCF event
1617  if ( scfMap_ != 0 ) {
1618  Double_t xCoord(0.0), yCoord(0.0);
1619  if ( kinematics->squareDP() ) {
1620  xCoord = kinematics->getmPrime();
1621  yCoord = kinematics->getThetaPrime();
1622  } else {
1623  xCoord = kinematics->getm13Sq();
1624  yCoord = kinematics->getm23Sq();
1625  }
1626 
1627  // Find the bin number where this event is generated
1628  Int_t binNo = scfMap_->binNumber( xCoord, yCoord );
1629 
1630  // Retrieve the migration histogram
1631  TH2* histo = scfMap_->trueHist( binNo );
1632 
1633  LauAbsEffModel * effModel = model->getEffModel();
1634  do {
1635  // Get a random point from the histogram
1636  histo->GetRandom2( xCoord, yCoord );
1637 
1638  // Update the kinematics
1639  if ( kinematics->squareDP() ) {
1640  kinematics->updateSqDPKinematics( xCoord, yCoord );
1641  } else {
1642  kinematics->updateKinematics( xCoord, yCoord );
1643  }
1644  } while ( ! effModel->passVeto( kinematics ) );
1645  }
1646  }
1647  }
1648  }
1649  } else {
1650  if (embeddedData) {
1651  embeddedData->getEmbeddedEvent(0);
1652  genSCF = this->storeSignalMCMatch( embeddedData );
1653  } else if ( useSCF_ ) {
1654  Double_t frac = scfFrac_.genValue();
1655  if ( frac < LauRandom::randomFun()->Rndm() ) {
1656  this->setGenNtupleIntegerBranchValue("genTMSig",1);
1657  this->setGenNtupleIntegerBranchValue("genSCFSig",0);
1658  genSCF = kFALSE;
1659  } else {
1660  this->setGenNtupleIntegerBranchValue("genTMSig",0);
1661  this->setGenNtupleIntegerBranchValue("genSCFSig",1);
1662  genSCF = kTRUE;
1663  }
1664  }
1665  }
1666  if (genOK) {
1667  if ( useSCF_ ) {
1668  if ( genSCF ) {
1669  this->generateExtraPdfValues(scfPdfs, embeddedData);
1670  } else {
1671  this->generateExtraPdfValues(sigPdfs, embeddedData);
1672  }
1673  } else {
1674  this->generateExtraPdfValues(sigPdfs, embeddedData);
1675  }
1676  }
1677  // Check for problems with the embedding
1678  if (embeddedData && (embeddedData->nEvents() == embeddedData->nUsedEvents())) {
1679  std::cerr << "WARNING in LauCPFitModel::generateSignalEvent : Source of embedded signal events used up, clearing the list of used events." << std::endl;
1680  embeddedData->clearUsedList();
1681  }
1682 
1683  return genOK;
1684 }
1685 
1687 {
1688  // Generate Bkgnd event
1689  Bool_t genOK(kTRUE);
1690 
1691  LauAbsBkgndDPModel* model(0);
1692  LauEmbeddedData* embeddedData(0);
1693  LauPdfList* extraPdfs(0);
1694  LauKinematics* kinematics(0);
1695 
1696  if (curEvtCharge_<0) {
1697  model = negBkgndDPModels_[bkgndID];
1698  if (this->enableEmbedding()) {
1699  embeddedData = negBkgndTree_[bkgndID];
1700  }
1701  extraPdfs = &negBkgndPdfs_[bkgndID];
1702  kinematics = negKinematics_;
1703  } else {
1704  model = posBkgndDPModels_[bkgndID];
1705  if (this->enableEmbedding()) {
1706  embeddedData = posBkgndTree_[bkgndID];
1707  }
1708  if ( tagged_ ) {
1709  extraPdfs = &posBkgndPdfs_[bkgndID];
1710  } else {
1711  extraPdfs = &negBkgndPdfs_[bkgndID];
1712  }
1713  kinematics = posKinematics_;
1714  }
1715 
1716  if (this->useDP()) {
1717  if (embeddedData) {
1718  embeddedData->getEmbeddedEvent(kinematics);
1719  } else {
1720  if (model == 0) {
1721  const TString& bkgndClass = this->bkgndClassName(bkgndID);
1722  std::cerr << "ERROR in LauCPFitModel::generateBkgndEvent : Can't find the DP model for background class \"" << bkgndClass << "\"." << std::endl;
1723  gSystem->Exit(EXIT_FAILURE);
1724  }
1725  genOK = model->generate();
1726  }
1727  } else {
1728  if (embeddedData) {
1729  embeddedData->getEmbeddedEvent(0);
1730  }
1731  }
1732  if (genOK) {
1733  this->generateExtraPdfValues(extraPdfs, embeddedData);
1734  }
1735 
1736  // Check for problems with the embedding
1737  if (embeddedData && (embeddedData->nEvents() == embeddedData->nUsedEvents())) {
1738  const TString& bkgndClass = this->bkgndClassName(bkgndID);
1739  std::cerr << "WARNING in LauCPFitModel::generateBkgndEvent : Source of embedded " << bkgndClass << " events used up, clearing the list of used events." << std::endl;
1740  embeddedData->clearUsedList();
1741  }
1742 
1743  return genOK;
1744 }
1745 
1747 {
1748  // Setup the required ntuple branches
1749  this->addGenNtupleDoubleBranch("evtWeight");
1750  this->addGenNtupleIntegerBranch("genSig");
1751  this->addGenNtupleDoubleBranch("efficiency");
1752  if ( useSCF_ || ( this->enableEmbedding() &&
1753  negSignalTree_ != 0 && negSignalTree_->haveBranch("mcMatch") &&
1754  posSignalTree_ != 0 && posSignalTree_->haveBranch("mcMatch") ) ) {
1755  this->addGenNtupleIntegerBranch("genTMSig");
1756  this->addGenNtupleIntegerBranch("genSCFSig");
1757  }
1758  const UInt_t nBkgnds = this->nBkgndClasses();
1759  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
1760  TString name( this->bkgndClassName(iBkgnd) );
1761  name.Prepend("gen");
1762  this->addGenNtupleIntegerBranch(name);
1763  }
1764  this->addGenNtupleIntegerBranch("charge");
1765  if (this->useDP() == kTRUE) {
1766  this->addGenNtupleDoubleBranch("m12");
1767  this->addGenNtupleDoubleBranch("m23");
1768  this->addGenNtupleDoubleBranch("m13");
1769  this->addGenNtupleDoubleBranch("m12Sq");
1770  this->addGenNtupleDoubleBranch("m23Sq");
1771  this->addGenNtupleDoubleBranch("m13Sq");
1772  this->addGenNtupleDoubleBranch("cosHel12");
1773  this->addGenNtupleDoubleBranch("cosHel23");
1774  this->addGenNtupleDoubleBranch("cosHel13");
1776  this->addGenNtupleDoubleBranch("mPrime");
1777  this->addGenNtupleDoubleBranch("thPrime");
1778  }
1779  }
1780  for (LauPdfList::const_iterator pdf_iter = negSignalPdfs_.begin(); pdf_iter != negSignalPdfs_.end(); ++pdf_iter) {
1781  std::vector<TString> varNames = (*pdf_iter)->varNames();
1782  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
1783  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
1784  this->addGenNtupleDoubleBranch( (*var_iter) );
1785  }
1786  }
1787  }
1788 }
1789 
1791 {
1792  LauKinematics* kinematics(0);
1793  if (curEvtCharge_<0) {
1794  kinematics = negKinematics_;
1795  } else {
1796  kinematics = posKinematics_;
1797  }
1798 
1799  // Store all the DP information
1800  this->setGenNtupleDoubleBranchValue("m12", kinematics->getm12());
1801  this->setGenNtupleDoubleBranchValue("m23", kinematics->getm23());
1802  this->setGenNtupleDoubleBranchValue("m13", kinematics->getm13());
1803  this->setGenNtupleDoubleBranchValue("m12Sq", kinematics->getm12Sq());
1804  this->setGenNtupleDoubleBranchValue("m23Sq", kinematics->getm23Sq());
1805  this->setGenNtupleDoubleBranchValue("m13Sq", kinematics->getm13Sq());
1806  this->setGenNtupleDoubleBranchValue("cosHel12", kinematics->getc12());
1807  this->setGenNtupleDoubleBranchValue("cosHel23", kinematics->getc23());
1808  this->setGenNtupleDoubleBranchValue("cosHel13", kinematics->getc13());
1809  if (kinematics->squareDP()) {
1810  this->setGenNtupleDoubleBranchValue("mPrime", kinematics->getmPrime());
1811  this->setGenNtupleDoubleBranchValue("thPrime", kinematics->getThetaPrime());
1812  }
1813 }
1814 
1816 {
1817  LauKinematics* kinematics(0);
1818  if (curEvtCharge_<0) {
1819  kinematics = negKinematics_;
1820  } else {
1821  kinematics = posKinematics_;
1822  }
1823 
1824  if (!extraPdfs) {
1825  std::cerr << "ERROR in LauCPFitModel::generateExtraPdfValues : Null pointer to PDF list." << std::endl;
1826  gSystem->Exit(EXIT_FAILURE);
1827  }
1828 
1829  if (extraPdfs->empty()) {
1830  //std::cerr << "WARNING in LauCPFitModel::generateExtraPdfValues : PDF list is empty." << std::endl;
1831  return;
1832  }
1833 
1834  // Generate from the extra PDFs
1835  for (LauPdfList::iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end(); ++pdf_iter) {
1836  LauFitData genValues;
1837  if (embeddedData) {
1838  genValues = embeddedData->getValues( (*pdf_iter)->varNames() );
1839  } else {
1840  genValues = (*pdf_iter)->generate(kinematics);
1841  }
1842  for ( LauFitData::const_iterator var_iter = genValues.begin(); var_iter != genValues.end(); ++var_iter ) {
1843  TString varName = var_iter->first;
1844  if ( varName != "m13Sq" && varName != "m23Sq" ) {
1845  Double_t value = var_iter->second;
1846  this->setGenNtupleDoubleBranchValue(varName,value);
1847  }
1848  }
1849  }
1850 }
1851 
1853 {
1854  // Default to TM
1855  Bool_t genSCF(kFALSE);
1856  Int_t match(1);
1857 
1858  // Check that we have a valid pointer and that embedded data has
1859  // the mcMatch branch. If so then get the match value.
1860  if ( embeddedData && embeddedData->haveBranch("mcMatch") ) {
1861  match = TMath::Nint( embeddedData->getValue("mcMatch") );
1862  }
1863 
1864  // Set the variables accordingly.
1865  if (match) {
1866  this->setGenNtupleIntegerBranchValue("genTMSig",1);
1867  this->setGenNtupleIntegerBranchValue("genSCFSig",0);
1868  genSCF = kFALSE;
1869  } else {
1870  this->setGenNtupleIntegerBranchValue("genTMSig",0);
1871  this->setGenNtupleIntegerBranchValue("genSCFSig",1);
1872  genSCF = kTRUE;
1873  }
1874 
1875  return genSCF;
1876 }
1877 
1879 {
1880  // Update the signal parameters and then the total normalisation for the signal likelihood
1881  if (this->useDP() == kTRUE) {
1882  this->updateCoeffs();
1885  }
1886 
1887  // Update the signal fraction from the background fractions if not doing an extended fit
1888  if ( !this->doEMLFit() && !signalEvents_->fixed() ) {
1889  this->updateSigEvents();
1890  }
1891 }
1892 
1894 {
1895  // The background parameters will have been set from Minuit.
1896  // We need to update the signal events using these.
1897  Double_t nTotEvts = this->eventsPerExpt();
1898 
1899  signalEvents_->range(-2.0*nTotEvts,2.0*nTotEvts);
1900  for (LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
1901  (*iter)->range(-2.0*nTotEvts,2.0*nTotEvts);
1902  }
1903 
1904  if (signalEvents_->fixed()) {
1905  return;
1906  }
1907 
1908  // Subtract background events (if any) from signal.
1909  Double_t signalEvents = nTotEvts;
1910  if (usingBkgnd_ == kTRUE) {
1911  for (LauBkgndYieldList::const_iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
1912  signalEvents -= (*iter)->value();
1913  }
1914  }
1915  signalEvents_->value(signalEvents);
1916 }
1917 
1919 {
1920  // Fill the internal data trees of the signal and background models.
1921  // Note that we store the events of both charges in both the
1922  // negative and the positive models. It's only later, at the stage
1923  // when the likelihood is being calculated, that we separate them.
1924 
1925  LauFitDataTree* inputFitData = this->fitData();
1926 
1927  // First the Dalitz plot variables (m_ij^2)
1928  if (this->useDP() == kTRUE) {
1929 
1930  // need to append SCF smearing bins before caching DP amplitudes
1931  if ( scfMap_ != 0 ) {
1932  this->appendBinCentres( inputFitData );
1933  }
1934 
1935  negSigModel_->fillDataTree(*inputFitData);
1936  posSigModel_->fillDataTree(*inputFitData);
1937 
1938  if (usingBkgnd_ == kTRUE) {
1939  for (LauBkgndDPModelList::iterator iter = negBkgndDPModels_.begin(); iter != negBkgndDPModels_.end(); ++iter) {
1940  (*iter)->fillDataTree(*inputFitData);
1941  }
1942  for (LauBkgndDPModelList::iterator iter = posBkgndDPModels_.begin(); iter != posBkgndDPModels_.end(); ++iter) {
1943  (*iter)->fillDataTree(*inputFitData);
1944  }
1945  }
1946  }
1947 
1948  // ...and then the extra PDFs
1949  this->cacheInfo(negSignalPdfs_, *inputFitData);
1950  this->cacheInfo(negScfPdfs_, *inputFitData);
1951  for (LauBkgndPdfsList::iterator iter = negBkgndPdfs_.begin(); iter != negBkgndPdfs_.end(); ++iter) {
1952  this->cacheInfo((*iter), *inputFitData);
1953  }
1954  if ( tagged_ ) {
1955  this->cacheInfo(posSignalPdfs_, *inputFitData);
1956  this->cacheInfo(posScfPdfs_, *inputFitData);
1957  for (LauBkgndPdfsList::iterator iter = posBkgndPdfs_.begin(); iter != posBkgndPdfs_.end(); ++iter) {
1958  this->cacheInfo((*iter), *inputFitData);
1959  }
1960  }
1961 
1962  // the SCF fractions and jacobians
1963  if ( useSCF_ && useSCFHist_ ) {
1964  if ( !inputFitData->haveBranch( "m13Sq" ) || !inputFitData->haveBranch( "m23Sq" ) ) {
1965  std::cerr << "ERROR in LauCPFitModel::cacheInputFitVars : Input data does not contain DP branches and so can't cache the SCF fraction." << std::endl;
1966  gSystem->Exit(EXIT_FAILURE);
1967  }
1968 
1969  UInt_t nEvents = inputFitData->nEvents();
1970  recoSCFFracs_.clear();
1971  recoSCFFracs_.reserve( nEvents );
1972  if ( negKinematics_->squareDP() ) {
1973  recoJacobians_.clear();
1974  recoJacobians_.reserve( nEvents );
1975  }
1976  for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) {
1977  const LauFitData& dataValues = inputFitData->getData(iEvt);
1978  LauFitData::const_iterator m13_iter = dataValues.find("m13Sq");
1979  LauFitData::const_iterator m23_iter = dataValues.find("m23Sq");
1980  negKinematics_->updateKinematics( m13_iter->second, m23_iter->second );
1981  Double_t scfFrac = scfFracHist_->calcEfficiency( negKinematics_ );
1982  recoSCFFracs_.push_back( scfFrac );
1983  if ( negKinematics_->squareDP() ) {
1985  }
1986  }
1987  }
1988 
1989  // finally cache the event charge
1990  evtCharges_.clear();
1991  if ( tagged_ ) {
1992  if ( !inputFitData->haveBranch( tagVarName_ ) ) {
1993  std::cerr << "ERROR in LauCPFitModel::cacheInputFitVars : Input data does not contain branch \"" << tagVarName_ << "\"." << std::endl;
1994  gSystem->Exit(EXIT_FAILURE);
1995  }
1996  UInt_t nEvents = inputFitData->nEvents();
1997  evtCharges_.reserve( nEvents );
1998  for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) {
1999  const LauFitData& dataValues = inputFitData->getData(iEvt);
2000  LauFitData::const_iterator iter = dataValues.find( tagVarName_ );
2001  curEvtCharge_ = static_cast<Int_t>( iter->second );
2002  evtCharges_.push_back( curEvtCharge_ );
2003  }
2004  }
2005 }
2006 
2008 {
2009  // We'll be caching the DP amplitudes and efficiencies of the centres of the true bins.
2010  // To do so, we attach some fake points at the end of inputData, the number of the entry
2011  // minus the total number of events corresponding to the number of the histogram for that
2012  // given true bin in the LauScfMap object. (What this means is that when Laura is provided with
2013  // the LauScfMap object by the user, it's the latter who has to make sure that it contains the
2014  // right number of histograms and in exactly the right order!)
2015 
2016  // Get the x and y co-ordinates of the bin centres
2017  std::vector<Double_t> binCentresXCoords;
2018  std::vector<Double_t> binCentresYCoords;
2019  scfMap_->listBinCentres(binCentresXCoords, binCentresYCoords);
2020 
2021  // The SCF histograms could be in square Dalitz plot histograms.
2022  // The dynamics takes normal Dalitz plot coords, so we might have to convert them back.
2023  Bool_t sqDP = negKinematics_->squareDP();
2024  UInt_t nBins = binCentresXCoords.size();
2025  fakeSCFFracs_.clear();
2026  fakeSCFFracs_.reserve( nBins );
2027  if ( sqDP ) {
2028  fakeJacobians_.clear();
2029  fakeJacobians_.reserve( nBins );
2030  }
2031 
2032  for (UInt_t iBin = 0; iBin < nBins; ++iBin) {
2033 
2034  if ( sqDP ) {
2035  negKinematics_->updateSqDPKinematics(binCentresXCoords[iBin],binCentresYCoords[iBin]);
2036  binCentresXCoords[iBin] = negKinematics_->getm13Sq();
2037  binCentresYCoords[iBin] = negKinematics_->getm23Sq();
2039  } else {
2040  negKinematics_->updateKinematics(binCentresXCoords[iBin],binCentresYCoords[iBin]);
2041  }
2042 
2044  }
2045 
2046  // Set up inputFitVars_ object to hold the fake events
2047  inputData->appendFakePoints(binCentresXCoords,binCentresYCoords);
2048 }
2049 
2051 {
2052  // Find out whether we have B- or B+
2053  if ( tagged_ ) {
2054  curEvtCharge_ = evtCharges_[iEvt];
2055 
2056  // check that the charge is either +1 or -1
2057  if (TMath::Abs(curEvtCharge_)!=1) {
2058  std::cerr << "ERROR in LauCPFitModel::getTotEvtLikelihood : Charge/tag not accepted value: " << curEvtCharge_ << std::endl;
2059  if (curEvtCharge_>0) {
2060  curEvtCharge_ = +1;
2061  } else {
2062  curEvtCharge_ = -1;
2063  }
2064  std::cerr << " : Making it: " << curEvtCharge_ << "." << std::endl;
2065  }
2066  }
2067 
2068  // Get the DP likelihood for signal and backgrounds
2069  this->getEvtDPLikelihood(iEvt);
2070 
2071  // Get the combined extra PDFs likelihood for signal and backgrounds
2072  this->getEvtExtraLikelihoods(iEvt);
2073 
2074  // If appropriate, combine the TM and SCF likelihoods
2075  Double_t sigLike = sigDPLike_ * sigExtraLike_;
2076  if ( useSCF_ ) {
2077  Double_t scfFrac(0.0);
2078  if (useSCFHist_) {
2079  scfFrac = recoSCFFracs_[iEvt];
2080  } else {
2081  scfFrac = scfFrac_.value();
2082  }
2083  sigLike *= (1.0 - scfFrac);
2084  if ( (scfMap_ != 0) && (this->useDP() == kTRUE) ) {
2085  // if we're smearing the SCF DP PDF then the SCF frac
2086  // is already included in the SCF DP likelihood
2087  sigLike += (scfDPLike_ * scfExtraLike_);
2088  } else {
2089  sigLike += (scfFrac * scfDPLike_ * scfExtraLike_);
2090  }
2091  }
2092 
2093  // Get the correct event fractions depending on the charge
2094  // Signal asymmetry is built into the DP model... but when the DP
2095  // isn't in the fit we need an explicit parameter
2096  Double_t signalEvents = signalEvents_->value() * 0.5;
2097  if (this->useDP() == kFALSE) {
2098  signalEvents *= (1.0 - curEvtCharge_ * signalAsym_->value());
2099  }
2100 
2101  // Construct the total event likelihood
2102  Double_t likelihood(0.0);
2103  if (usingBkgnd_) {
2104  likelihood = sigLike*signalEvents;
2105  const UInt_t nBkgnds = this->nBkgndClasses();
2106  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
2107  Double_t bkgndEvents = bkgndEvents_[bkgndID]->value() * 0.5 * (1.0 - curEvtCharge_ * bkgndAsym_[bkgndID]->value());
2108  likelihood += bkgndEvents*bkgndDPLike_[bkgndID]*bkgndExtraLike_[bkgndID];
2109  }
2110  } else {
2111  likelihood = sigLike*0.5;
2112  }
2113  return likelihood;
2114 }
2115 
2117 {
2118  Double_t eventSum(0.0);
2119  eventSum += signalEvents_->value();
2120  if (usingBkgnd_) {
2121  for (LauBkgndYieldList::const_iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
2122  eventSum += (*iter)->value();
2123  }
2124  }
2125  return eventSum;
2126 }
2127 
2129 {
2130  // Function to return the signal and background likelihoods for the
2131  // Dalitz plot for the given event evtNo.
2132 
2133  if ( ! this->useDP() ) {
2134  // There's always going to be a term in the likelihood for the
2135  // signal, so we'd better not zero it.
2136  sigDPLike_ = 1.0;
2137  scfDPLike_ = 1.0;
2138 
2139  const UInt_t nBkgnds = this->nBkgndClasses();
2140  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
2141  if (usingBkgnd_ == kTRUE) {
2142  bkgndDPLike_[bkgndID] = 1.0;
2143  } else {
2144  bkgndDPLike_[bkgndID] = 0.0;
2145  }
2146  }
2147 
2148  return;
2149  }
2150 
2151  const UInt_t nBkgnds = this->nBkgndClasses();
2152  if ( tagged_ ) {
2153  if (curEvtCharge_==+1) {
2157 
2158  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
2159  if (usingBkgnd_ == kTRUE) {
2160  bkgndDPLike_[bkgndID] = posBkgndDPModels_[bkgndID]->getLikelihood(iEvt);
2161  } else {
2162  bkgndDPLike_[bkgndID] = 0.0;
2163  }
2164  }
2165  } else {
2169 
2170  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
2171  if (usingBkgnd_ == kTRUE) {
2172  bkgndDPLike_[bkgndID] = negBkgndDPModels_[bkgndID]->getLikelihood(iEvt);
2173  } else {
2174  bkgndDPLike_[bkgndID] = 0.0;
2175  }
2176  }
2177  }
2178  } else {
2181 
2184 
2185  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
2186  if (usingBkgnd_ == kTRUE) {
2187  bkgndDPLike_[bkgndID] = 0.5 * ( posBkgndDPModels_[bkgndID]->getLikelihood(iEvt) +
2188  negBkgndDPModels_[bkgndID]->getLikelihood(iEvt) );
2189  } else {
2190  bkgndDPLike_[bkgndID] = 0.0;
2191  }
2192  }
2193  }
2194 
2195  if ( useSCF_ == kTRUE ) {
2196  if ( scfMap_ == 0 ) {
2197  // we're not smearing the SCF DP position
2198  // so the likelihood is the same as the TM
2200  } else {
2201  // calculate the smeared SCF DP likelihood
2202  scfDPLike_ = this->getEvtSCFDPLikelihood(iEvt);
2203  }
2204  }
2205 
2206  // Calculate the signal normalisation
2207  // NB the 2.0 is there so that the 0.5 factor is applied to
2208  // signal and background in the same place otherwise you get
2209  // normalisation problems when you switch off the DP in the fit
2210  Double_t norm = negSigModel_->getDPNorm() + posSigModel_->getDPNorm();
2211  sigDPLike_ *= 2.0/norm;
2212  scfDPLike_ *= 2.0/norm;
2213 }
2214 
2216 {
2217  Double_t scfDPLike(0.0);
2218 
2219  Double_t recoJacobian(1.0);
2220  Double_t xCoord(0.0);
2221  Double_t yCoord(0.0);
2222 
2223  Bool_t squareDP = negKinematics_->squareDP();
2224  if ( squareDP ) {
2225  xCoord = negSigModel_->getEvtmPrime();
2226  yCoord = negSigModel_->getEvtthPrime();
2227  recoJacobian = recoJacobians_[iEvt];
2228  } else {
2229  xCoord = negSigModel_->getEvtm13Sq();
2230  yCoord = negSigModel_->getEvtm23Sq();
2231  }
2232 
2233  // Find the bin that our reco event falls in
2234  Int_t recoBin = scfMap_->binNumber( xCoord, yCoord );
2235 
2236  // Find out which true Bins contribute to the given reco bin
2237  const std::vector<Int_t>* trueBins = scfMap_->trueBins(recoBin);
2238 
2239  const Int_t nDataEvents = this->eventsPerExpt();
2240 
2241  // Loop over the true bins
2242  for (std::vector<Int_t>::const_iterator iter = trueBins->begin(); iter != trueBins->end(); ++iter)
2243  {
2244  Int_t trueBin = (*iter);
2245 
2246  // prob of a true event in the given true bin migrating to the reco bin
2247  Double_t pRecoGivenTrue = scfMap_->prob( recoBin, trueBin );
2248  Double_t pTrue(0.0);
2249 
2250  // We've cached the DP amplitudes and the efficiency for the
2251  // true bin centres, just after the data points
2252  if ( tagged_ ) {
2253  LauIsobarDynamics* sigModel(0);
2254  if (curEvtCharge_<0) {
2255  sigModel = negSigModel_;
2256  } else {
2257  sigModel = posSigModel_;
2258  }
2259 
2260  sigModel->calcLikelihoodInfo( nDataEvents + trueBin );
2261 
2262  pTrue = sigModel->getEvtDPAmp().abs2() * sigModel->getEvtEff();
2263  } else {
2264  posSigModel_->calcLikelihoodInfo( nDataEvents + trueBin );
2265  negSigModel_->calcLikelihoodInfo( nDataEvents + trueBin );
2266 
2267  pTrue = 0.5 * ( posSigModel_->getEvtDPAmp().abs2() * posSigModel_->getEvtEff() +
2269  }
2270 
2271  // Get the cached SCF fraction (and jacobian if we're using the square DP)
2272  Double_t scfFraction = fakeSCFFracs_[ trueBin ];
2273  Double_t jacobian(1.0);
2274  if ( squareDP ) {
2275  jacobian = fakeJacobians_[ trueBin ];
2276  }
2277 
2278  scfDPLike += pTrue * jacobian * scfFraction * pRecoGivenTrue;
2279  }
2280 
2281  // Divide by the reco jacobian
2282  scfDPLike /= recoJacobian;
2283 
2284  return scfDPLike;
2285 }
2286 
2288 {
2289  // Function to return the signal and background likelihoods for the
2290  // extra variables for the given event evtNo.
2291 
2292  sigExtraLike_ = 1.0;
2293 
2294  const UInt_t nBkgnds = this->nBkgndClasses();
2295 
2296  if ( ! tagged_ || curEvtCharge_ < 0 ) {
2297  sigExtraLike_ = this->prodPdfValue( negSignalPdfs_, iEvt );
2298 
2299  if (useSCF_) {
2300  scfExtraLike_ = this->prodPdfValue( negScfPdfs_, iEvt );
2301  }
2302 
2303  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
2304  if (usingBkgnd_) {
2305  bkgndExtraLike_[bkgndID] = this->prodPdfValue( negBkgndPdfs_[bkgndID], iEvt );
2306  } else {
2307  bkgndExtraLike_[bkgndID] = 0.0;
2308  }
2309  }
2310  } else {
2311  sigExtraLike_ = this->prodPdfValue( posSignalPdfs_, iEvt );
2312 
2313  if (useSCF_) {
2314  scfExtraLike_ = this->prodPdfValue( posScfPdfs_, iEvt );
2315  }
2316 
2317  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
2318  if (usingBkgnd_) {
2319  bkgndExtraLike_[bkgndID] = this->prodPdfValue( posBkgndPdfs_[bkgndID], iEvt );
2320  } else {
2321  bkgndExtraLike_[bkgndID] = 0.0;
2322  }
2323  }
2324  }
2325 }
2326 
2328 {
2329  negCoeffs_.clear(); posCoeffs_.clear();
2330  negCoeffs_.reserve(nSigComp_); posCoeffs_.reserve(nSigComp_);
2331  for (UInt_t i = 0; i < nSigComp_; i++) {
2332  negCoeffs_.push_back(coeffPars_[i]->antiparticleCoeff());
2333  posCoeffs_.push_back(coeffPars_[i]->particleCoeff());
2334  }
2335 }
2336 
2338 {
2339  // add branches for storing the experiment number and the number of
2340  // the event within the current experiment
2341  this->addSPlotNtupleIntegerBranch("iExpt");
2342  this->addSPlotNtupleIntegerBranch("iEvtWithinExpt");
2343 
2344  // Store the efficiency of the event (for inclusive BF calculations).
2345  if (this->storeDPEff()) {
2346  this->addSPlotNtupleDoubleBranch("efficiency");
2348  this->addSPlotNtupleDoubleBranch("scffraction");
2349  }
2350  }
2351 
2352  // Store the total event likelihood for each species.
2353  if (useSCF_) {
2354  this->addSPlotNtupleDoubleBranch("sigTMTotalLike");
2355  this->addSPlotNtupleDoubleBranch("sigSCFTotalLike");
2356  this->addSPlotNtupleDoubleBranch("sigSCFFrac");
2357  } else {
2358  this->addSPlotNtupleDoubleBranch("sigTotalLike");
2359  }
2360  if (usingBkgnd_) {
2361  const UInt_t nBkgnds = this->nBkgndClasses();
2362  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
2363  TString name( this->bkgndClassName(iBkgnd) );
2364  name += "TotalLike";
2365  this->addSPlotNtupleDoubleBranch(name);
2366  }
2367  }
2368 
2369  // Store the DP likelihoods
2370  if (this->useDP()) {
2371  if (useSCF_) {
2372  this->addSPlotNtupleDoubleBranch("sigTMDPLike");
2373  this->addSPlotNtupleDoubleBranch("sigSCFDPLike");
2374  } else {
2375  this->addSPlotNtupleDoubleBranch("sigDPLike");
2376  }
2377  if (usingBkgnd_) {
2378  const UInt_t nBkgnds = this->nBkgndClasses();
2379  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
2380  TString name( this->bkgndClassName(iBkgnd) );
2381  name += "DPLike";
2382  this->addSPlotNtupleDoubleBranch(name);
2383  }
2384  }
2385  }
2386 
2387  // Store the likelihoods for each extra PDF
2388  if (useSCF_) {
2389  this->addSPlotNtupleBranches(&negSignalPdfs_, "sigTM");
2390  this->addSPlotNtupleBranches(&negScfPdfs_, "sigSCF");
2391  } else {
2392  this->addSPlotNtupleBranches(&negSignalPdfs_, "sig");
2393  }
2394  if (usingBkgnd_) {
2395  const UInt_t nBkgnds = this->nBkgndClasses();
2396  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
2397  const TString& bkgndClass = this->bkgndClassName(iBkgnd);
2398  const LauPdfList* pdfList = &(negBkgndPdfs_[iBkgnd]);
2399  this->addSPlotNtupleBranches(pdfList, bkgndClass);
2400  }
2401  }
2402 }
2403 
2404 void LauCPFitModel::addSPlotNtupleBranches(const LauPdfList* extraPdfs, const TString& prefix)
2405 {
2406  if (extraPdfs) {
2407  // Loop through each of the PDFs
2408  for (LauPdfList::const_iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end(); ++pdf_iter) {
2409 
2410  // Count the number of input variables that are not
2411  // DP variables (used in the case where there is DP
2412  // dependence for e.g. MVA)
2413  UInt_t nVars(0);
2414  std::vector<TString> varNames = (*pdf_iter)->varNames();
2415  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
2416  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
2417  ++nVars;
2418  }
2419  }
2420 
2421  if ( nVars == 1 ) {
2422  // If the PDF only has one variable then
2423  // simply add one branch for that variable
2424  TString varName = (*pdf_iter)->varName();
2425  TString name(prefix);
2426  name += varName;
2427  name += "Like";
2428  this->addSPlotNtupleDoubleBranch(name);
2429  } else if ( nVars == 2 ) {
2430  // If the PDF has two variables then we
2431  // need a branch for them both together and
2432  // branches for each
2433  TString allVars("");
2434  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
2435  allVars += (*var_iter);
2436  TString name(prefix);
2437  name += (*var_iter);
2438  name += "Like";
2439  this->addSPlotNtupleDoubleBranch(name);
2440  }
2441  TString name(prefix);
2442  name += allVars;
2443  name += "Like";
2444  this->addSPlotNtupleDoubleBranch(name);
2445  } else {
2446  std::cerr << "WARNING in LauCPFitModel::addSPlotNtupleBranches : Can't yet deal with 3D PDFs." << std::endl;
2447  }
2448  }
2449  }
2450 }
2451 
2452 Double_t LauCPFitModel::setSPlotNtupleBranchValues(LauPdfList* extraPdfs, const TString& prefix, UInt_t iEvt)
2453 {
2454  // Store the PDF value for each variable in the list
2455  Double_t totalLike(1.0);
2456  Double_t extraLike(0.0);
2457  if (extraPdfs) {
2458  for (LauPdfList::iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end(); ++pdf_iter) {
2459 
2460  // calculate the likelihood for this event
2461  (*pdf_iter)->calcLikelihoodInfo(iEvt);
2462  extraLike = (*pdf_iter)->getLikelihood();
2463  totalLike *= extraLike;
2464 
2465  // Count the number of input variables that are not
2466  // DP variables (used in the case where there is DP
2467  // dependence for e.g. MVA)
2468  UInt_t nVars(0);
2469  std::vector<TString> varNames = (*pdf_iter)->varNames();
2470  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
2471  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
2472  ++nVars;
2473  }
2474  }
2475 
2476  if ( nVars == 1 ) {
2477  // If the PDF only has one variable then
2478  // simply store the value for that variable
2479  TString varName = (*pdf_iter)->varName();
2480  TString name(prefix);
2481  name += varName;
2482  name += "Like";
2483  this->setSPlotNtupleDoubleBranchValue(name, extraLike);
2484  } else if ( nVars == 2 ) {
2485  // If the PDF has two variables then we
2486  // store the value for them both together
2487  // and for each on their own
2488  TString allVars("");
2489  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
2490  allVars += (*var_iter);
2491  TString name(prefix);
2492  name += (*var_iter);
2493  name += "Like";
2494  Double_t indivLike = (*pdf_iter)->getLikelihood( (*var_iter) );
2495  this->setSPlotNtupleDoubleBranchValue(name, indivLike);
2496  }
2497  TString name(prefix);
2498  name += allVars;
2499  name += "Like";
2500  this->setSPlotNtupleDoubleBranchValue(name, extraLike);
2501  } else {
2502  std::cerr << "WARNING in LauCPFitModel::setSPlotNtupleBranchValues : Can't yet deal with 3D PDFs." << std::endl;
2503  }
2504  }
2505  }
2506  return totalLike;
2507 }
2508 
2510 {
2511  LauSPlot::NameSet nameSet;
2512  if (this->useDP()) {
2513  nameSet.insert("DP");
2514  }
2515  // Loop through all the signal PDFs
2516  for (LauPdfList::const_iterator pdf_iter = negSignalPdfs_.begin(); pdf_iter != negSignalPdfs_.end(); ++pdf_iter) {
2517  // Loop over the variables involved in each PDF
2518  std::vector<TString> varNames = (*pdf_iter)->varNames();
2519  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
2520  // If they are not DP coordinates then add them
2521  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
2522  nameSet.insert( (*var_iter) );
2523  }
2524  }
2525  }
2526  return nameSet;
2527 }
2528 
2530 {
2531  LauSPlot::NumbMap numbMap;
2532  if (!signalEvents_->fixed() && this->doEMLFit()) {
2533  numbMap["sig"] = signalEvents_->genValue();
2534  }
2535  if ( usingBkgnd_ ) {
2536  const UInt_t nBkgnds = this->nBkgndClasses();
2537  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
2538  const TString& bkgndClass = this->bkgndClassName(iBkgnd);
2539  const LauParameter* par = bkgndEvents_[iBkgnd];
2540  if (!par->fixed()) {
2541  numbMap[bkgndClass] = par->genValue();
2542  }
2543  }
2544  }
2545  return numbMap;
2546 }
2547 
2549 {
2550  LauSPlot::NumbMap numbMap;
2551  if (signalEvents_->fixed() && this->doEMLFit()) {
2552  numbMap["sig"] = signalEvents_->genValue();
2553  }
2554  if ( usingBkgnd_ ) {
2555  const UInt_t nBkgnds = this->nBkgndClasses();
2556  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
2557  const TString& bkgndClass = this->bkgndClassName(iBkgnd);
2558  const LauParameter* par = bkgndEvents_[iBkgnd];
2559  if (par->fixed()) {
2560  numbMap[bkgndClass] = par->genValue();
2561  }
2562  }
2563  }
2564  return numbMap;
2565 }
2566 
2568 {
2569  // This makes the assumption that the form of the positive and
2570  // negative PDFs are the same, which seems reasonable to me
2571 
2572  LauSPlot::TwoDMap twodimMap;
2573 
2574  for (LauPdfList::const_iterator pdf_iter = negSignalPdfs_.begin(); pdf_iter != negSignalPdfs_.end(); ++pdf_iter) {
2575  // Count the number of input variables that are not DP variables
2576  UInt_t nVars(0);
2577  std::vector<TString> varNames = (*pdf_iter)->varNames();
2578  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
2579  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
2580  ++nVars;
2581  }
2582  }
2583  if ( nVars == 2 ) {
2584  if (useSCF_) {
2585  twodimMap.insert( std::make_pair( "sigTM", std::make_pair( varNames[0], varNames[1] ) ) );
2586  } else {
2587  twodimMap.insert( std::make_pair( "sig", std::make_pair( varNames[0], varNames[1] ) ) );
2588  }
2589  }
2590  }
2591 
2592  if ( useSCF_ ) {
2593  for (LauPdfList::const_iterator pdf_iter = negScfPdfs_.begin(); pdf_iter != negScfPdfs_.end(); ++pdf_iter) {
2594  // Count the number of input variables that are not DP variables
2595  UInt_t nVars(0);
2596  std::vector<TString> varNames = (*pdf_iter)->varNames();
2597  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
2598  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
2599  ++nVars;
2600  }
2601  }
2602  if ( nVars == 2 ) {
2603  twodimMap.insert( std::make_pair( "sigSCF", std::make_pair( varNames[0], varNames[1] ) ) );
2604  }
2605  }
2606  }
2607 
2608  if (usingBkgnd_) {
2609  const UInt_t nBkgnds = this->nBkgndClasses();
2610  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
2611  const TString& bkgndClass = this->bkgndClassName(iBkgnd);
2612  const LauPdfList& pdfList = negBkgndPdfs_[iBkgnd];
2613  for (LauPdfList::const_iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) {
2614  // Count the number of input variables that are not DP variables
2615  UInt_t nVars(0);
2616  std::vector<TString> varNames = (*pdf_iter)->varNames();
2617  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
2618  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
2619  ++nVars;
2620  }
2621  }
2622  if ( nVars == 2 ) {
2623  twodimMap.insert( std::make_pair( bkgndClass, std::make_pair( varNames[0], varNames[1] ) ) );
2624  }
2625  }
2626  }
2627  }
2628 
2629  return twodimMap;
2630 }
2631 
2633 {
2634  std::cout << "INFO in LauCPFitModel::storePerEvtLlhds : Storing per-event likelihood values..." << std::endl;
2635 
2636  // if we've not been using the DP model then we need to cache all
2637  // the info here so that we can get the efficiency from it
2638 
2639  LauFitDataTree* inputFitData = this->fitData();
2640 
2641  if (!this->useDP() && this->storeDPEff()) {
2644  negSigModel_->fillDataTree(*inputFitData);
2645  posSigModel_->fillDataTree(*inputFitData);
2646  }
2647 
2648  UInt_t evtsPerExpt(this->eventsPerExpt());
2649 
2650  LauIsobarDynamics* sigModel(0);
2651  LauPdfList* sigPdfs(0);
2652  LauPdfList* scfPdfs(0);
2653  LauBkgndPdfsList* bkgndPdfs(0);
2654 
2655  for (UInt_t iEvt = 0; iEvt < evtsPerExpt; ++iEvt) {
2656 
2657  this->setSPlotNtupleIntegerBranchValue("iExpt",this->iExpt());
2658  this->setSPlotNtupleIntegerBranchValue("iEvtWithinExpt",iEvt);
2659 
2660  // Find out whether we have B- or B+
2661  if ( tagged_ ) {
2662  const LauFitData& dataValues = inputFitData->getData(iEvt);
2663  LauFitData::const_iterator iter = dataValues.find("charge");
2664  curEvtCharge_ = static_cast<Int_t>(iter->second);
2665 
2666  if (curEvtCharge_==+1) {
2667  sigModel = posSigModel_;
2668  sigPdfs = &posSignalPdfs_;
2669  scfPdfs = &posScfPdfs_;
2670  bkgndPdfs = &posBkgndPdfs_;
2671  } else {
2672  sigModel = negSigModel_;
2673  sigPdfs = &negSignalPdfs_;
2674  scfPdfs = &negScfPdfs_;
2675  bkgndPdfs = &negBkgndPdfs_;
2676  }
2677  } else {
2678  sigPdfs = &negSignalPdfs_;
2679  scfPdfs = &negScfPdfs_;
2680  bkgndPdfs = &negBkgndPdfs_;
2681  }
2682 
2683  // the DP information
2684  this->getEvtDPLikelihood(iEvt);
2685  if (this->storeDPEff()) {
2686  if (!this->useDP()) {
2689  }
2690  if ( tagged_ ) {
2691  this->setSPlotNtupleDoubleBranchValue("efficiency",sigModel->getEvtEff());
2693  this->setSPlotNtupleDoubleBranchValue("scffraction",sigModel->getEvtScfFraction());
2694  }
2695  } else {
2699  }
2700  }
2701  }
2702  if (this->useDP()) {
2704  if (useSCF_) {
2706  this->setSPlotNtupleDoubleBranchValue("sigTMDPLike",sigDPLike_);
2707  this->setSPlotNtupleDoubleBranchValue("sigSCFDPLike",scfDPLike_);
2708  } else {
2709  this->setSPlotNtupleDoubleBranchValue("sigDPLike",sigDPLike_);
2710  }
2711  if (usingBkgnd_) {
2712  const UInt_t nBkgnds = this->nBkgndClasses();
2713  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
2714  TString name = this->bkgndClassName(iBkgnd);
2715  name += "DPLike";
2716  this->setSPlotNtupleDoubleBranchValue(name,bkgndDPLike_[iBkgnd]);
2717  }
2718  }
2719  } else {
2720  sigTotalLike_ = 1.0;
2721  if (useSCF_) {
2722  scfTotalLike_ = 1.0;
2723  }
2724  if (usingBkgnd_) {
2725  const UInt_t nBkgnds = this->nBkgndClasses();
2726  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
2727  bkgndTotalLike_[iBkgnd] = 1.0;
2728  }
2729  }
2730  }
2731 
2732  // the signal PDF values
2733  if ( useSCF_ ) {
2734  sigTotalLike_ *= this->setSPlotNtupleBranchValues(sigPdfs, "sigTM", iEvt);
2735  scfTotalLike_ *= this->setSPlotNtupleBranchValues(scfPdfs, "sigSCF", iEvt);
2736  } else {
2737  sigTotalLike_ *= this->setSPlotNtupleBranchValues(sigPdfs, "sig", iEvt);
2738  }
2739 
2740  // the background PDF values
2741  if (usingBkgnd_) {
2742  const UInt_t nBkgnds = this->nBkgndClasses();
2743  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
2744  const TString& bkgndClass = this->bkgndClassName(iBkgnd);
2745  LauPdfList& pdfs = (*bkgndPdfs)[iBkgnd];
2746  bkgndTotalLike_[iBkgnd] *= this->setSPlotNtupleBranchValues(&(pdfs), bkgndClass, iEvt);
2747  }
2748  }
2749 
2750  // the total likelihoods
2751  if (useSCF_) {
2752  Double_t scfFrac(0.0);
2753  if ( useSCFHist_ ) {
2754  scfFrac = recoSCFFracs_[iEvt];
2755  } else {
2756  scfFrac = scfFrac_.value();
2757  }
2758  this->setSPlotNtupleDoubleBranchValue("sigSCFFrac",scfFrac);
2759  sigTotalLike_ *= ( 1.0 - scfFrac );
2760  if ( scfMap_ == 0 ) {
2761  scfTotalLike_ *= scfFrac;
2762  }
2763  this->setSPlotNtupleDoubleBranchValue("sigTMTotalLike",sigTotalLike_);
2764  this->setSPlotNtupleDoubleBranchValue("sigSCFTotalLike",scfTotalLike_);
2765  } else {
2766  this->setSPlotNtupleDoubleBranchValue("sigTotalLike",sigTotalLike_);
2767  }
2768  if (usingBkgnd_) {
2769  const UInt_t nBkgnds = this->nBkgndClasses();
2770  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
2771  TString name = this->bkgndClassName(iBkgnd);
2772  name += "TotalLike";
2774  }
2775  }
2776  // fill the tree
2777  this->fillSPlotNtupleBranches();
2778  }
2779  std::cout << "INFO in LauCPFitModel::storePerEvtLlhds : Finished storing per-event likelihood values." << std::endl;
2780 }
2781 
2782 void LauCPFitModel::embedNegSignal(const TString& fileName, const TString& treeName,
2783  Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment,
2784  Bool_t useReweighting)
2785 {
2786  if (negSignalTree_) {
2787  std::cerr << "ERROR in LauCPFitModel::embedNegSignal : Already embedding signal from a file." << std::endl;
2788  return;
2789  }
2790 
2791  if (!reuseEventsWithinEnsemble && reuseEventsWithinExperiment) {
2792  std::cerr << "WARNING in LauCPFitModel::embedNegSignal : Conflicting options provided, will not reuse events at all." << std::endl;
2793  reuseEventsWithinExperiment = kFALSE;
2794  }
2795 
2796  negSignalTree_ = new LauEmbeddedData(fileName,treeName,reuseEventsWithinExperiment);
2797  Bool_t dataOK = negSignalTree_->findBranches();
2798  if (!dataOK) {
2799  delete negSignalTree_; negSignalTree_ = 0;
2800  std::cerr << "ERROR in LauCPFitModel::embedNegSignal : Problem creating data tree for embedding." << std::endl;
2801  return;
2802  }
2803  reuseSignal_ = reuseEventsWithinEnsemble;
2804  useNegReweighting_ = useReweighting;
2805  if (this->enableEmbedding() == kFALSE) {this->enableEmbedding(kTRUE);}
2806 }
2807 
2808 void LauCPFitModel::embedNegBkgnd(const TString& bkgndClass, const TString& fileName, const TString& treeName,
2809  Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment)
2810 {
2811  if ( ! this->validBkgndClass( bkgndClass ) ) {
2812  std::cerr << "ERROR in LauCPFitModel::embedBkgnd : Invalid background class \"" << bkgndClass << "\"." << std::endl;
2813  std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
2814  return;
2815  }
2816 
2817  UInt_t bkgndID = this->bkgndClassID( bkgndClass );
2818 
2819  if (negBkgndTree_[bkgndID]) {
2820  std::cerr << "ERROR in LauCPFitModel::embedNegBkgnd : Already embedding background from a file." << std::endl;
2821  return;
2822  }
2823 
2824  if (!reuseEventsWithinEnsemble && reuseEventsWithinExperiment) {
2825  std::cerr << "WARNING in LauCPFitModel::embedNegBkgnd : Conflicting options provided, will not reuse events at all." << std::endl;
2826  reuseEventsWithinExperiment = kFALSE;
2827  }
2828 
2829  negBkgndTree_[bkgndID] = new LauEmbeddedData(fileName,treeName,reuseEventsWithinExperiment);
2830  Bool_t dataOK = negBkgndTree_[bkgndID]->findBranches();
2831  if (!dataOK) {
2832  delete negBkgndTree_[bkgndID]; negBkgndTree_[bkgndID] = 0;
2833  std::cerr << "ERROR in LauCPFitModel::embedNegBkgnd : Problem creating data tree for embedding." << std::endl;
2834  return;
2835  }
2836  reuseBkgnd_[bkgndID] = reuseEventsWithinEnsemble;
2837  if (this->enableEmbedding() == kFALSE) {this->enableEmbedding(kTRUE);}
2838 }
2839 
2840 void LauCPFitModel::embedPosSignal(const TString& fileName, const TString& treeName,
2841  Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment,
2842  Bool_t useReweighting)
2843 {
2844  if (posSignalTree_) {
2845  std::cerr << "ERROR in LauCPFitModel::embedPosSignal : Already embedding signal from a file." << std::endl;
2846  return;
2847  }
2848 
2849  if (!reuseEventsWithinEnsemble && reuseEventsWithinExperiment) {
2850  std::cerr << "WARNING in LauCPFitModel::embedPosSignal : Conflicting options provided, will not reuse events at all." << std::endl;
2851  reuseEventsWithinExperiment = kFALSE;
2852  }
2853 
2854  posSignalTree_ = new LauEmbeddedData(fileName,treeName,reuseEventsWithinExperiment);
2855  Bool_t dataOK = posSignalTree_->findBranches();
2856  if (!dataOK) {
2857  delete posSignalTree_; posSignalTree_ = 0;
2858  std::cerr << "ERROR in LauCPFitModel::embedPosSignal : Problem creating data tree for embedding." << std::endl;
2859  return;
2860  }
2861  reuseSignal_ = reuseEventsWithinEnsemble;
2862  usePosReweighting_ = useReweighting;
2863  if (this->enableEmbedding() == kFALSE) {this->enableEmbedding(kTRUE);}
2864 }
2865 
2866 void LauCPFitModel::embedPosBkgnd(const TString& bkgndClass, const TString& fileName, const TString& treeName,
2867  Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment)
2868 {
2869  if ( ! this->validBkgndClass( bkgndClass ) ) {
2870  std::cerr << "ERROR in LauCPFitModel::embedBkgnd : Invalid background class \"" << bkgndClass << "\"." << std::endl;
2871  std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
2872  return;
2873  }
2874 
2875  UInt_t bkgndID = this->bkgndClassID( bkgndClass );
2876 
2877  if (posBkgndTree_[bkgndID]) {
2878  std::cerr << "ERROR in LauCPFitModel::embedPosBkgnd : Already embedding background from a file." << std::endl;
2879  return;
2880  }
2881 
2882  if (!reuseEventsWithinEnsemble && reuseEventsWithinExperiment) {
2883  std::cerr << "WARNING in LauCPFitModel::embedPosBkgnd : Conflicting options provided, will not reuse events at all." << std::endl;
2884  reuseEventsWithinExperiment = kFALSE;
2885  }
2886 
2887  posBkgndTree_[bkgndID] = new LauEmbeddedData(fileName,treeName,reuseEventsWithinExperiment);
2888  Bool_t dataOK = posBkgndTree_[bkgndID]->findBranches();
2889  if (!dataOK) {
2890  delete posBkgndTree_[bkgndID]; posBkgndTree_[bkgndID] = 0;
2891  std::cerr << "ERROR in LauCPFitModel::embedPosBkgnd : Problem creating data tree for embedding." << std::endl;
2892  return;
2893  }
2894  reuseBkgnd_[bkgndID] = reuseEventsWithinEnsemble;
2895  if (this->enableEmbedding() == kFALSE) {this->enableEmbedding(kTRUE);}
2896 }
2897 
2898 void LauCPFitModel::weightEvents( const TString& /*dataFileName*/, const TString& /*dataTreeName*/ )
2899 {
2900  std::cerr << "ERROR in LauCPFitModel::weightEvents : Method not available for this fit model." << std::endl;
2901  return;
2902 }
2903 
Double_t range() const
The range allowed for the parameter.
void generateExtraPdfValues(LauPdfList *extraPdfs, LauEmbeddedData *embeddedData)
Generate from the extra PDFs.
void modifyDataTree()
Recache the amplitude values for those that have changed.
Bool_t generateSignalEvent()
Generate signal event.
virtual void addSPlotNtupleIntegerBranch(const TString &name)
Add a branch to the sPlot tree for storing an integer.
The abstract interface for a background Dalitz plot model.
Bool_t squareDP() const
Are the square Dalitz plot co-ordinates being calculated?
std::vector< LauParameter > acp_
A_CP parameter.
virtual LauSPlot::NameSet variableNames() const
Returns the names of all variables in the fit.
Double_t getc23() const
Get the cosine of the helicity angle theta23.
TRandom * randomFun()
Access the singleton random number generator with a particular seed.
Definition: LauRandom.cc:20
std::vector< LauParameter > LauParameterList
List of parameters.
Bool_t fixed() const
Check whether the parameter is fixed or floated.
LauKinematics * negKinematics_
The B- Dalitz plot kinematics object.
UInt_t nEvents() const
Retrieve the number of events.
UInt_t eventsPerExpt() const
Obtain the total number of events in the current experiment.
virtual void printAsymmetries(std::ostream &output)
Print the asymmetries.
Bool_t writeLatexTable() const
Determine whether writing out of the latex table is enabled.
LauPdfList negSignalPdfs_
The B- signal PDFs.
void addSPlotNtupleBranches(const LauPdfList *extraPdfs, const TString &prefix)
Add sPlot branches for the extra PDFs.
Bool_t forceAsym_
Option to force an asymmetry.
Double_t scfTotalLike_
Total SCF likelihood.
void setExtraPdfParameters()
Set the fit parameters for the extra PDFs.
void updateCoeffs(const std::vector< LauComplex > &coeffs)
Update the complex coefficients for the resonances.
void updateSqDPKinematics(Double_t mPrime, Double_t thetaPrime)
Update all kinematic quantities based on the square DP co-ordinates m&#39; and theta&#39;.
LauBkgndDPModelList posBkgndDPModels_
The B+ background Dalitz plot models.
void setEffHisto(const TH2 *effHisto, Bool_t useInterpolation=kTRUE, Bool_t fluctuateBins=kFALSE, Double_t avEff=-1.0, Double_t absError=-1.0, Bool_t useUpperHalfOnly=kFALSE, Bool_t squareDP=kFALSE)
Set the efficiency variation across the phase space using a predetermined 2D histogram.
Definition: LauEffModel.cc:55
void cacheInfo(LauPdfList &pdfList, const LauFitDataTree &theData)
Have all PDFs in the list cache the data.
virtual void setGenNtupleIntegerBranchValue(const TString &name, Int_t value)
Set the value of an integer branch in the gen tree.
File containing declaration of LauAsymmCalc class.
ClassImp(LauAbsCoeffSet)
UInt_t nNormPar_
Number of normalisation parameters (yields, asymmetries)
LauEmbeddedData * posSignalTree_
The B+ signal event tree.
Bool_t getReweightedEvent(LauIsobarDynamics *dynamics)
Retrieve an event from the data sample, applying an accept/reject based on the given DP model...
LauParameter()
Default constructor.
Definition: LauParameter.cc:30
Bool_t hasResonance(const TString &resName) const
Check whether this model includes a named resonance.
Double_t getm23() const
Get the m23 invariant mass.
Int_t curEvtCharge_
Current event charge.
void calcExtraFractions(Bool_t initValues=kFALSE)
Calculate the CP-conserving and CP-violating fit fractions.
const TString & name() const
The parameter name.
virtual LauParameter acp()=0
Calculate the CP asymmetry.
Double_t getmPrime() const
Get m&#39; value.
Class that defines the particular 3-body decay under study.
Definition: LauDaughters.hh:33
UInt_t addFitParameters(LauPdfList &pdfList)
Add parameters of the PDFs in the list to the list of all fit parameters.
Class for calculating the asymmetry between two variables.
Definition: LauAsymmCalc.hh:25
LauEmbeddedData * negSignalTree_
The B- signal event tree.
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.
File containing declaration of LauAbsCoeffSet class.
void clearUsedList()
Clear the list of used events.
Bool_t useNegReweighting_
Boolean to use reweighting for B-.
void calcExtraInfo(const Bool_t init=kFALSE)
Calculate the fit fractions, mean efficiency and total DP rate.
TH2 * trueHist(Int_t trueBin)
Retrieve the migration histogram for trueBin.
Definition: LauScfMap.cc:176
std::vector< LauComplex > negCoeffs_
The complex coefficients for B-.
LauBkgndPdfsList posBkgndPdfs_
The B+ background PDFs.
std::multimap< TString, std::pair< TString, TString > > TwoDMap
Type to associate the name of the species that have 2D PDFs with the names of the two variables invol...
Definition: LauSPlot.hh:68
LauParArray negFitFrac_
The B- fit fractions.
void listBinCentres(std::vector< Double_t > &xCoords, std::vector< Double_t > &yCoords) const
Create lists of the co-ordinates of the centres of true bins.
Definition: LauScfMap.cc:114
virtual void getEvtExtraLikelihoods(UInt_t iEvt)
Determine the signal and background likelihood for the extra variables for a given event...
Double_t sigDPLike_
Signal DP likelihood value.
File containing declaration of LauDaughters class.
const std::vector< Int_t > * trueBins(Int_t recoBin) const
Find which true bins contribute to the given reco bin.
Definition: LauScfMap.cc:155
LauPdfList negScfPdfs_
The B- SCF PDFs.
Bool_t useSCF_
Is the signal split into TM and SCF.
virtual Bool_t genExpt()
Toy MC generation and fitting overloaded functions.
virtual void setGenNtupleDoubleBranchValue(const TString &name, Double_t value)
Set the value of a double branch in the gen tree.
TString negParent_
Name of the parent particle.
TString posParent_
Name of the parent particle.
void embedNegSignal(const TString &fileName, const TString &treeName, Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment=kFALSE, Bool_t useReweighting=kFALSE)
Embed full simulation events for the B- signal, rather than generating toy from the PDFs...
File containing declaration of LauScfMap class.
Double_t getEvtthPrime() const
Retrieve the square Dalitz plot coordinate, theta&#39;, for the current event.
std::vector< Double_t > fakeSCFFracs_
The cached values of the SCF fraction for each bin centre.
TString getConjResName(const TString &resName) const
Retrieve the name of the charge conjugate of a named resonance.
std::vector< std::vector< LauParameter > > LauParArray
Type to define an array of parameters.
LauAbsEffModel * getEffModel()
Retrieve the model for the efficiency across the Dalitz plot.
Int_t fitStatus() const
Access the fit status information.
virtual void printParValues() const =0
Print the current values of the parameters.
Bool_t storeSignalMCMatch(LauEmbeddedData *embeddedData)
Store the MC truth info on the TM/SCF nature of the embedded signal event.
Bool_t usingBkgnd_
Background boolean.
File containing declaration of LauPrint class.
Pure abstract base class for defining the efficiency description across the signal Dalitz plot...
LauGenInfo eventsToGenerate()
Determine the number of events to generate for each hypothesis.
Double_t prodPdfValue(LauPdfList &pdfList, UInt_t iEvt)
Calculate the product of the per-event likelihoods of the PDFs in the list.
std::vector< LauAbsPdf * > LauPdfList
List of Pdfs.
std::map< TString, Double_t > LauFitData
Type for holding event data.
LauBkgndYieldList bkgndAsym_
Background asymmetries(s)
Double_t getm23Sq() const
Get the m23 invariant mass square.
LauIsobarDynamics * negSigModel_
The B- signal Dalitz plot model.
const LauParArray & getFitFractions() const
Retrieve the fit fractions for the amplitude components.
void recalculateNormalisation()
recalculate Normalization
Bool_t squareDP() const
Determine to use or not the square Dalitz plot.
Definition: LauDaughters.hh:78
LauParameter negDPRate_
The average DP rate for B-.
Double_t getEvtScfFraction() const
Retrieve the fraction of events that are poorly reconstructed (the self cross feed fraction) for the ...
Class to define various output print commands.
Definition: LauPrint.hh:29
Double_t getValue(const TString &name) const
Get the value of a specified branch.
LauDaughters * getDaughters()
Retrieve the daughters.
std::vector< Double_t > bkgndTotalLike_
Total background likelihood(s)
void calcLikelihoodInfo(const UInt_t iEvt)
Calculate the likelihood (and all associated information) for the given event number.
virtual void fillGenNtupleBranches()
Fill the gen tuple branches.
LauParameter getMeanEff() const
Retrieve the mean efficiency across the Dalitz plot.
std::vector< Double_t > recoSCFFracs_
The cached values of the SCF fraction for each event.
LauParameter posDPRate_
The average DP rate for B+.
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.
std::vector< Double_t > bkgndExtraLike_
Background likelihood value(s) from extra PDFs.
void setupGenNtupleBranches()
Setup the required ntuple branches.
void printFitParameters(const LauPdfList &pdfList, std::ostream &fout) const
Print the fit parameters for all PDFs in the list.
Double_t scfExtraLike_
SCF likelihood from extra PDFs.
Bool_t haveBranch(const TString &name) const
Boolean to determine whether branch exists.
File containing declaration of LauCPFitModel class.
Bool_t reuseSignal_
Boolean to reuse signal events.
LauParArray posFitFrac_
The B+ fit fractions.
virtual void cacheInputFitVars()
Read in the input fit data variables, e.g. m13Sq and m23Sq.
void setSignalDPParameters()
Set the fit parameters for the DP model.
Double_t getEvtmPrime() const
Retrieve the square Dalitz plot coordinate, m&#39;, for the current event.
void randomiseInitFitPars()
Randomise the initial fit parameters.
void appendBinCentres(LauFitDataTree *inputData)
Append fake data points to the inputData for each bin in the SCF smearing matrix. ...
File containing declaration of LauIsobarDynamics class.
virtual UInt_t index() const
Retrieve the index number of the coefficient set.
void updateKinematics(Double_t m13Sq, Double_t m23Sq)
Update all kinematic quantities based on the DP co-ordinates m13Sq and m23Sq.
void fillDataTree(const LauFitDataTree &fitDataTree)
Fill the internal data structure that caches the resonance dynamics.
virtual void checkInitFitParams()
Check the initial fit parameters.
const LauParameterList & extraPars() const
Access the extra variables.
Double_t abs2() const
Obtain the square of the absolute value of the complex number.
Definition: LauComplex.hh:232
virtual void recalculateNormalisation()
Recalculate Normalization the signal DP models.
std::vector< LauParameter * > LauParameterPList
List of parameter pointers.
const Bool_t tagged_
IS the analysis tagged?
UInt_t nExtraPdfPar_
Number of extra PDF parameters.
virtual void finaliseFitResults(const TString &tablePrefixName)
Get the fit results and store them.
Double_t getEvtm23Sq() const
Retrieve the invariant mass squared of the second and third daughters in the current event...
virtual void printFitFractions(std::ostream &output)
Print the fit fractions, total DP rate and mean efficiency.
Double_t calcEfficiency(const LauKinematics *kinematics) const
Determine the efficiency for a given point in the Dalitz plot.
Definition: LauEffModel.cc:261
LauParArray negFitFracEffUnCorr_
Fit B- fractions (uncorrected for the efficiency)
Class for defining a CP fit model.
Abstract interface to the fitting and toy MC model.
std::vector< LauParameter * > & getFloatingParameters()
Retrieve the floating parameters of the resonance models.
UInt_t iExpt() const
Obtain the number of the current experiment.
const TString tagVarName_
Event charge.
LauParameter * signalEvents_
Signal yield.
LauParArray posFitFracEffUnCorr_
Fit B+ fractions (uncorrected for the efficiency)
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< Int_t > evtCharges_
Vector to store event charges.
const LauFitNtuple * fitNtuple() const
Access the fit ntuple.
Int_t binNumber(Double_t xCoord, Double_t yCoord) const
Find the global bin number for the given co-ordinates.
Definition: LauScfMap.cc:144
Bool_t doPoissonSmearing() const
Determine whether Poisson smearing is enabled for the toy MC generation.
Double_t getc12() const
Get the cosine of the helicity angle theta12.
Double_t getm13() const
Get the m13 invariant mass.
std::vector< LauPdfList > LauBkgndPdfsList
Typedef for a vector of background PDFs.
void calcAsymmetries(Bool_t initValues=kFALSE)
Calculate the CP asymmetries.
std::map< TString, Double_t > NumbMap
Type to associate a category name with a double precision number, e.g. a yield or PDF value for a giv...
Definition: LauSPlot.hh:62
Bool_t doEMLFit() const
Determine whether an extended maximum likelihood fit it being performed.
LauParameter getDPRate() const
Retrieve the overall Dalitz plot rate.
LauScfMap * scfMap_
The smearing matrix for the SCF DP PDF.
File containing declaration of LauComplex class.
Double_t sigExtraLike_
Signal likelihood from extra PDFs.
LauFitData getValues(const std::vector< TString > &names) const
Get values of specified branches.
const LauParArray & getFitFractionsEfficiencyUncorrected() const
Retrieve the fit fractions for the amplitude components.
virtual void setNBkgndEvents(LauParameter *nBkgndEvents)
Set the background event yield(s)
virtual void initialise()
Initialise the fit.
virtual LauSPlot::NumbMap fixdSpeciesNames() const
Returns the names and yields of species that are fixed in the fit.
Class for defining the abstract interface for complex coefficient classes.
Double_t nll() const
Access the current NLL value.
Double_t sigTotalLike_
Total signal likelihood.
void storeCorrMatrix(UInt_t iExpt, Double_t NLL, Int_t fitStatus, const TMatrixD &covMatrix)
Store the correlation matrix and other fit information.
Definition: LauFitNtuple.cc:64
void clearFitParVectors()
Clear the vectors containing fit parameters.
virtual void setAmpCoeffSet(LauAbsCoeffSet *coeffSet)
Set the DP amplitude coefficients.
Class for defining the fit parameter objects.
Definition: LauParameter.hh:34
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.
UInt_t nSigComp_
Number of signal components.
virtual void setupSPlotNtupleBranches()
Add branches to store experiment number and the event number within the experiment.
void clearExtraVarVectors()
Clear the vectors containing extra ntuple variables.
LauParameter * signalAsym_
Signal asymmetry.
void getEmbeddedEvent(LauKinematics *kinematics)
Retrieve an event from the data sample.
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::vector< LauParameter > getExtraParameters()
Retrieve any extra parameters/quantities (e.g. K-matrix total fit fractions)
std::vector< Double_t > bkgndDPLike_
Background DP likelihood value(s)
void embedNegBkgnd(const TString &bgClass, 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...
LauBkgndEmbDataList negBkgndTree_
The B- background event tree.
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.
LauBkgndReuseEventsList reuseBkgnd_
Vector of booleans to reuse background events.
std::vector< LauAbsCoeffSet * > coeffPars_
Magnitudes and Phases.
LauParameter posMeanEff_
The mean efficiency for B+ model.
UInt_t bkgndClassID(const TString &className) const
The number assigned to a background class.
const LauFitData & getData(UInt_t iEvt) const
Retrieve the data for a given event.
File containing LauRandom namespace.
Bool_t generate()
Generate a toy MC signal event.
Bool_t useRandomInitFitPars() const
Determine whether the initial values of the fit parameters, in particular the isobar coefficient para...
File containing declaration of LauEffModel class.
const LauFitDataTree * fitData() const
Access the data store.
LauParArray CPCFitFrac_
The CP conserving fit fraction.
virtual void setNSigEvents(LauParameter *nSigEvents)
Set the signal event yield.
virtual void weightEvents(const TString &dataFileName, const TString &dataTreeName)
Weight events based on the DP model.
Bool_t validBkgndClass(const TString &className) const
Check if the given background class is in the list.
Bool_t usePosReweighting_
Boolean to use reweighting for B+.
Bool_t generateBkgndEvent(UInt_t bgID)
Generate background event.
virtual void writeOutTable(const TString &outputFile)
Write the fit results in latex table format.
Double_t getm12Sq() const
Get the m12 invariant mass square.
Class that implements the efficiency description across the signal Dalitz plot.
Definition: LauEffModel.hh:37
TString getNameParent() const
Get name of the parent particle.
const LauParameterPSet & resPars() const
Access the fit variables which affect the DP normalisation.
UInt_t nUsedEvents() const
Retrieve the number of events that have already been sampled.
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.
virtual void updateCoeffs()
Update the coefficients.
virtual void setSPlotNtupleIntegerBranchValue(const TString &name, Int_t value)
Set the value of an integer branch in the sPlot tree.
Double_t setSPlotNtupleBranchValues(LauPdfList *extraPdfs, const TString &prefix, UInt_t iEvt)
Set the branches for the sPlot ntuple with extra PDFs.
Double_t initValue() const
The initial value of the parameter.
virtual ~LauCPFitModel()
Destructor.
std::vector< LauComplex > posCoeffs_
The complex coefficients for B+.
virtual void initialiseDPModels()
Initialise the signal DP models.
void updateSigEvents()
Update the signal events after Minuit sets background parameters.
Class for defining signal dynamics using the isobar model.
Double_t getm12() const
Get the m12 invariant mass.
virtual void fillSPlotNtupleBranches()
Fill the sPlot tuple.
virtual void setupBkgndVectors()
Define the length of the background vectors.
File containing LauConstants namespace.
void setSignalPdfs(LauAbsPdf *negPdf, LauAbsPdf *posPdf)
Set the signal PDFs.
void setFitNEvents()
Set the initial yields.
virtual LauSPlot::TwoDMap twodimPDFs() const
Returns the species and variables for all 2D PDFs in the fit.
virtual LauSPlot::NumbMap freeSpeciesNames() const
Returns the names and yields of species that are free in the fit.
LauEffModel * scfFracHist_
The histogram giving the DP-dependence of the SCF fraction.
void printFormat(std::ostream &stream, Double_t value) const
Method to choose the printing format to a specified level of precision.
Definition: LauPrint.cc:32
ToyMCStatus checkToyMC(Bool_t printErrorMessages=kTRUE, Bool_t printInfoMessages=kFALSE)
Check the status of the toy MC generation.
Class to store the results from the fit into an ntuple.
Definition: LauFitNtuple.hh:42
UInt_t nSigDPPar_
Number of signal DP parameters.
Bool_t haveBranch(const TString &name) const
Check if the named branch is stored.
virtual void propagateParUpdates()
Calculate things that depend on the fit parameters after they have been updated by Minuit...
Class for representing the 4D smearing matrix for mis-reconstructed signal (self cross feed) ...
Definition: LauScfMap.hh:29
virtual void storePerEvtLlhds()
Store the per event likelihood values.
Double_t scfDPLike_
SCF DP likelihood value.
Double_t getThetaPrime() const
Get theta&#39; value.
void updatePull()
Call to update the bias and pull values.
Double_t getEvtEff() const
Retrieve the efficiency for the current event.
UInt_t nBkgndClasses() const
Returns the number of background classes.
const TString & bkgndClassName(UInt_t classID) const
Get the name of a background class from the number.
std::vector< Double_t > recoJacobians_
The cached values of the sqDP jacobians for each event.
Bool_t enableEmbedding() const
Determine whether embedding of events is enabled in the generation.
virtual Double_t getEventSum() const
Get the total number of events.
virtual TString name() const
Retrieve the name of the coefficient set.
Class for defining the abstract interface for PDF classes.
Definition: LauAbsPdf.hh:41
LauIsobarDynamics * posSigModel_
The B+ signal Dalitz plot model.
const TMatrixD & covarianceMatrix() const
Access the fit covariance matrix.
const LauParameterPList & fitPars() const
Access the fit variables.
LauPdfList posScfPdfs_
The B+ SCF PDFs.
Bool_t usingScfModel()
Check whether a self cross feed fraction model is being used.
LauBkgndYieldList bkgndEvents_
Background yield(s)
Class for calculating 3-body kinematic quantities.
Double_t value() const
The value of the parameter.
LauBkgndDPModelList negBkgndDPModels_
The B- background Dalitz plot models.
UInt_t nEvents() const
Retrieve the number of events.
LauParArray CPVFitFrac_
The CP violating fit fraction.
virtual void addGenNtupleDoubleBranch(const TString &name)
Add a branch to the gen tree for storing a double.
LauPdfList posSignalPdfs_
The B+ signal PDFs.
Double_t calcSqDPJacobian()
Calculate the Jacobian for the transformation m23^2, m13^2 -&gt; m&#39;, theta&#39; (square DP) ...
void updateFitParameters(LauPdfList &pdfList)
Update the fit parameters for all PDFs in the list.
void embedPosBkgnd(const TString &bgClass, 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...
std::vector< Double_t > fakeJacobians_
The cached values of the sqDP jacobians for each true bin.
virtual const TString & baseName() const
Retrieve the base name of the coefficient set.
Double_t getAsymmetry()
Retrieve the asymmetry.
Definition: LauAsymmCalc.hh:48
virtual Double_t getEvtSCFDPLikelihood(UInt_t iEvt)
Calculate the SCF likelihood for the DP for a given event.
Bool_t findBranches()
Find and read the branches in data tree.
File containing declaration of LauAbsBkgndDPModel class.
void setSCFPdfs(LauAbsPdf *negPdf, LauAbsPdf *posPdf)
Set the SCF PDF for a given variable.
std::map< std::pair< TString, Int_t >, std::pair< Int_t, Double_t > > LauGenInfo
Define a map to be used to store a category name and numbers.
Bool_t useSCFHist_
Is the SCF fraction DP-dependent.
LauParameter negMeanEff_
The mean efficiency for B- model.
void setBkgndPdfs(const TString &bkgndClass, LauAbsPdf *negPdf, LauAbsPdf *posPdf)
Set the background PDFs.
Double_t prob(Int_t recoBin, Int_t trueBin) const
Probability of a true event in the given true bin migrating to the reco bin.
Definition: LauScfMap.cc:165
void initialise(const std::vector< LauComplex > &coeffs)
Initialise the Dalitz plot dynamics.
virtual Double_t getTotEvtLikelihood(UInt_t iEvt)
Get the total likelihood for each event.
virtual void getEvtDPLikelihood(UInt_t iEvt)
Calculate the signal and background likelihoods for the DP for a given event.
Double_t getEvtm13Sq() const
Retrieve the invariant mass squared of the first and third daughters in the current event...
Double_t genValue() const
The value generated for the parameter.
LauBkgndEmbDataList posBkgndTree_
The B+ background event tree.
void setBkgndDPModels(const TString &bkgndClass, LauAbsBkgndDPModel *negModel, LauAbsBkgndDPModel *posModel)
Set the background DP models.
LauBkgndPdfsList negBkgndPdfs_
The B- background PDFs.
Class to store the input fit variables.
LauParameter scfFrac_
The (global) SCF fraction parameter.
std::set< TString > NameSet
Type to store names, e.g. of the discriminating/control variables.
Definition: LauSPlot.hh:59
File containing declaration of LauFitNtuple class.
void setDPBranchValues()
Store all of the DP information.
LauKinematics * posKinematics_
The B+ Dalitz plot kinematics object.
void embedPosSignal(const TString &fileName, const TString &treeName, Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment=kFALSE, Bool_t useReweighting=kFALSE)
Embed full simulation events for the B+ signal, rather than generating toy from the PDFs...
std::vector< LauParameter > fitFracAsymm_
The fit fraction asymmetries.
virtual Bool_t generate()=0
Generate a toy MC event from the model.
void setExtraNtupleVars()
Set-up other parameters that are derived from the fit results, e.g. fit fractions.
virtual Bool_t passVeto(const LauKinematics *kinematics) const =0
Determine whether the given DP position is outside the vetoes.