laura is hosted by Hepforge, IPPP Durham
Laura++  v2r2
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 - 2013.
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 "LauAbsDPDynamics.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(LauAbsDPDynamics* negModel, LauAbsDPDynamics* 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  // First of all check that, if a given component is being used,
422  // we've got the PDFs for all the variables involved
423  if ( this->useDP() ) {
424  if ((negSigModel_ == 0) || (posSigModel_ == 0)) {
425  std::cerr << "ERROR in LauCPFitModel::initialise : the pointer to one (neg or pos) of the signal DP models is null.\n";
426  std::cerr << " : Removing the Dalitz Plot from the model." << std::endl;
427  this->useDP(kFALSE);
428  }
429  if ( usingBkgnd_ ) {
430  if ( negBkgndDPModels_.empty() || posBkgndDPModels_.empty() ) {
431  std::cerr << "ERROR in LauCPFitModel::initialise : No background DP models found.\n";
432  std::cerr << " : Removing the Dalitz plot from the model." << std::endl;
433  this->useDP(kFALSE);
434  }
435  for (LauBkgndDPModelList::const_iterator dpmodel_iter = negBkgndDPModels_.begin(); dpmodel_iter != negBkgndDPModels_.end(); ++dpmodel_iter ) {
436  if ( (*dpmodel_iter) == 0 ) {
437  std::cerr << "ERROR in LauCPFitModel::initialise : The pointer to one of the background DP models is null.\n";
438  std::cerr << " : Removing the Dalitz Plot from the model." << std::endl;
439  this->useDP(kFALSE);
440  break;
441  }
442  }
443  for (LauBkgndDPModelList::const_iterator dpmodel_iter = posBkgndDPModels_.begin(); dpmodel_iter != posBkgndDPModels_.end(); ++dpmodel_iter ) {
444  if ( (*dpmodel_iter) == 0 ) {
445  std::cerr << "ERROR in LauCPFitModel::initialise : The pointer to one of the background DP models is null.\n";
446  std::cerr << " : Removing the Dalitz Plot from the model." << std::endl;
447  this->useDP(kFALSE);
448  break;
449  }
450  }
451  }
452  }
453 
454  // Next check that, if a given component is being used we've got the
455  // right number of PDFs for all the variables involved
456  // TODO - should probably check variable names and so on as well
457 
458  UInt_t nsigpdfvars(0);
459  for ( LauPdfList::const_iterator pdf_iter = negSignalPdfs_.begin(); pdf_iter != negSignalPdfs_.end(); ++pdf_iter ) {
460  std::vector<TString> varNames = (*pdf_iter)->varNames();
461  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
462  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
463  ++nsigpdfvars;
464  }
465  }
466  }
467  if (useSCF_) {
468  UInt_t nscfpdfvars(0);
469  for ( LauPdfList::const_iterator pdf_iter = negScfPdfs_.begin(); pdf_iter != negScfPdfs_.end(); ++pdf_iter ) {
470  std::vector<TString> varNames = (*pdf_iter)->varNames();
471  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
472  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
473  ++nscfpdfvars;
474  }
475  }
476  }
477  if (nscfpdfvars != nsigpdfvars) {
478  std::cerr << "ERROR in LauCPFitModel::initialise : There are " << nsigpdfvars << " TM signal PDF variables but " << nscfpdfvars << " SCF signal PDF variables." << std::endl;
479  gSystem->Exit(EXIT_FAILURE);
480  }
481  }
482  if (usingBkgnd_) {
483  for (LauBkgndPdfsList::const_iterator bgclass_iter = negBkgndPdfs_.begin(); bgclass_iter != negBkgndPdfs_.end(); ++bgclass_iter) {
484  UInt_t nbkgndpdfvars(0);
485  const LauPdfList& pdfList = (*bgclass_iter);
486  for ( LauPdfList::const_iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter ) {
487  std::vector<TString> varNames = (*pdf_iter)->varNames();
488  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
489  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
490  ++nbkgndpdfvars;
491  }
492  }
493  }
494  if (nbkgndpdfvars != nsigpdfvars) {
495  std::cerr << "ERROR in LauCPFitModel::initialise : There are " << nsigpdfvars << " signal PDF variables but " << nbkgndpdfvars << " bkgnd PDF variables." << std::endl;
496  gSystem->Exit(EXIT_FAILURE);
497  }
498  }
499  }
500 
501  // Clear the vectors of parameter information so we can start from scratch
502  this->clearFitParVectors();
503 
504  // Set the fit parameters for signal and background models
505  this->setSignalDPParameters();
506 
507  // Set the fit parameters for the various extra PDFs
508  this->setExtraPdfParameters();
509 
510  // Set the initial bg and signal events
511  this->setFitNEvents();
512 
513  // Check that we have the expected number of fit variables
514  const LauParameterPList& fitVars = this->fitPars();
515  if (fitVars.size() != (nSigDPPar_ + nExtraPdfPar_ + nNormPar_)) {
516  std::cerr << "ERROR in LauCPFitModel::initialise : Number of fit parameters not of expected size. Exiting" << std::endl;
517  gSystem->Exit(EXIT_FAILURE);
518  }
519 
520  // From the initial parameter values calculate the coefficients
521  // so they can be passed to the signal model
522  this->updateCoeffs();
523 
524  // Initialisation
525  if (this->useDP() == kTRUE) {
526  this->initialiseDPModels();
527  }
528 
529  if (!this->useDP() && negSignalPdfs_.empty()) {
530  std::cerr << "ERROR in LauCPFitModel::initialise : Signal model doesn't exist for any variable." << std::endl;
531  gSystem->Exit(EXIT_FAILURE);
532  }
533 
534  this->setExtraNtupleVars();
535 }
536 
538 {
539  std::cout << "INFO in LauCPFitModel::initialiseDPModels : Initialising signal DP model" << std::endl;
542 
543  if (usingBkgnd_ == kTRUE) {
544  for (LauBkgndDPModelList::iterator iter = negBkgndDPModels_.begin(); iter != negBkgndDPModels_.end(); ++iter) {
545  (*iter)->initialise();
546  }
547  for (LauBkgndDPModelList::iterator iter = posBkgndDPModels_.begin(); iter != posBkgndDPModels_.end(); ++iter) {
548  (*iter)->initialise();
549  }
550  }
551 }
552 
554 {
555  // Set the fit parameters for the signal model.
556 
557  nSigDPPar_ = 0;
558  if ( ! this->useDP() ) {
559  return;
560  }
561 
562  std::cout << "INFO in LauCPFitModel::setSignalDPParameters : Setting the initial fit parameters for the signal DP model." << std::endl;
563 
564  // Need to check that the number of components we have and that the dynamics has matches up
565  UInt_t nNegAmp = negSigModel_->getnAmp();
566  UInt_t nPosAmp = posSigModel_->getnAmp();
567  if ( nNegAmp != nPosAmp ) {
568  std::cerr << "ERROR in LauCPFitModel::setSignalDPParameters : Unequal number of signal DP components in the negative and positive models: " << nNegAmp << " != " << nPosAmp << std::endl;
569  gSystem->Exit(EXIT_FAILURE);
570  }
571  if ( nNegAmp != nSigComp_ ) {
572  std::cerr << "ERROR in LauCPFitModel::setSignalDPParameters : Number of signal DP components in the model (" << nNegAmp << ") not equal to number of coefficients supplied (" << nSigComp_ << ")." << std::endl;
573  gSystem->Exit(EXIT_FAILURE);
574  }
575 
576 
577  // Place signal model parameters in vector of fit variables
578  LauParameterPList& fitVars = this->fitPars();
579  for (UInt_t i = 0; i < nSigComp_; i++) {
580  LauParameterPList pars = coeffPars_[i]->getParameters();
581  for (LauParameterPList::iterator iter = pars.begin(); iter != pars.end(); ++iter) {
582  if ( !(*iter)->clone() ) {
583  fitVars.push_back(*iter);
584  ++nSigDPPar_;
585  }
586  }
587  }
588 }
589 
591 {
592  // Include all the parameters of the PDF in the fit
593  // NB all of them are passed to the fit, even though some have been fixed through parameter.fixed(kTRUE)
594  // With the new "cloned parameter" scheme only "original" parameters are passed to the fit.
595  // Their clones are updated automatically when the originals are updated.
596 
597  nExtraPdfPar_ = 0;
598 
600  if ( tagged_ ) {
602  }
603 
604  if (useSCF_ == kTRUE) {
606  if ( tagged_ ) {
608  }
609  }
610 
611  if (usingBkgnd_ == kTRUE) {
612  for (LauBkgndPdfsList::iterator iter = negBkgndPdfs_.begin(); iter != negBkgndPdfs_.end(); ++iter) {
613  nExtraPdfPar_ += this->addFitParameters(*iter);
614  }
615  if ( tagged_ ) {
616  for (LauBkgndPdfsList::iterator iter = posBkgndPdfs_.begin(); iter != posBkgndPdfs_.end(); ++iter) {
617  nExtraPdfPar_ += this->addFitParameters(*iter);
618  }
619  }
620  }
621 }
622 
624 {
625  if ( signalEvents_ == 0 ) {
626  std::cerr << "ERROR in LauCPFitModel::setFitNEvents : Signal yield not defined." << std::endl;
627  return;
628  }
629  nNormPar_ = 0;
630 
631  // initialise the total number of events to be the sum of all the hypotheses
632  Double_t nTotEvts = signalEvents_->value();
633  for (LauBkgndYieldList::const_iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
634  nTotEvts += (*iter)->value();
635  if ( (*iter) == 0 ) {
636  std::cerr << "ERROR in LauCPFitModel::setFitNEvents : Background yield not defined." << std::endl;
637  return;
638  }
639  }
640  this->eventsPerExpt(TMath::FloorNint(nTotEvts));
641 
642  LauParameterPList& fitVars = this->fitPars();
643 
644  // if doing an extended ML fit add the number of signal events into the fit parameters
645  if (this->doEMLFit()) {
646  std::cout << "INFO in LauCPFitModel::setFitNEvents : Initialising number of events for signal and background components..." << std::endl;
647  // add the signal fraction to the list of fit parameters
648  fitVars.push_back(signalEvents_);
649  ++nNormPar_;
650  } else {
651  std::cout << "INFO in LauCPFitModel::setFitNEvents : Initialising number of events for background components (and hence signal)..." << std::endl;
652  }
653 
654  // if not using the DP in the model we need an explicit signal asymmetry parameter
655  if (this->useDP() == kFALSE) {
656  fitVars.push_back(signalAsym_);
657  ++nNormPar_;
658  }
659 
660  if (useSCF_ && !useSCFHist_) {
661  fitVars.push_back(&scfFrac_);
662  ++nNormPar_;
663  }
664 
665  if (usingBkgnd_ == kTRUE) {
666  for (LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
667  LauParameter* parameter = (*iter);
668  fitVars.push_back(parameter);
669  ++nNormPar_;
670  }
671  for (LauBkgndYieldList::iterator iter = bkgndAsym_.begin(); iter != bkgndAsym_.end(); ++iter) {
672  LauParameter* parameter = (*iter);
673  fitVars.push_back(parameter);
674  ++nNormPar_;
675  }
676  }
677 }
678 
680 {
681  // Set-up other parameters derived from the fit results, e.g. fit fractions.
682 
683  if (this->useDP() != kTRUE) {
684  return;
685  }
686 
687  // First clear the vectors so we start from scratch
688  this->clearExtraVarVectors();
689 
690  LauParameterList& extraVars = this->extraPars();
691 
692  // Add the positive and negative fit fractions for each signal component
694  if (negFitFrac_.size() != nSigComp_) {
695  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << negFitFrac_.size() << std::endl;
696  gSystem->Exit(EXIT_FAILURE);
697  }
698  for (UInt_t i(0); i<nSigComp_; ++i) {
699  if (negFitFrac_[i].size() != nSigComp_) {
700  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << negFitFrac_[i].size() << std::endl;
701  gSystem->Exit(EXIT_FAILURE);
702  }
703  }
704 
706  if (posFitFrac_.size() != nSigComp_) {
707  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << posFitFrac_.size() << std::endl;
708  gSystem->Exit(EXIT_FAILURE);
709  }
710  for (UInt_t i(0); i<nSigComp_; ++i) {
711  if (posFitFrac_[i].size() != nSigComp_) {
712  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << posFitFrac_[i].size() << std::endl;
713  gSystem->Exit(EXIT_FAILURE);
714  }
715  }
716 
717  // Add the positive and negative fit fractions that have not been corrected for the efficiency for each signal component
719  if (negFitFracEffUnCorr_.size() != nSigComp_) {
720  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << negFitFracEffUnCorr_.size() << std::endl;
721  gSystem->Exit(EXIT_FAILURE);
722  }
723  for (UInt_t i(0); i<nSigComp_; ++i) {
724  if (negFitFracEffUnCorr_[i].size() != nSigComp_) {
725  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << negFitFracEffUnCorr_[i].size() << std::endl;
726  gSystem->Exit(EXIT_FAILURE);
727  }
728  }
729 
731  if (posFitFracEffUnCorr_.size() != nSigComp_) {
732  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << posFitFracEffUnCorr_.size() << std::endl;
733  gSystem->Exit(EXIT_FAILURE);
734  }
735  for (UInt_t i(0); i<nSigComp_; ++i) {
736  if (posFitFracEffUnCorr_[i].size() != nSigComp_) {
737  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << posFitFracEffUnCorr_[i].size() << std::endl;
738  gSystem->Exit(EXIT_FAILURE);
739  }
740  }
741 
742  for (UInt_t i(0); i<nSigComp_; ++i) {
743  for (UInt_t j(i); j<nSigComp_; ++j) {
744  TString name = negFitFrac_[i][j].name();
745  name.Insert( name.Index("FitFrac"), "Neg" );
746  negFitFrac_[i][j].name(name);
747  extraVars.push_back(negFitFrac_[i][j]);
748 
749  name = posFitFrac_[i][j].name();
750  name.Insert( name.Index("FitFrac"), "Pos" );
751  posFitFrac_[i][j].name(name);
752  extraVars.push_back(posFitFrac_[i][j]);
753 
754  name = negFitFracEffUnCorr_[i][j].name();
755  name.Insert( name.Index("FitFrac"), "Neg" );
756  negFitFracEffUnCorr_[i][j].name(name);
757  extraVars.push_back(negFitFracEffUnCorr_[i][j]);
758 
759  name = posFitFracEffUnCorr_[i][j].name();
760  name.Insert( name.Index("FitFrac"), "Pos" );
761  posFitFracEffUnCorr_[i][j].name(name);
762  extraVars.push_back(posFitFracEffUnCorr_[i][j]);
763  }
764  }
765 
766  // Store any extra parameters/quantities from the DP model (e.g. K-matrix total fit fractions)
767  std::vector<LauParameter> negExtraPars = negSigModel_->getExtraParameters();
768  std::vector<LauParameter>::iterator negExtraIter;
769  for (negExtraIter = negExtraPars.begin(); negExtraIter != negExtraPars.end(); ++negExtraIter) {
770  LauParameter negExtraParameter = (*negExtraIter);
771  extraVars.push_back(negExtraParameter);
772  }
773 
774  std::vector<LauParameter> posExtraPars = posSigModel_->getExtraParameters();
775  std::vector<LauParameter>::iterator posExtraIter;
776  for (posExtraIter = posExtraPars.begin(); posExtraIter != posExtraPars.end(); ++posExtraIter) {
777  LauParameter posExtraParameter = (*posExtraIter);
778  extraVars.push_back(posExtraParameter);
779  }
780 
781  // Now add in the DP efficiency value
782  Double_t initMeanEff = negSigModel_->getMeanEff().initValue();
783  negMeanEff_.value(initMeanEff);
784  negMeanEff_.genValue(initMeanEff);
785  negMeanEff_.initValue(initMeanEff);
786  extraVars.push_back(negMeanEff_);
787 
788  initMeanEff = posSigModel_->getMeanEff().initValue();
789  posMeanEff_.value(initMeanEff);
790  posMeanEff_.genValue(initMeanEff);
791  posMeanEff_.initValue(initMeanEff);
792  extraVars.push_back(posMeanEff_);
793 
794  // Also add in the DP rates
795  Double_t initDPRate = negSigModel_->getDPRate().initValue();
796  negDPRate_.value(initDPRate);
797  negDPRate_.genValue(initDPRate);
798  negDPRate_.initValue(initDPRate);
799  extraVars.push_back(negDPRate_);
800 
801  initDPRate = posSigModel_->getDPRate().initValue();
802  posDPRate_.value(initDPRate);
803  posDPRate_.genValue(initDPRate);
804  posDPRate_.initValue(initDPRate);
805  extraVars.push_back(posDPRate_);
806 
807  // Calculate the CPC and CPV Fit Fractions, ACPs and FitFrac asymmetries
808  this->calcExtraFractions(kTRUE);
809  this->calcAsymmetries(kTRUE);
810 
811  // Add the CP violating and CP conserving fit fractions for each signal component
812  for (UInt_t i = 0; i < nSigComp_; i++) {
813  for (UInt_t j = i; j < nSigComp_; j++) {
814  extraVars.push_back(CPVFitFrac_[i][j]);
815  }
816  }
817  for (UInt_t i = 0; i < nSigComp_; i++) {
818  for (UInt_t j = i; j < nSigComp_; j++) {
819  extraVars.push_back(CPCFitFrac_[i][j]);
820  }
821  }
822 
823  // Add the Fit Fraction asymmetry for each signal component
824  for (UInt_t i = 0; i < nSigComp_; i++) {
825  extraVars.push_back(fitFracAsymm_[i]);
826  }
827 
828  // Add the calculated CP asymmetry for each signal component
829  for (UInt_t i = 0; i < nSigComp_; i++) {
830  extraVars.push_back(acp_[i]);
831  }
832 }
833 
834 void LauCPFitModel::calcExtraFractions(Bool_t initValues)
835 {
836  // Calculate the CP-conserving and CP-violating fit fractions
837 
838  if (initValues) {
839  // create the structure
840  CPCFitFrac_.clear();
841  CPVFitFrac_.clear();
842  CPCFitFrac_.resize(nSigComp_);
843  CPVFitFrac_.resize(nSigComp_);
844  for (UInt_t i(0); i<nSigComp_; ++i) {
845  CPCFitFrac_[i].resize(nSigComp_);
846  CPVFitFrac_[i].resize(nSigComp_);
847 
848  for (UInt_t j(i); j<nSigComp_; ++j) {
849  TString name = negFitFrac_[i][j].name();
850  name.Replace( name.Index("Neg"), 3, "CPC" );
851  CPCFitFrac_[i][j].name( name );
852  CPCFitFrac_[i][j].valueAndRange( 0.0, -100.0, 100.0 );
853 
854  name = negFitFrac_[i][j].name();
855  name.Replace( name.Index("Neg"), 3, "CPV" );
856  CPVFitFrac_[i][j].name( name );
857  CPVFitFrac_[i][j].valueAndRange( 0.0, -100.0, 100.0 );
858  }
859  }
860  }
861 
862  Double_t denom = negDPRate_ + posDPRate_;
863 
864  for (UInt_t i(0); i<nSigComp_; ++i) {
865  for (UInt_t j(i); j<nSigComp_; ++j) {
866 
867  Double_t negTerm = negFitFrac_[i][j]*negDPRate_;
868  Double_t posTerm = posFitFrac_[i][j]*posDPRate_;
869 
870  Double_t cpcFitFrac = (negTerm + posTerm)/denom;
871  Double_t cpvFitFrac = (negTerm - posTerm)/denom;
872 
873  CPCFitFrac_[i][j] = cpcFitFrac;
874  CPVFitFrac_[i][j] = cpvFitFrac;
875 
876  if (initValues) {
877  CPCFitFrac_[i][j].genValue(cpcFitFrac);
878  CPCFitFrac_[i][j].initValue(cpcFitFrac);
879 
880  CPVFitFrac_[i][j].genValue(cpvFitFrac);
881  CPVFitFrac_[i][j].initValue(cpvFitFrac);
882  }
883  }
884  }
885 }
886 
887 void LauCPFitModel::calcAsymmetries(Bool_t initValues)
888 {
889  // Calculate the CP asymmetries
890  // Also calculate the fit fraction asymmetries
891 
892  for (UInt_t i = 0; i < nSigComp_; i++) {
893 
894  acp_[i] = coeffPars_[i]->acp();
895 
896  LauAsymmCalc asymmCalc(negFitFrac_[i][i].value(), posFitFrac_[i][i].value());
897  Double_t asym = asymmCalc.getAsymmetry();
898  fitFracAsymm_[i] = asym;
899  if (initValues) {
900  fitFracAsymm_[i].genValue(asym);
901  fitFracAsymm_[i].initValue(asym);
902  }
903  }
904 }
905 
906 void LauCPFitModel::finaliseFitResults(const TString& tablePrefixName)
907 {
908  // Retrieve parameters from the fit results for calculations and toy generation
909  // and eventually store these in output root ntuples/text files
910 
911  // Now take the fit parameters and update them as necessary
912  // i.e. to make mag > 0.0, phase in the right range.
913  // This function will also calculate any other values, such as the
914  // fit fractions, using any errors provided by fitParErrors as appropriate.
915  // Also obtain the pull values: (measured - generated)/(average error)
916 
917  if (this->useDP() == kTRUE) {
918  for (UInt_t i = 0; i < nSigComp_; ++i) {
919  // Check whether we have "a/b > 0.0", and phases in the right range
920  coeffPars_[i]->finaliseValues();
921  }
922  }
923 
924  // update the pulls on the event fractions and asymmetries
925  if (this->doEMLFit()) {
927  }
928  if (this->useDP() == kFALSE) {
930  }
931  if (useSCF_ && !useSCFHist_) {
933  }
934  if (usingBkgnd_ == kTRUE) {
935  for (LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
936  (*iter)->updatePull();
937  }
938  for (LauBkgndYieldList::iterator iter = bkgndAsym_.begin(); iter != bkgndAsym_.end(); ++iter) {
939  (*iter)->updatePull();
940  }
941  }
942 
943  // Update the pulls on all the extra PDFs' parameters
946  if (useSCF_ == kTRUE) {
949  }
950  if (usingBkgnd_ == kTRUE) {
951  for (LauBkgndPdfsList::iterator iter = negBkgndPdfs_.begin(); iter != negBkgndPdfs_.end(); ++iter) {
952  this->updateFitParameters(*iter);
953  }
954  for (LauBkgndPdfsList::iterator iter = posBkgndPdfs_.begin(); iter != posBkgndPdfs_.end(); ++iter) {
955  this->updateFitParameters(*iter);
956  }
957  }
958 
959  // Fill the fit results to the ntuple
960 
961  // update the coefficients and then calculate the fit fractions and ACP's
962  if (this->useDP() == kTRUE) {
963  this->updateCoeffs();
966 
967  LauParArray negFitFrac = negSigModel_->getFitFractions();
968  if (negFitFrac.size() != nSigComp_) {
969  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << negFitFrac.size() << std::endl;
970  gSystem->Exit(EXIT_FAILURE);
971  }
972  for (UInt_t i(0); i<nSigComp_; ++i) {
973  if (negFitFrac[i].size() != nSigComp_) {
974  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << negFitFrac[i].size() << std::endl;
975  gSystem->Exit(EXIT_FAILURE);
976  }
977  }
978 
979  LauParArray posFitFrac = posSigModel_->getFitFractions();
980  if (posFitFrac.size() != nSigComp_) {
981  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << posFitFrac.size() << std::endl;
982  gSystem->Exit(EXIT_FAILURE);
983  }
984  for (UInt_t i(0); i<nSigComp_; ++i) {
985  if (posFitFrac[i].size() != nSigComp_) {
986  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << posFitFrac[i].size() << std::endl;
987  gSystem->Exit(EXIT_FAILURE);
988  }
989  }
990 
992  if (negFitFracEffUnCorr.size() != nSigComp_) {
993  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << negFitFracEffUnCorr.size() << std::endl;
994  gSystem->Exit(EXIT_FAILURE);
995  }
996  for (UInt_t i(0); i<nSigComp_; ++i) {
997  if (negFitFracEffUnCorr[i].size() != nSigComp_) {
998  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << negFitFracEffUnCorr[i].size() << std::endl;
999  gSystem->Exit(EXIT_FAILURE);
1000  }
1001  }
1002 
1004  if (posFitFracEffUnCorr.size() != nSigComp_) {
1005  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << posFitFracEffUnCorr.size() << std::endl;
1006  gSystem->Exit(EXIT_FAILURE);
1007  }
1008  for (UInt_t i(0); i<nSigComp_; ++i) {
1009  if (posFitFracEffUnCorr[i].size() != nSigComp_) {
1010  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << posFitFracEffUnCorr[i].size() << std::endl;
1011  gSystem->Exit(EXIT_FAILURE);
1012  }
1013  }
1014 
1015  for (UInt_t i(0); i<nSigComp_; ++i) {
1016  for (UInt_t j(i); j<nSigComp_; ++j) {
1017  negFitFrac_[i][j].value(negFitFrac[i][j].value());
1018  posFitFrac_[i][j].value(posFitFrac[i][j].value());
1019  negFitFracEffUnCorr_[i][j].value(negFitFracEffUnCorr[i][j].value());
1020  posFitFracEffUnCorr_[i][j].value(posFitFracEffUnCorr[i][j].value());
1021  }
1022  }
1023 
1028 
1029  this->calcExtraFractions();
1030  this->calcAsymmetries();
1031 
1032  // Then store the final fit parameters, and any extra parameters for
1033  // the signal model (e.g. fit fractions, FF asymmetries, ACPs, mean efficiency and DP rate)
1034 
1035  this->clearExtraVarVectors();
1036  LauParameterList& extraVars = this->extraPars();
1037 
1038  // Add the positive and negative fit fractions for each signal component
1039  for (UInt_t i(0); i<nSigComp_; ++i) {
1040  for (UInt_t j(i); j<nSigComp_; ++j) {
1041  extraVars.push_back(negFitFrac_[i][j]);
1042  extraVars.push_back(posFitFrac_[i][j]);
1043  extraVars.push_back(negFitFracEffUnCorr_[i][j]);
1044  extraVars.push_back(posFitFracEffUnCorr_[i][j]);
1045  }
1046  }
1047 
1048  // Store any extra parameters/quantities from the DP model (e.g. K-matrix total fit fractions)
1049  std::vector<LauParameter> negExtraPars = negSigModel_->getExtraParameters();
1050  std::vector<LauParameter>::iterator negExtraIter;
1051  for (negExtraIter = negExtraPars.begin(); negExtraIter != negExtraPars.end(); ++negExtraIter) {
1052  LauParameter negExtraParameter = (*negExtraIter);
1053  extraVars.push_back(negExtraParameter);
1054  }
1055 
1056  std::vector<LauParameter> posExtraPars = posSigModel_->getExtraParameters();
1057  std::vector<LauParameter>::iterator posExtraIter;
1058  for (posExtraIter = posExtraPars.begin(); posExtraIter != posExtraPars.end(); ++posExtraIter) {
1059  LauParameter posExtraParameter = (*posExtraIter);
1060  extraVars.push_back(posExtraParameter);
1061  }
1062 
1063  extraVars.push_back(negMeanEff_);
1064  extraVars.push_back(posMeanEff_);
1065  extraVars.push_back(negDPRate_);
1066  extraVars.push_back(posDPRate_);
1067 
1068  for (UInt_t i = 0; i < nSigComp_; i++) {
1069  for (UInt_t j(i); j<nSigComp_; ++j) {
1070  extraVars.push_back(CPVFitFrac_[i][j]);
1071  }
1072  }
1073  for (UInt_t i = 0; i < nSigComp_; i++) {
1074  for (UInt_t j(i); j<nSigComp_; ++j) {
1075  extraVars.push_back(CPCFitFrac_[i][j]);
1076  }
1077  }
1078  for (UInt_t i(0); i<nSigComp_; ++i) {
1079  extraVars.push_back(fitFracAsymm_[i]);
1080  }
1081  for (UInt_t i(0); i<nSigComp_; ++i) {
1082  extraVars.push_back(acp_[i]);
1083  }
1084 
1085  this->printFitFractions(std::cout);
1086  this->printAsymmetries(std::cout);
1087  }
1088 
1089  const LauParameterPList& fitVars = this->fitPars();
1090  const LauParameterList& extraVars = this->extraPars();
1091  LauFitNtuple* ntuple = this->fitNtuple();
1092  ntuple->storeParsAndErrors(fitVars, extraVars);
1093 
1094  // find out the correlation matrix for the parameters
1095  ntuple->storeCorrMatrix(this->iExpt(), this->nll(), this->fitStatus(), this->covarianceMatrix());
1096 
1097  // Fill the data into ntuple
1098  ntuple->updateFitNtuple();
1099 
1100  // Print out the partial fit fractions, phases and the
1101  // averaged efficiency, reweighted by the dynamics (and anything else)
1102  if (this->writeLatexTable()) {
1103  TString sigOutFileName(tablePrefixName);
1104  sigOutFileName += "_"; sigOutFileName += this->iExpt(); sigOutFileName += "Expt.tex";
1105  this->writeOutTable(sigOutFileName);
1106  }
1107 }
1108 
1109 void LauCPFitModel::printFitFractions(std::ostream& output)
1110 {
1111  // Print out Fit Fractions, total DP rate and mean efficiency
1112  // First for the B- events
1113  for (UInt_t i = 0; i < nSigComp_; i++) {
1114  const TString compName(coeffPars_[i]->name());
1115  output << negParent_ << " FitFraction for component " << i << " (" << compName << ") = " << negFitFrac_[i][i] << std::endl;
1116  }
1117  output << negParent_ << " overall DP rate (integral of matrix element squared) = " << negDPRate_ << std::endl;
1118  output << negParent_ << " average efficiency weighted by whole DP dynamics = " << negMeanEff_ << std::endl;
1119 
1120  // Then for the positive sample
1121  for (UInt_t i = 0; i < nSigComp_; i++) {
1122  const TString compName(coeffPars_[i]->name());
1123  const TString conjName(negSigModel_->getConjResName(compName));
1124  output << posParent_ << " FitFraction for component " << i << " (" << conjName << ") = " << posFitFrac_[i][i] << std::endl;
1125  }
1126  output << posParent_ << " overall DP rate (integral of matrix element squared) = " << posDPRate_ << std::endl;
1127  output << posParent_ << " average efficiency weighted by whole DP dynamics = " << posMeanEff_ << std::endl;
1128 }
1129 
1130 void LauCPFitModel::printAsymmetries(std::ostream& output)
1131 {
1132  for (UInt_t i = 0; i < nSigComp_; i++) {
1133  const TString compName(coeffPars_[i]->name());
1134  output << "Fit Fraction asymmetry for component " << i << " (" << compName << ") = " << fitFracAsymm_[i] << std::endl;
1135  }
1136  for (UInt_t i = 0; i < nSigComp_; i++) {
1137  const TString compName(coeffPars_[i]->name());
1138  output << "ACP for component " << i << " (" << compName << ") = " << acp_[i] << std::endl;
1139  }
1140 }
1141 
1142 void LauCPFitModel::writeOutTable(const TString& outputFile)
1143 {
1144  // Write out the results of the fit to a tex-readable table
1145  // TODO - need to include the yields in this table
1146  std::ofstream fout(outputFile);
1147  LauPrint print;
1148 
1149  std::cout << "INFO in LauCPFitModel::writeOutTable : Writing out results of the fit to the tex file " << outputFile << std::endl;
1150 
1151  if (this->useDP() == kTRUE) {
1152  // print the fit coefficients in one table
1153  coeffPars_.front()->printTableHeading(fout);
1154  for (UInt_t i = 0; i < nSigComp_; i++) {
1155  coeffPars_[i]->printTableRow(fout);
1156  }
1157  fout << "\\hline" << std::endl;
1158  fout << "\\end{tabular}" << std::endl << std::endl;
1159 
1160  // print the fit fractions and asymmetries in another
1161  fout << "\\begin{tabular}{|l|c|c|c|c|}" << std::endl;
1162  fout << "\\hline" << std::endl;
1163  fout << "Component & " << negParent_ << " Fit Fraction & " << posParent_ << " Fit Fraction & Fit Fraction Asymmetry & ACP \\\\" << std::endl;
1164  fout << "\\hline" << std::endl;
1165 
1166  Double_t negFitFracSum(0.0);
1167  Double_t posFitFracSum(0.0);
1168 
1169  for (UInt_t i = 0; i < nSigComp_; i++) {
1170  TString resName = coeffPars_[i]->name();
1171  resName = resName.ReplaceAll("_", "\\_");
1172 
1173  Double_t negFitFrac = negFitFrac_[i][i].value();
1174  Double_t posFitFrac = posFitFrac_[i][i].value();
1175  negFitFracSum += negFitFrac;
1176  posFitFracSum += posFitFrac;
1177 
1178  Double_t fitFracAsymm = fitFracAsymm_[i].value();
1179 
1180  Double_t acp = acp_[i].value();
1181  Double_t acpErr = acp_[i].error();
1182 
1183  fout << resName << " & $";
1184 
1185  print.printFormat(fout, negFitFrac);
1186  fout << "$ & $";
1187  print.printFormat(fout, posFitFrac);
1188  fout << "$ & $";
1189 
1190  print.printFormat(fout, fitFracAsymm);
1191  fout << "$ & $";
1192 
1193  print.printFormat(fout, acp);
1194  fout << " \\pm ";
1195  print.printFormat(fout, acpErr);
1196  fout << "$ \\\\" << std::endl;
1197  }
1198  fout << "\\hline" << std::endl;
1199 
1200  // Also print out sum of fit fractions
1201  fout << "Fit Fraction Sum & $";
1202  print.printFormat(fout, negFitFracSum);
1203  fout << "$ & $";
1204  print.printFormat(fout, posFitFracSum);
1205  fout << "$ & & \\\\" << std::endl;
1206  fout << "\\hline" << std::endl;
1207 
1208  fout << "DP rate & $";
1209  print.printFormat(fout, negDPRate_.value());
1210  fout << "$ & $";
1211  print.printFormat(fout, posDPRate_.value());
1212  fout << "$ & & \\\\" << std::endl;
1213 
1214  fout << "$< \\varepsilon > $ & $";
1215  print.printFormat(fout, negMeanEff_.value());
1216  fout << "$ & $";
1217  print.printFormat(fout, posMeanEff_.value());
1218  fout << "$ & & \\\\" << std::endl;
1219  fout << "\\hline" << std::endl;
1220  fout << "\\end{tabular}" << std::endl << std::endl;
1221  }
1222 
1223  if (!negSignalPdfs_.empty()) {
1224  fout << "\\begin{tabular}{|l|c|}" << std::endl;
1225  fout << "\\hline" << std::endl;
1226  if (useSCF_ == kTRUE) {
1227  fout << "\\Extra TM Signal PDFs' Parameters: & \\\\" << std::endl;
1228  } else {
1229  fout << "\\Extra Signal PDFs' Parameters: & \\\\" << std::endl;
1230  }
1231  this->printFitParameters(negSignalPdfs_, fout);
1232  if ( tagged_ ) {
1233  this->printFitParameters(posSignalPdfs_, fout);
1234  }
1235  if (useSCF_ == kTRUE && !negScfPdfs_.empty()) {
1236  fout << "\\hline" << std::endl;
1237  fout << "\\Extra SCF Signal PDFs' Parameters: & \\\\" << std::endl;
1238  this->printFitParameters(negScfPdfs_, fout);
1239  if ( tagged_ ) {
1240  this->printFitParameters(posScfPdfs_, fout);
1241  }
1242  }
1243  if (usingBkgnd_ == kTRUE && !negBkgndPdfs_.empty()) {
1244  fout << "\\hline" << std::endl;
1245  fout << "\\Extra Background PDFs' Parameters: & \\\\" << std::endl;
1246  for (LauBkgndPdfsList::const_iterator iter = negBkgndPdfs_.begin(); iter != negBkgndPdfs_.end(); ++iter) {
1247  this->printFitParameters(*iter, fout);
1248  }
1249  if ( tagged_ ) {
1250  for (LauBkgndPdfsList::const_iterator iter = posBkgndPdfs_.begin(); iter != posBkgndPdfs_.end(); ++iter) {
1251  this->printFitParameters(*iter, fout);
1252  }
1253  }
1254  }
1255  fout << "\\hline \n\\end{tabular}" << std::endl << std::endl;
1256  }
1257 }
1258 
1260 {
1261  // Update the number of signal events to be total-sum(background events)
1262  this->updateSigEvents();
1263 
1264  // Check whether we want to have randomised initial fit parameters for the signal model
1265  if (this->useRandomInitFitPars() == kTRUE) {
1266  std::cout << "INFO in LauCPFitModel::checkInitFitParams : Setting random parameters for the signal model" << std::endl;
1267  this->randomiseInitFitPars();
1268  }
1269 }
1270 
1272 {
1273  // Only randomise those parameters that are not fixed!
1274  std::cout << "INFO in LauCPFitModel::randomiseInitFitPars : Randomising the initial fit magnitudes and phases of the components..." << std::endl;
1275 
1276  for (UInt_t i = 0; i < nSigComp_; i++) {
1277  coeffPars_[i]->randomiseInitValues();
1278  }
1279 }
1280 
1282 {
1283  // Determine the number of events to generate for each hypothesis
1284  // If we're smearing then smear each one individually
1285 
1286  LauGenInfo nEvtsGen;
1287 
1288  // Signal
1289  Double_t evtWeight(1.0);
1290  Double_t nEvts = signalEvents_->genValue();
1291  if ( nEvts < 0.0 ) {
1292  evtWeight = -1.0;
1293  nEvts = TMath::Abs( nEvts );
1294  }
1295  Double_t asym(0.0);
1296  Double_t sigAsym(0.0);
1297  // need to include this as an alternative in case the DP isn't in the model
1298  if ( !this->useDP() || forceAsym_ ) {
1299  sigAsym = signalAsym_->genValue();
1300  } else {
1301  Double_t negRate = negSigModel_->getDPNorm();
1302  Double_t posRate = posSigModel_->getDPNorm();
1303  if (negRate+posRate>1e-30) {
1304  sigAsym = (negRate-posRate)/(negRate+posRate);
1305  }
1306  }
1307  asym = sigAsym;
1308  Int_t nPosEvts = static_cast<Int_t>((nEvts/2.0 * (1.0 - asym)) + 0.5);
1309  Int_t nNegEvts = static_cast<Int_t>((nEvts/2.0 * (1.0 + asym)) + 0.5);
1310  if (this->doPoissonSmearing()) {
1311  nNegEvts = LauRandom::randomFun()->Poisson(nNegEvts);
1312  nPosEvts = LauRandom::randomFun()->Poisson(nPosEvts);
1313  }
1314  nEvtsGen[std::make_pair("signal",-1)] = std::make_pair(nNegEvts,evtWeight);
1315  nEvtsGen[std::make_pair("signal",+1)] = std::make_pair(nPosEvts,evtWeight);
1316 
1317  // backgrounds
1318  const UInt_t nBkgnds = this->nBkgndClasses();
1319  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
1320  const TString& bkgndClass = this->bkgndClassName(bkgndID);
1321  const LauParameter* evtsPar = bkgndEvents_[bkgndID];
1322  const LauParameter* asymPar = bkgndAsym_[bkgndID];
1323  evtWeight = 1.0;
1324  nEvts = TMath::FloorNint( evtsPar->genValue() );
1325  if ( nEvts < 0 ) {
1326  evtWeight = -1.0;
1327  nEvts = TMath::Abs( nEvts );
1328  }
1329  asym = asymPar->genValue();
1330  nPosEvts = static_cast<Int_t>((nEvts/2.0 * (1.0 - asym)) + 0.5);
1331  nNegEvts = static_cast<Int_t>((nEvts/2.0 * (1.0 + asym)) + 0.5);
1332  if (this->doPoissonSmearing()) {
1333  nNegEvts = LauRandom::randomFun()->Poisson(nNegEvts);
1334  nPosEvts = LauRandom::randomFun()->Poisson(nPosEvts);
1335  }
1336  nEvtsGen[std::make_pair(bkgndClass,-1)] = std::make_pair(nNegEvts,evtWeight);
1337  nEvtsGen[std::make_pair(bkgndClass,+1)] = std::make_pair(nPosEvts,evtWeight);
1338  }
1339 
1340  std::cout << "INFO in LauCPFitModel::eventsToGenerate : Generating toy MC with:" << std::endl;
1341  std::cout << " : Signal asymmetry = " << sigAsym << " and number of signal events = " << signalEvents_->genValue() << std::endl;
1342  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
1343  const TString& bkgndClass = this->bkgndClassName(bkgndID);
1344  const LauParameter* evtsPar = bkgndEvents_[bkgndID];
1345  const LauParameter* asymPar = bkgndAsym_[bkgndID];
1346  std::cout << " : " << bkgndClass << " asymmetry = " << asymPar->genValue() << " and number of " << bkgndClass << " events = " << evtsPar->genValue() << std::endl;
1347  }
1348 
1349  return nEvtsGen;
1350 }
1351 
1353 {
1354  // Routine to generate toy Monte Carlo events according to the various models we have defined.
1355 
1356  // Determine the number of events to generate for each hypothesis
1357  LauGenInfo nEvts = this->eventsToGenerate();
1358 
1359  Bool_t genOK(kTRUE);
1360  Int_t evtNum(0);
1361 
1362  const UInt_t nBkgnds = this->nBkgndClasses();
1363  std::vector<TString> bkgndClassNames(nBkgnds);
1364  std::vector<TString> bkgndClassNamesGen(nBkgnds);
1365  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
1366  TString name( this->bkgndClassName(iBkgnd) );
1367  bkgndClassNames[iBkgnd] = name;
1368  bkgndClassNamesGen[iBkgnd] = "gen"+name;
1369  }
1370 
1371  const Bool_t storeSCFTruthInfo = ( useSCF_ || ( this->enableEmbedding() &&
1372  negSignalTree_ != 0 && negSignalTree_->haveBranch("mcMatch") &&
1373  posSignalTree_ != 0 && posSignalTree_->haveBranch("mcMatch") ) );
1374 
1375  // Loop over the hypotheses and generate the requested number of events for each one
1376  for (LauGenInfo::const_iterator iter = nEvts.begin(); iter != nEvts.end(); ++iter) {
1377 
1378  const TString& type(iter->first.first);
1379  curEvtCharge_ = iter->first.second;
1380  Double_t evtWeight( iter->second.second );
1381  Int_t nEvtsGen( iter->second.first );
1382 
1383  for (Int_t iEvt(0); iEvt<nEvtsGen; ++iEvt) {
1384 
1385  this->setGenNtupleDoubleBranchValue( "evtWeight", evtWeight );
1386  this->setGenNtupleDoubleBranchValue( "efficiency", 1.0 );
1387 
1388  if (type == "signal") {
1389  this->setGenNtupleIntegerBranchValue("genSig",1);
1390  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
1391  this->setGenNtupleIntegerBranchValue( bkgndClassNamesGen[iBkgnd], 0 );
1392  }
1393  genOK = this->generateSignalEvent();
1394  if ( curEvtCharge_ > 0 ){
1395  this->setGenNtupleDoubleBranchValue( "efficiency", posSigModel_->getEvtEff() );
1396  } else {
1397  this->setGenNtupleDoubleBranchValue( "efficiency", negSigModel_->getEvtEff() );
1398  }
1399 
1400  } else {
1401  this->setGenNtupleIntegerBranchValue("genSig",0);
1402  if ( storeSCFTruthInfo ) {
1403  this->setGenNtupleIntegerBranchValue("genTMSig",0);
1404  this->setGenNtupleIntegerBranchValue("genSCFSig",0);
1405  }
1406  UInt_t bkgndID(0);
1407  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
1408  Int_t gen(0);
1409  if ( bkgndClassNames[iBkgnd] == type ) {
1410  gen = 1;
1411  bkgndID = iBkgnd;
1412  }
1413  this->setGenNtupleIntegerBranchValue( bkgndClassNamesGen[iBkgnd], gen );
1414  }
1415  genOK = this->generateBkgndEvent(bkgndID);
1416  }
1417  if (!genOK) {
1418  // If there was a problem with the generation then break out and return.
1419  // The problem model will have adjusted itself so that all should be OK next time.
1420  break;
1421  }
1422  if (this->useDP() == kTRUE) {
1423  this->setDPBranchValues();
1424  }
1425 
1426  // Store the event charge
1428 
1429  // Store the event number (within this experiment)
1430  // and then increment it
1431  this->setGenNtupleIntegerBranchValue("iEvtWithinExpt",evtNum);
1432  ++evtNum;
1433 
1434  this->fillGenNtupleBranches();
1435  if (iEvt%500 == 0) {std::cout << "INFO in LauCPFitModel::genExpt : Generated event number " << iEvt << " out of " << nEvtsGen << " " << type << " events." << std::endl;}
1436  }
1437 
1438  if (!genOK) {
1439  break;
1440  }
1441  }
1442 
1443  if (this->useDP() && genOK) {
1444 
1445  negSigModel_->checkToyMC(kTRUE,kTRUE);
1446  posSigModel_->checkToyMC(kTRUE,kTRUE);
1447 
1448  // Get the fit fractions if they're to be written into the latex table
1449  if (this->writeLatexTable()) {
1450  LauParArray negFitFrac = negSigModel_->getFitFractions();
1451  if (negFitFrac.size() != nSigComp_) {
1452  std::cerr << "ERROR in LauCPFitModel::genExpt : Fit Fraction array of unexpected dimension: " << negFitFrac.size() << std::endl;
1453  gSystem->Exit(EXIT_FAILURE);
1454  }
1455  for (UInt_t i(0); i<nSigComp_; ++i) {
1456  if (negFitFrac[i].size() != nSigComp_) {
1457  std::cerr << "ERROR in LauCPFitModel::genExpt : Fit Fraction array of unexpected dimension: " << negFitFrac[i].size() << std::endl;
1458  gSystem->Exit(EXIT_FAILURE);
1459  }
1460  }
1461  LauParArray posFitFrac = posSigModel_->getFitFractions();
1462  if (posFitFrac.size() != nSigComp_) {
1463  std::cerr << "ERROR in LauCPFitModel::genExpt : Fit Fraction array of unexpected dimension: " << posFitFrac.size() << std::endl;
1464  gSystem->Exit(EXIT_FAILURE);
1465  }
1466  for (UInt_t i(0); i<nSigComp_; ++i) {
1467  if (posFitFrac[i].size() != nSigComp_) {
1468  std::cerr << "ERROR in LauCPFitModel::genExpt : Fit Fraction array of unexpected dimension: " << posFitFrac[i].size() << std::endl;
1469  gSystem->Exit(EXIT_FAILURE);
1470  }
1471  }
1472 
1473  for (UInt_t i(0); i<nSigComp_; ++i) {
1474  for (UInt_t j(i); j<nSigComp_; ++j) {
1475  negFitFrac_[i][j].value(negFitFrac[i][j].value());
1476  posFitFrac_[i][j].value(posFitFrac[i][j].value());
1477  }
1478  }
1483  }
1484  }
1485 
1486  // If we're reusing embedded events or if the generation is being
1487  // reset then clear the lists of used events
1488  if (reuseSignal_ || !genOK) {
1489  if (negSignalTree_) {
1491  }
1492  if (posSignalTree_) {
1494  }
1495  }
1496  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
1497  LauEmbeddedData* data = negBkgndTree_[bkgndID];
1498  if (reuseBkgnd_[bkgndID] || !genOK) {
1499  if (data) {
1500  data->clearUsedList();
1501  }
1502  }
1503  }
1504  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
1505  LauEmbeddedData* data = posBkgndTree_[bkgndID];
1506  if (reuseBkgnd_[bkgndID] || !genOK) {
1507  if (data) {
1508  data->clearUsedList();
1509  }
1510  }
1511  }
1512 
1513  return genOK;
1514 }
1515 
1517 {
1518  // Generate signal event
1519  Bool_t genOK(kTRUE);
1520  Bool_t genSCF(kFALSE);
1521 
1522  LauAbsDPDynamics* model(0);
1523  LauKinematics* kinematics(0);
1524  LauEmbeddedData* embeddedData(0);
1525  LauPdfList* sigPdfs(0);
1526  LauPdfList* scfPdfs(0);
1527 
1528  Bool_t doReweighting(kFALSE);
1529 
1530  if (curEvtCharge_<0) {
1531  model = negSigModel_;
1532  kinematics = negKinematics_;
1533  sigPdfs = &negSignalPdfs_;
1534  scfPdfs = &negScfPdfs_;
1535  if (this->enableEmbedding()) {
1536  embeddedData = negSignalTree_;
1537  doReweighting = useNegReweighting_;
1538  }
1539  } else {
1540  model = posSigModel_;
1541  kinematics = posKinematics_;
1542  if ( tagged_ ) {
1543  sigPdfs = &posSignalPdfs_;
1544  scfPdfs = &posScfPdfs_;
1545  } else {
1546  sigPdfs = &negSignalPdfs_;
1547  scfPdfs = &negScfPdfs_;
1548  }
1549  if (this->enableEmbedding()) {
1550  embeddedData = posSignalTree_;
1551  doReweighting = usePosReweighting_;
1552  }
1553  }
1554 
1555  if (this->useDP()) {
1556  if (embeddedData) {
1557  if (doReweighting) {
1558  // Select a (random) event from the generated data. Then store the
1559  // reconstructed DP co-ords, together with other pdf information,
1560  // as the embedded data.
1561  genOK = embeddedData->getReweightedEvent(model);
1562  } else {
1563  // Just get the information of a (randomly) selected event in the
1564  // embedded data
1565  embeddedData->getEmbeddedEvent(kinematics);
1566  }
1567  genSCF = this->storeSignalMCMatch( embeddedData );
1568  } else {
1569  genOK = model->generate();
1570  if ( genOK && useSCF_ ) {
1571  Double_t frac(0.0);
1572  if ( useSCFHist_ ) {
1573  frac = scfFracHist_->calcEfficiency( kinematics );
1574  } else {
1575  frac = scfFrac_.genValue();
1576  }
1577  if ( frac < LauRandom::randomFun()->Rndm() ) {
1578  this->setGenNtupleIntegerBranchValue("genTMSig",1);
1579  this->setGenNtupleIntegerBranchValue("genSCFSig",0);
1580  genSCF = kFALSE;
1581  } else {
1582  this->setGenNtupleIntegerBranchValue("genTMSig",0);
1583  this->setGenNtupleIntegerBranchValue("genSCFSig",1);
1584  genSCF = kTRUE;
1585 
1586  // Optionally smear the DP position
1587  // of the SCF event
1588  if ( scfMap_ != 0 ) {
1589  Double_t xCoord(0.0), yCoord(0.0);
1590  if ( kinematics->squareDP() ) {
1591  xCoord = kinematics->getmPrime();
1592  yCoord = kinematics->getThetaPrime();
1593  } else {
1594  xCoord = kinematics->getm13Sq();
1595  yCoord = kinematics->getm23Sq();
1596  }
1597 
1598  // Find the bin number where this event is generated
1599  Int_t binNo = scfMap_->binNumber( xCoord, yCoord );
1600 
1601  // Retrieve the migration histogram
1602  TH2* histo = scfMap_->trueHist( binNo );
1603 
1604  LauAbsEffModel * effModel = model->getEffModel();
1605  do {
1606  // Get a random point from the histogram
1607  histo->GetRandom2( xCoord, yCoord );
1608 
1609  // Update the kinematics
1610  if ( kinematics->squareDP() ) {
1611  kinematics->updateSqDPKinematics( xCoord, yCoord );
1612  } else {
1613  kinematics->updateKinematics( xCoord, yCoord );
1614  }
1615  } while ( ! effModel->passVeto( kinematics ) );
1616  }
1617  }
1618  }
1619  }
1620  } else {
1621  if (embeddedData) {
1622  embeddedData->getEmbeddedEvent(0);
1623  genSCF = this->storeSignalMCMatch( embeddedData );
1624  } else if ( useSCF_ ) {
1625  Double_t frac = scfFrac_.genValue();
1626  if ( frac < LauRandom::randomFun()->Rndm() ) {
1627  this->setGenNtupleIntegerBranchValue("genTMSig",1);
1628  this->setGenNtupleIntegerBranchValue("genSCFSig",0);
1629  genSCF = kFALSE;
1630  } else {
1631  this->setGenNtupleIntegerBranchValue("genTMSig",0);
1632  this->setGenNtupleIntegerBranchValue("genSCFSig",1);
1633  genSCF = kTRUE;
1634  }
1635  }
1636  }
1637  if (genOK) {
1638  if ( useSCF_ ) {
1639  if ( genSCF ) {
1640  this->generateExtraPdfValues(scfPdfs, embeddedData);
1641  } else {
1642  this->generateExtraPdfValues(sigPdfs, embeddedData);
1643  }
1644  } else {
1645  this->generateExtraPdfValues(sigPdfs, embeddedData);
1646  }
1647  }
1648  // Check for problems with the embedding
1649  if (embeddedData && (embeddedData->nEvents() == embeddedData->nUsedEvents())) {
1650  std::cerr << "WARNING in LauCPFitModel::generateSignalEvent : Source of embedded signal events used up, clearing the list of used events." << std::endl;
1651  embeddedData->clearUsedList();
1652  }
1653 
1654  return genOK;
1655 }
1656 
1658 {
1659  // Generate Bkgnd event
1660  Bool_t genOK(kTRUE);
1661 
1662  LauAbsBkgndDPModel* model(0);
1663  LauEmbeddedData* embeddedData(0);
1664  LauPdfList* extraPdfs(0);
1665  LauKinematics* kinematics(0);
1666 
1667  if (curEvtCharge_<0) {
1668  model = negBkgndDPModels_[bkgndID];
1669  if (this->enableEmbedding()) {
1670  embeddedData = negBkgndTree_[bkgndID];
1671  }
1672  extraPdfs = &negBkgndPdfs_[bkgndID];
1673  kinematics = negKinematics_;
1674  } else {
1675  model = posBkgndDPModels_[bkgndID];
1676  if (this->enableEmbedding()) {
1677  embeddedData = posBkgndTree_[bkgndID];
1678  }
1679  if ( tagged_ ) {
1680  extraPdfs = &posBkgndPdfs_[bkgndID];
1681  } else {
1682  extraPdfs = &negBkgndPdfs_[bkgndID];
1683  }
1684  kinematics = posKinematics_;
1685  }
1686 
1687  if (this->useDP()) {
1688  if (embeddedData) {
1689  embeddedData->getEmbeddedEvent(kinematics);
1690  } else {
1691  if (model == 0) {
1692  const TString& bkgndClass = this->bkgndClassName(bkgndID);
1693  std::cerr << "ERROR in LauCPFitModel::generateBkgndEvent : Can't find the DP model for background class \"" << bkgndClass << "\"." << std::endl;
1694  gSystem->Exit(EXIT_FAILURE);
1695  }
1696  genOK = model->generate();
1697  }
1698  } else {
1699  if (embeddedData) {
1700  embeddedData->getEmbeddedEvent(0);
1701  }
1702  }
1703  if (genOK) {
1704  this->generateExtraPdfValues(extraPdfs, embeddedData);
1705  }
1706 
1707  // Check for problems with the embedding
1708  if (embeddedData && (embeddedData->nEvents() == embeddedData->nUsedEvents())) {
1709  const TString& bkgndClass = this->bkgndClassName(bkgndID);
1710  std::cerr << "WARNING in LauCPFitModel::generateBkgndEvent : Source of embedded " << bkgndClass << " events used up, clearing the list of used events." << std::endl;
1711  embeddedData->clearUsedList();
1712  }
1713 
1714  return genOK;
1715 }
1716 
1718 {
1719  // Setup the required ntuple branches
1720  this->addGenNtupleDoubleBranch("evtWeight");
1721  this->addGenNtupleIntegerBranch("genSig");
1722  this->addGenNtupleDoubleBranch("efficiency");
1723  if ( useSCF_ || ( this->enableEmbedding() &&
1724  negSignalTree_ != 0 && negSignalTree_->haveBranch("mcMatch") &&
1725  posSignalTree_ != 0 && posSignalTree_->haveBranch("mcMatch") ) ) {
1726  this->addGenNtupleIntegerBranch("genTMSig");
1727  this->addGenNtupleIntegerBranch("genSCFSig");
1728  }
1729  const UInt_t nBkgnds = this->nBkgndClasses();
1730  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
1731  TString name( this->bkgndClassName(iBkgnd) );
1732  name.Prepend("gen");
1733  this->addGenNtupleIntegerBranch(name);
1734  }
1735  this->addGenNtupleIntegerBranch("charge");
1736  if (this->useDP() == kTRUE) {
1737  this->addGenNtupleDoubleBranch("m12");
1738  this->addGenNtupleDoubleBranch("m23");
1739  this->addGenNtupleDoubleBranch("m13");
1740  this->addGenNtupleDoubleBranch("m12Sq");
1741  this->addGenNtupleDoubleBranch("m23Sq");
1742  this->addGenNtupleDoubleBranch("m13Sq");
1743  this->addGenNtupleDoubleBranch("cosHel12");
1744  this->addGenNtupleDoubleBranch("cosHel23");
1745  this->addGenNtupleDoubleBranch("cosHel13");
1747  this->addGenNtupleDoubleBranch("mPrime");
1748  this->addGenNtupleDoubleBranch("thPrime");
1749  }
1750  }
1751  for (LauPdfList::const_iterator pdf_iter = negSignalPdfs_.begin(); pdf_iter != negSignalPdfs_.end(); ++pdf_iter) {
1752  std::vector<TString> varNames = (*pdf_iter)->varNames();
1753  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
1754  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
1755  this->addGenNtupleDoubleBranch( (*var_iter) );
1756  }
1757  }
1758  }
1759 }
1760 
1762 {
1763  LauKinematics* kinematics(0);
1764  if (curEvtCharge_<0) {
1765  kinematics = negKinematics_;
1766  } else {
1767  kinematics = posKinematics_;
1768  }
1769 
1770  // Store all the DP information
1771  this->setGenNtupleDoubleBranchValue("m12", kinematics->getm12());
1772  this->setGenNtupleDoubleBranchValue("m23", kinematics->getm23());
1773  this->setGenNtupleDoubleBranchValue("m13", kinematics->getm13());
1774  this->setGenNtupleDoubleBranchValue("m12Sq", kinematics->getm12Sq());
1775  this->setGenNtupleDoubleBranchValue("m23Sq", kinematics->getm23Sq());
1776  this->setGenNtupleDoubleBranchValue("m13Sq", kinematics->getm13Sq());
1777  this->setGenNtupleDoubleBranchValue("cosHel12", kinematics->getc12());
1778  this->setGenNtupleDoubleBranchValue("cosHel23", kinematics->getc23());
1779  this->setGenNtupleDoubleBranchValue("cosHel13", kinematics->getc13());
1780  if (kinematics->squareDP()) {
1781  this->setGenNtupleDoubleBranchValue("mPrime", kinematics->getmPrime());
1782  this->setGenNtupleDoubleBranchValue("thPrime", kinematics->getThetaPrime());
1783  }
1784 }
1785 
1787 {
1788  LauKinematics* kinematics(0);
1789  if (curEvtCharge_<0) {
1790  kinematics = negKinematics_;
1791  } else {
1792  kinematics = posKinematics_;
1793  }
1794 
1795  if (!extraPdfs) {
1796  std::cerr << "ERROR in LauCPFitModel::generateExtraPdfValues : Null pointer to PDF list." << std::endl;
1797  gSystem->Exit(EXIT_FAILURE);
1798  }
1799 
1800  if (extraPdfs->empty()) {
1801  //std::cerr << "WARNING in LauCPFitModel::generateExtraPdfValues : PDF list is empty." << std::endl;
1802  return;
1803  }
1804 
1805  // Generate from the extra PDFs
1806  for (LauPdfList::iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end(); ++pdf_iter) {
1807  LauFitData genValues;
1808  if (embeddedData) {
1809  genValues = embeddedData->getValues( (*pdf_iter)->varNames() );
1810  } else {
1811  genValues = (*pdf_iter)->generate(kinematics);
1812  }
1813  for ( LauFitData::const_iterator var_iter = genValues.begin(); var_iter != genValues.end(); ++var_iter ) {
1814  TString varName = var_iter->first;
1815  if ( varName != "m13Sq" && varName != "m23Sq" ) {
1816  Double_t value = var_iter->second;
1817  this->setGenNtupleDoubleBranchValue(varName,value);
1818  }
1819  }
1820  }
1821 }
1822 
1824 {
1825  // Default to TM
1826  Bool_t genSCF(kFALSE);
1827  Int_t match(1);
1828 
1829  // Check that we have a valid pointer and that embedded data has
1830  // the mcMatch branch. If so then get the match value.
1831  if ( embeddedData && embeddedData->haveBranch("mcMatch") ) {
1832  match = TMath::Nint( embeddedData->getValue("mcMatch") );
1833  }
1834 
1835  // Set the variables accordingly.
1836  if (match) {
1837  this->setGenNtupleIntegerBranchValue("genTMSig",1);
1838  this->setGenNtupleIntegerBranchValue("genSCFSig",0);
1839  genSCF = kFALSE;
1840  } else {
1841  this->setGenNtupleIntegerBranchValue("genTMSig",0);
1842  this->setGenNtupleIntegerBranchValue("genSCFSig",1);
1843  genSCF = kTRUE;
1844  }
1845 
1846  return genSCF;
1847 }
1848 
1850 {
1851  // Update the signal parameters and then the total normalisation for the signal likelihood
1852  if (this->useDP() == kTRUE) {
1853  this->updateCoeffs();
1856  }
1857 
1858  // Update the signal fraction from the background fractions if not doing an extended fit
1859  if ( !this->doEMLFit() && !signalEvents_->fixed() ) {
1860  this->updateSigEvents();
1861  }
1862 }
1863 
1865 {
1866  // The background parameters will have been set from Minuit.
1867  // We need to update the signal events using these.
1868  Double_t nTotEvts = this->eventsPerExpt();
1869 
1870  signalEvents_->range(-2.0*nTotEvts,2.0*nTotEvts);
1871  for (LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
1872  (*iter)->range(-2.0*nTotEvts,2.0*nTotEvts);
1873  }
1874 
1875  if (signalEvents_->fixed()) {
1876  return;
1877  }
1878 
1879  // Subtract background events (if any) from signal.
1880  Double_t signalEvents = nTotEvts;
1881  if (usingBkgnd_ == kTRUE) {
1882  for (LauBkgndYieldList::const_iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
1883  signalEvents -= (*iter)->value();
1884  }
1885  }
1886  signalEvents_->value(signalEvents);
1887 }
1888 
1890 {
1891  // Fill the internal data trees of the signal and background models.
1892  // Note that we store the events of both charges in both the
1893  // negative and the positive models. It's only later, at the stage
1894  // when the likelihood is being calculated, that we separate them.
1895 
1896  LauFitDataTree* inputFitData = this->fitData();
1897 
1898  // First the Dalitz plot variables (m_ij^2)
1899  if (this->useDP() == kTRUE) {
1900 
1901  // need to append SCF smearing bins before caching DP amplitudes
1902  if ( scfMap_ != 0 ) {
1903  this->appendBinCentres( inputFitData );
1904  }
1905 
1906  negSigModel_->fillDataTree(*inputFitData);
1907  posSigModel_->fillDataTree(*inputFitData);
1908 
1909  if (usingBkgnd_ == kTRUE) {
1910  for (LauBkgndDPModelList::iterator iter = negBkgndDPModels_.begin(); iter != negBkgndDPModels_.end(); ++iter) {
1911  (*iter)->fillDataTree(*inputFitData);
1912  }
1913  for (LauBkgndDPModelList::iterator iter = posBkgndDPModels_.begin(); iter != posBkgndDPModels_.end(); ++iter) {
1914  (*iter)->fillDataTree(*inputFitData);
1915  }
1916  }
1917  }
1918 
1919  // ...and then the extra PDFs
1920  this->cacheInfo(negSignalPdfs_, *inputFitData);
1921  this->cacheInfo(negScfPdfs_, *inputFitData);
1922  for (LauBkgndPdfsList::iterator iter = negBkgndPdfs_.begin(); iter != negBkgndPdfs_.end(); ++iter) {
1923  this->cacheInfo((*iter), *inputFitData);
1924  }
1925  if ( tagged_ ) {
1926  this->cacheInfo(posSignalPdfs_, *inputFitData);
1927  this->cacheInfo(posScfPdfs_, *inputFitData);
1928  for (LauBkgndPdfsList::iterator iter = posBkgndPdfs_.begin(); iter != posBkgndPdfs_.end(); ++iter) {
1929  this->cacheInfo((*iter), *inputFitData);
1930  }
1931  }
1932 
1933  // the SCF fractions and jacobians
1934  if ( useSCF_ && useSCFHist_ ) {
1935  if ( !inputFitData->haveBranch( "m13Sq" ) || !inputFitData->haveBranch( "m23Sq" ) ) {
1936  std::cerr << "ERROR in LauCPFitModel::cacheInputFitVars : Input data does not contain DP branches and so can't cache the SCF fraction." << std::endl;
1937  gSystem->Exit(EXIT_FAILURE);
1938  }
1939 
1940  UInt_t nEvents = inputFitData->nEvents();
1941  recoSCFFracs_.clear();
1942  recoSCFFracs_.reserve( nEvents );
1943  if ( negKinematics_->squareDP() ) {
1944  recoJacobians_.clear();
1945  recoJacobians_.reserve( nEvents );
1946  }
1947  for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) {
1948  const LauFitData& dataValues = inputFitData->getData(iEvt);
1949  LauFitData::const_iterator m13_iter = dataValues.find("m13Sq");
1950  LauFitData::const_iterator m23_iter = dataValues.find("m23Sq");
1951  negKinematics_->updateKinematics( m13_iter->second, m23_iter->second );
1952  Double_t scfFrac = scfFracHist_->calcEfficiency( negKinematics_ );
1953  recoSCFFracs_.push_back( scfFrac );
1954  if ( negKinematics_->squareDP() ) {
1956  }
1957  }
1958  }
1959 
1960  // finally cache the event charge
1961  evtCharges_.clear();
1962  if ( tagged_ ) {
1963  if ( !inputFitData->haveBranch( tagVarName_ ) ) {
1964  std::cerr << "ERROR in LauCPFitModel::cacheInputFitVars : Input data does not contain branch \"" << tagVarName_ << "\"." << std::endl;
1965  gSystem->Exit(EXIT_FAILURE);
1966  }
1967  UInt_t nEvents = inputFitData->nEvents();
1968  evtCharges_.reserve( nEvents );
1969  for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) {
1970  const LauFitData& dataValues = inputFitData->getData(iEvt);
1971  LauFitData::const_iterator iter = dataValues.find( tagVarName_ );
1972  curEvtCharge_ = static_cast<Int_t>( iter->second );
1973  evtCharges_.push_back( curEvtCharge_ );
1974  }
1975  }
1976 }
1977 
1979 {
1980  // We'll be caching the DP amplitudes and efficiencies of the centres of the true bins.
1981  // To do so, we attach some fake points at the end of inputData, the number of the entry
1982  // minus the total number of events corresponding to the number of the histogram for that
1983  // given true bin in the LauScfMap object. (What this means is that when Laura is provided with
1984  // the LauScfMap object by the user, it's the latter who has to make sure that it contains the
1985  // right number of histograms and in exactly the right order!)
1986 
1987  // Get the x and y co-ordinates of the bin centres
1988  std::vector<Double_t> binCentresXCoords;
1989  std::vector<Double_t> binCentresYCoords;
1990  scfMap_->listBinCentres(binCentresXCoords, binCentresYCoords);
1991 
1992  // The SCF histograms could be in square Dalitz plot histograms.
1993  // The dynamics takes normal Dalitz plot coords, so we might have to convert them back.
1994  Bool_t sqDP = negKinematics_->squareDP();
1995  UInt_t nBins = binCentresXCoords.size();
1996  fakeSCFFracs_.clear();
1997  fakeSCFFracs_.reserve( nBins );
1998  if ( sqDP ) {
1999  fakeJacobians_.clear();
2000  fakeJacobians_.reserve( nBins );
2001  }
2002 
2003  for (UInt_t iBin = 0; iBin < nBins; ++iBin) {
2004 
2005  if ( sqDP ) {
2006  negKinematics_->updateSqDPKinematics(binCentresXCoords[iBin],binCentresYCoords[iBin]);
2007  binCentresXCoords[iBin] = negKinematics_->getm13Sq();
2008  binCentresYCoords[iBin] = negKinematics_->getm23Sq();
2010  } else {
2011  negKinematics_->updateKinematics(binCentresXCoords[iBin],binCentresYCoords[iBin]);
2012  }
2013 
2015  }
2016 
2017  // Set up inputFitVars_ object to hold the fake events
2018  inputData->appendFakePoints(binCentresXCoords,binCentresYCoords);
2019 }
2020 
2022 {
2023  // Find out whether we have B- or B+
2024  if ( tagged_ ) {
2025  curEvtCharge_ = evtCharges_[iEvt];
2026 
2027  // check that the charge is either +1 or -1
2028  if (TMath::Abs(curEvtCharge_)!=1) {
2029  std::cerr << "ERROR in LauCPFitModel::getTotEvtLikelihood : Charge/tag not accepted value: " << curEvtCharge_ << std::endl;
2030  if (curEvtCharge_>0) {
2031  curEvtCharge_ = +1;
2032  } else {
2033  curEvtCharge_ = -1;
2034  }
2035  std::cerr << " : Making it: " << curEvtCharge_ << "." << std::endl;
2036  }
2037  }
2038 
2039  // Get the DP likelihood for signal and backgrounds
2040  this->getEvtDPLikelihood(iEvt);
2041 
2042  // Get the combined extra PDFs likelihood for signal and backgrounds
2043  this->getEvtExtraLikelihoods(iEvt);
2044 
2045  // If appropriate, combine the TM and SCF likelihoods
2046  Double_t sigLike = sigDPLike_ * sigExtraLike_;
2047  if ( useSCF_ ) {
2048  Double_t scfFrac(0.0);
2049  if (useSCFHist_) {
2050  scfFrac = recoSCFFracs_[iEvt];
2051  } else {
2052  scfFrac = scfFrac_.value();
2053  }
2054  sigLike *= (1.0 - scfFrac);
2055  if ( (scfMap_ != 0) && (this->useDP() == kTRUE) ) {
2056  // if we're smearing the SCF DP PDF then the SCF frac
2057  // is already included in the SCF DP likelihood
2058  sigLike += (scfDPLike_ * scfExtraLike_);
2059  } else {
2060  sigLike += (scfFrac * scfDPLike_ * scfExtraLike_);
2061  }
2062  }
2063 
2064  // Get the correct event fractions depending on the charge
2065  // Signal asymmetry is built into the DP model... but when the DP
2066  // isn't in the fit we need an explicit parameter
2067  Double_t signalEvents = signalEvents_->value() * 0.5;
2068  if (this->useDP() == kFALSE) {
2069  signalEvents *= (1.0 - curEvtCharge_ * signalAsym_->value());
2070  }
2071 
2072  // Construct the total event likelihood
2073  Double_t likelihood(0.0);
2074  if (usingBkgnd_) {
2075  likelihood = sigLike*signalEvents;
2076  const UInt_t nBkgnds = this->nBkgndClasses();
2077  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
2078  Double_t bkgndEvents = bkgndEvents_[bkgndID]->value() * 0.5 * (1.0 - curEvtCharge_ * bkgndAsym_[bkgndID]->value());
2079  likelihood += bkgndEvents*bkgndDPLike_[bkgndID]*bkgndExtraLike_[bkgndID];
2080  }
2081  } else {
2082  likelihood = sigLike*0.5;
2083  }
2084  return likelihood;
2085 }
2086 
2088 {
2089  Double_t eventSum(0.0);
2090  eventSum += signalEvents_->value();
2091  if (usingBkgnd_) {
2092  for (LauBkgndYieldList::const_iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
2093  eventSum += (*iter)->value();
2094  }
2095  }
2096  return eventSum;
2097 }
2098 
2100 {
2101  // Function to return the signal and background likelihoods for the
2102  // Dalitz plot for the given event evtNo.
2103 
2104  if ( ! this->useDP() ) {
2105  // There's always going to be a term in the likelihood for the
2106  // signal, so we'd better not zero it.
2107  sigDPLike_ = 1.0;
2108  scfDPLike_ = 1.0;
2109 
2110  const UInt_t nBkgnds = this->nBkgndClasses();
2111  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
2112  if (usingBkgnd_ == kTRUE) {
2113  bkgndDPLike_[bkgndID] = 1.0;
2114  } else {
2115  bkgndDPLike_[bkgndID] = 0.0;
2116  }
2117  }
2118 
2119  return;
2120  }
2121 
2122  const UInt_t nBkgnds = this->nBkgndClasses();
2123  if ( tagged_ ) {
2124  if (curEvtCharge_==+1) {
2128 
2129  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
2130  if (usingBkgnd_ == kTRUE) {
2131  bkgndDPLike_[bkgndID] = posBkgndDPModels_[bkgndID]->getLikelihood(iEvt);
2132  } else {
2133  bkgndDPLike_[bkgndID] = 0.0;
2134  }
2135  }
2136  } else {
2140 
2141  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
2142  if (usingBkgnd_ == kTRUE) {
2143  bkgndDPLike_[bkgndID] = negBkgndDPModels_[bkgndID]->getLikelihood(iEvt);
2144  } else {
2145  bkgndDPLike_[bkgndID] = 0.0;
2146  }
2147  }
2148  }
2149  } else {
2152 
2155 
2156  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
2157  if (usingBkgnd_ == kTRUE) {
2158  bkgndDPLike_[bkgndID] = 0.5 * ( posBkgndDPModels_[bkgndID]->getLikelihood(iEvt) +
2159  negBkgndDPModels_[bkgndID]->getLikelihood(iEvt) );
2160  } else {
2161  bkgndDPLike_[bkgndID] = 0.0;
2162  }
2163  }
2164  }
2165 
2166  if ( useSCF_ == kTRUE ) {
2167  if ( scfMap_ == 0 ) {
2168  // we're not smearing the SCF DP position
2169  // so the likelihood is the same as the TM
2171  } else {
2172  // calculate the smeared SCF DP likelihood
2173  scfDPLike_ = this->getEvtSCFDPLikelihood(iEvt);
2174  }
2175  }
2176 
2177  // Calculate the signal normalisation
2178  // NB the 2.0 is there so that the 0.5 factor is applied to
2179  // signal and background in the same place otherwise you get
2180  // normalisation problems when you switch off the DP in the fit
2181  Double_t norm = negSigModel_->getDPNorm() + posSigModel_->getDPNorm();
2182  sigDPLike_ *= 2.0/norm;
2183  scfDPLike_ *= 2.0/norm;
2184 }
2185 
2187 {
2188  Double_t scfDPLike(0.0);
2189 
2190  Double_t recoJacobian(1.0);
2191  Double_t xCoord(0.0);
2192  Double_t yCoord(0.0);
2193 
2194  Bool_t squareDP = negKinematics_->squareDP();
2195  if ( squareDP ) {
2196  xCoord = negSigModel_->getEvtmPrime();
2197  yCoord = negSigModel_->getEvtthPrime();
2198  recoJacobian = recoJacobians_[iEvt];
2199  } else {
2200  xCoord = negSigModel_->getEvtm13Sq();
2201  yCoord = negSigModel_->getEvtm23Sq();
2202  }
2203 
2204  // Find the bin that our reco event falls in
2205  Int_t recoBin = scfMap_->binNumber( xCoord, yCoord );
2206 
2207  // Find out which true Bins contribute to the given reco bin
2208  const std::vector<Int_t>* trueBins = scfMap_->trueBins(recoBin);
2209 
2210  const Int_t nDataEvents = this->eventsPerExpt();
2211 
2212  // Loop over the true bins
2213  for (std::vector<Int_t>::const_iterator iter = trueBins->begin(); iter != trueBins->end(); ++iter)
2214  {
2215  Int_t trueBin = (*iter);
2216 
2217  // prob of a true event in the given true bin migrating to the reco bin
2218  Double_t pRecoGivenTrue = scfMap_->prob( recoBin, trueBin );
2219  Double_t pTrue(0.0);
2220 
2221  // We've cached the DP amplitudes and the efficiency for the
2222  // true bin centres, just after the data points
2223  if ( tagged_ ) {
2224  LauAbsDPDynamics* sigModel(0);
2225  if (curEvtCharge_<0) {
2226  sigModel = negSigModel_;
2227  } else {
2228  sigModel = posSigModel_;
2229  }
2230 
2231  sigModel->calcLikelihoodInfo( nDataEvents + trueBin );
2232 
2233  pTrue = sigModel->getEvtDPAmp().abs2() * sigModel->getEvtEff();
2234  } else {
2235  posSigModel_->calcLikelihoodInfo( nDataEvents + trueBin );
2236  negSigModel_->calcLikelihoodInfo( nDataEvents + trueBin );
2237 
2238  pTrue = 0.5 * ( posSigModel_->getEvtDPAmp().abs2() * posSigModel_->getEvtEff() +
2240  }
2241 
2242  // Get the cached SCF fraction (and jacobian if we're using the square DP)
2243  Double_t scfFraction = fakeSCFFracs_[ trueBin ];
2244  Double_t jacobian(1.0);
2245  if ( squareDP ) {
2246  jacobian = fakeJacobians_[ trueBin ];
2247  }
2248 
2249  scfDPLike += pTrue * jacobian * scfFraction * pRecoGivenTrue;
2250  }
2251 
2252  // Divide by the reco jacobian
2253  scfDPLike /= recoJacobian;
2254 
2255  return scfDPLike;
2256 }
2257 
2259 {
2260  // Function to return the signal and background likelihoods for the
2261  // extra variables for the given event evtNo.
2262 
2263  sigExtraLike_ = 1.0;
2264 
2265  const UInt_t nBkgnds = this->nBkgndClasses();
2266 
2267  if ( ! tagged_ || curEvtCharge_ < 0 ) {
2268  sigExtraLike_ = this->prodPdfValue( negSignalPdfs_, iEvt );
2269 
2270  if (useSCF_) {
2271  scfExtraLike_ = this->prodPdfValue( negScfPdfs_, iEvt );
2272  }
2273 
2274  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
2275  if (usingBkgnd_) {
2276  bkgndExtraLike_[bkgndID] = this->prodPdfValue( negBkgndPdfs_[bkgndID], iEvt );
2277  } else {
2278  bkgndExtraLike_[bkgndID] = 0.0;
2279  }
2280  }
2281  } else {
2282  sigExtraLike_ = this->prodPdfValue( posSignalPdfs_, iEvt );
2283 
2284  if (useSCF_) {
2285  scfExtraLike_ = this->prodPdfValue( posScfPdfs_, iEvt );
2286  }
2287 
2288  for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
2289  if (usingBkgnd_) {
2290  bkgndExtraLike_[bkgndID] = this->prodPdfValue( posBkgndPdfs_[bkgndID], iEvt );
2291  } else {
2292  bkgndExtraLike_[bkgndID] = 0.0;
2293  }
2294  }
2295  }
2296 }
2297 
2299 {
2300  negCoeffs_.clear(); posCoeffs_.clear();
2301  negCoeffs_.reserve(nSigComp_); posCoeffs_.reserve(nSigComp_);
2302  for (UInt_t i = 0; i < nSigComp_; i++) {
2303  negCoeffs_.push_back(coeffPars_[i]->antiparticleCoeff());
2304  posCoeffs_.push_back(coeffPars_[i]->particleCoeff());
2305  }
2306 }
2307 
2309 {
2310  // add branches for storing the experiment number and the number of
2311  // the event within the current experiment
2312  this->addSPlotNtupleIntegerBranch("iExpt");
2313  this->addSPlotNtupleIntegerBranch("iEvtWithinExpt");
2314 
2315  // Store the efficiency of the event (for inclusive BF calculations).
2316  if (this->storeDPEff()) {
2317  this->addSPlotNtupleDoubleBranch("efficiency");
2319  this->addSPlotNtupleDoubleBranch("scffraction");
2320  }
2321  }
2322 
2323  // Store the total event likelihood for each species.
2324  if (useSCF_) {
2325  this->addSPlotNtupleDoubleBranch("sigTMTotalLike");
2326  this->addSPlotNtupleDoubleBranch("sigSCFTotalLike");
2327  this->addSPlotNtupleDoubleBranch("sigSCFFrac");
2328  } else {
2329  this->addSPlotNtupleDoubleBranch("sigTotalLike");
2330  }
2331  if (usingBkgnd_) {
2332  const UInt_t nBkgnds = this->nBkgndClasses();
2333  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
2334  TString name( this->bkgndClassName(iBkgnd) );
2335  name += "TotalLike";
2336  this->addSPlotNtupleDoubleBranch(name);
2337  }
2338  }
2339 
2340  // Store the DP likelihoods
2341  if (this->useDP()) {
2342  if (useSCF_) {
2343  this->addSPlotNtupleDoubleBranch("sigTMDPLike");
2344  this->addSPlotNtupleDoubleBranch("sigSCFDPLike");
2345  } else {
2346  this->addSPlotNtupleDoubleBranch("sigDPLike");
2347  }
2348  if (usingBkgnd_) {
2349  const UInt_t nBkgnds = this->nBkgndClasses();
2350  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
2351  TString name( this->bkgndClassName(iBkgnd) );
2352  name += "DPLike";
2353  this->addSPlotNtupleDoubleBranch(name);
2354  }
2355  }
2356  }
2357 
2358  // Store the likelihoods for each extra PDF
2359  if (useSCF_) {
2360  this->addSPlotNtupleBranches(&negSignalPdfs_, "sigTM");
2361  this->addSPlotNtupleBranches(&negScfPdfs_, "sigSCF");
2362  } else {
2363  this->addSPlotNtupleBranches(&negSignalPdfs_, "sig");
2364  }
2365  if (usingBkgnd_) {
2366  const UInt_t nBkgnds = this->nBkgndClasses();
2367  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
2368  const TString& bkgndClass = this->bkgndClassName(iBkgnd);
2369  const LauPdfList* pdfList = &(negBkgndPdfs_[iBkgnd]);
2370  this->addSPlotNtupleBranches(pdfList, bkgndClass);
2371  }
2372  }
2373 }
2374 
2375 void LauCPFitModel::addSPlotNtupleBranches(const LauPdfList* extraPdfs, const TString& prefix)
2376 {
2377  if (extraPdfs) {
2378  // Loop through each of the PDFs
2379  for (LauPdfList::const_iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end(); ++pdf_iter) {
2380 
2381  // Count the number of input variables that are not
2382  // DP variables (used in the case where there is DP
2383  // dependence for e.g. MVA)
2384  UInt_t nVars(0);
2385  std::vector<TString> varNames = (*pdf_iter)->varNames();
2386  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
2387  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
2388  ++nVars;
2389  }
2390  }
2391 
2392  if ( nVars == 1 ) {
2393  // If the PDF only has one variable then
2394  // simply add one branch for that variable
2395  TString varName = (*pdf_iter)->varName();
2396  TString name(prefix);
2397  name += varName;
2398  name += "Like";
2399  this->addSPlotNtupleDoubleBranch(name);
2400  } else if ( nVars == 2 ) {
2401  // If the PDF has two variables then we
2402  // need a branch for them both together and
2403  // branches for each
2404  TString allVars("");
2405  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
2406  allVars += (*var_iter);
2407  TString name(prefix);
2408  name += (*var_iter);
2409  name += "Like";
2410  this->addSPlotNtupleDoubleBranch(name);
2411  }
2412  TString name(prefix);
2413  name += allVars;
2414  name += "Like";
2415  this->addSPlotNtupleDoubleBranch(name);
2416  } else {
2417  std::cerr << "WARNING in LauCPFitModel::addSPlotNtupleBranches : Can't yet deal with 3D PDFs." << std::endl;
2418  }
2419  }
2420  }
2421 }
2422 
2423 Double_t LauCPFitModel::setSPlotNtupleBranchValues(LauPdfList* extraPdfs, const TString& prefix, UInt_t iEvt)
2424 {
2425  // Store the PDF value for each variable in the list
2426  Double_t totalLike(1.0);
2427  Double_t extraLike(0.0);
2428  if (extraPdfs) {
2429  for (LauPdfList::iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end(); ++pdf_iter) {
2430 
2431  // calculate the likelihood for this event
2432  (*pdf_iter)->calcLikelihoodInfo(iEvt);
2433  extraLike = (*pdf_iter)->getLikelihood();
2434  totalLike *= extraLike;
2435 
2436  // Count the number of input variables that are not
2437  // DP variables (used in the case where there is DP
2438  // dependence for e.g. MVA)
2439  UInt_t nVars(0);
2440  std::vector<TString> varNames = (*pdf_iter)->varNames();
2441  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
2442  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
2443  ++nVars;
2444  }
2445  }
2446 
2447  if ( nVars == 1 ) {
2448  // If the PDF only has one variable then
2449  // simply store the value for that variable
2450  TString varName = (*pdf_iter)->varName();
2451  TString name(prefix);
2452  name += varName;
2453  name += "Like";
2454  this->setSPlotNtupleDoubleBranchValue(name, extraLike);
2455  } else if ( nVars == 2 ) {
2456  // If the PDF has two variables then we
2457  // store the value for them both together
2458  // and for each on their own
2459  TString allVars("");
2460  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
2461  allVars += (*var_iter);
2462  TString name(prefix);
2463  name += (*var_iter);
2464  name += "Like";
2465  Double_t indivLike = (*pdf_iter)->getLikelihood( (*var_iter) );
2466  this->setSPlotNtupleDoubleBranchValue(name, indivLike);
2467  }
2468  TString name(prefix);
2469  name += allVars;
2470  name += "Like";
2471  this->setSPlotNtupleDoubleBranchValue(name, extraLike);
2472  } else {
2473  std::cerr << "WARNING in LauCPFitModel::setSPlotNtupleBranchValues : Can't yet deal with 3D PDFs." << std::endl;
2474  }
2475  }
2476  }
2477  return totalLike;
2478 }
2479 
2481 {
2482  LauSPlot::NameSet nameSet;
2483  if (this->useDP()) {
2484  nameSet.insert("DP");
2485  }
2486  // Loop through all the signal PDFs
2487  for (LauPdfList::const_iterator pdf_iter = negSignalPdfs_.begin(); pdf_iter != negSignalPdfs_.end(); ++pdf_iter) {
2488  // Loop over the variables involved in each PDF
2489  std::vector<TString> varNames = (*pdf_iter)->varNames();
2490  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
2491  // If they are not DP coordinates then add them
2492  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
2493  nameSet.insert( (*var_iter) );
2494  }
2495  }
2496  }
2497  return nameSet;
2498 }
2499 
2501 {
2502  LauSPlot::NumbMap numbMap;
2503  if (!signalEvents_->fixed() && this->doEMLFit()) {
2504  numbMap["sig"] = signalEvents_->genValue();
2505  }
2506  if ( usingBkgnd_ ) {
2507  const UInt_t nBkgnds = this->nBkgndClasses();
2508  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
2509  const TString& bkgndClass = this->bkgndClassName(iBkgnd);
2510  const LauParameter* par = bkgndEvents_[iBkgnd];
2511  if (!par->fixed()) {
2512  numbMap[bkgndClass] = par->genValue();
2513  }
2514  }
2515  }
2516  return numbMap;
2517 }
2518 
2520 {
2521  LauSPlot::NumbMap numbMap;
2522  if (signalEvents_->fixed() && this->doEMLFit()) {
2523  numbMap["sig"] = signalEvents_->genValue();
2524  }
2525  if ( usingBkgnd_ ) {
2526  const UInt_t nBkgnds = this->nBkgndClasses();
2527  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
2528  const TString& bkgndClass = this->bkgndClassName(iBkgnd);
2529  const LauParameter* par = bkgndEvents_[iBkgnd];
2530  if (par->fixed()) {
2531  numbMap[bkgndClass] = par->genValue();
2532  }
2533  }
2534  }
2535  return numbMap;
2536 }
2537 
2539 {
2540  // This makes the assumption that the form of the positive and
2541  // negative PDFs are the same, which seems reasonable to me
2542 
2543  LauSPlot::TwoDMap twodimMap;
2544 
2545  for (LauPdfList::const_iterator pdf_iter = negSignalPdfs_.begin(); pdf_iter != negSignalPdfs_.end(); ++pdf_iter) {
2546  // Count the number of input variables that are not DP variables
2547  UInt_t nVars(0);
2548  std::vector<TString> varNames = (*pdf_iter)->varNames();
2549  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
2550  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
2551  ++nVars;
2552  }
2553  }
2554  if ( nVars == 2 ) {
2555  if (useSCF_) {
2556  twodimMap.insert( std::make_pair( "sigTM", std::make_pair( varNames[0], varNames[1] ) ) );
2557  } else {
2558  twodimMap.insert( std::make_pair( "sig", std::make_pair( varNames[0], varNames[1] ) ) );
2559  }
2560  }
2561  }
2562 
2563  if ( useSCF_ ) {
2564  for (LauPdfList::const_iterator pdf_iter = negScfPdfs_.begin(); pdf_iter != negScfPdfs_.end(); ++pdf_iter) {
2565  // Count the number of input variables that are not DP variables
2566  UInt_t nVars(0);
2567  std::vector<TString> varNames = (*pdf_iter)->varNames();
2568  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
2569  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
2570  ++nVars;
2571  }
2572  }
2573  if ( nVars == 2 ) {
2574  twodimMap.insert( std::make_pair( "sigSCF", std::make_pair( varNames[0], varNames[1] ) ) );
2575  }
2576  }
2577  }
2578 
2579  if (usingBkgnd_) {
2580  const UInt_t nBkgnds = this->nBkgndClasses();
2581  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
2582  const TString& bkgndClass = this->bkgndClassName(iBkgnd);
2583  const LauPdfList& pdfList = negBkgndPdfs_[iBkgnd];
2584  for (LauPdfList::const_iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) {
2585  // Count the number of input variables that are not DP variables
2586  UInt_t nVars(0);
2587  std::vector<TString> varNames = (*pdf_iter)->varNames();
2588  for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
2589  if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
2590  ++nVars;
2591  }
2592  }
2593  if ( nVars == 2 ) {
2594  twodimMap.insert( std::make_pair( bkgndClass, std::make_pair( varNames[0], varNames[1] ) ) );
2595  }
2596  }
2597  }
2598  }
2599 
2600  return twodimMap;
2601 }
2602 
2604 {
2605  std::cout << "INFO in LauCPFitModel::storePerEvtLlhds : Storing per-event likelihood values..." << std::endl;
2606 
2607  // if we've not been using the DP model then we need to cache all
2608  // the info here so that we can get the efficiency from it
2609 
2610  LauFitDataTree* inputFitData = this->fitData();
2611 
2612  if (!this->useDP() && this->storeDPEff()) {
2615  negSigModel_->fillDataTree(*inputFitData);
2616  posSigModel_->fillDataTree(*inputFitData);
2617  }
2618 
2619  UInt_t evtsPerExpt(this->eventsPerExpt());
2620 
2621  LauAbsDPDynamics* sigModel(0);
2622  LauPdfList* sigPdfs(0);
2623  LauPdfList* scfPdfs(0);
2624  LauBkgndPdfsList* bkgndPdfs(0);
2625 
2626  for (UInt_t iEvt = 0; iEvt < evtsPerExpt; ++iEvt) {
2627 
2628  this->setSPlotNtupleIntegerBranchValue("iExpt",this->iExpt());
2629  this->setSPlotNtupleIntegerBranchValue("iEvtWithinExpt",iEvt);
2630 
2631  // Find out whether we have B- or B+
2632  if ( tagged_ ) {
2633  const LauFitData& dataValues = inputFitData->getData(iEvt);
2634  LauFitData::const_iterator iter = dataValues.find("charge");
2635  curEvtCharge_ = static_cast<Int_t>(iter->second);
2636 
2637  if (curEvtCharge_==+1) {
2638  sigModel = posSigModel_;
2639  sigPdfs = &posSignalPdfs_;
2640  scfPdfs = &posScfPdfs_;
2641  bkgndPdfs = &posBkgndPdfs_;
2642  } else {
2643  sigModel = negSigModel_;
2644  sigPdfs = &negSignalPdfs_;
2645  scfPdfs = &negScfPdfs_;
2646  bkgndPdfs = &negBkgndPdfs_;
2647  }
2648  } else {
2649  sigPdfs = &negSignalPdfs_;
2650  scfPdfs = &negScfPdfs_;
2651  bkgndPdfs = &negBkgndPdfs_;
2652  }
2653 
2654  // the DP information
2655  this->getEvtDPLikelihood(iEvt);
2656  if (this->storeDPEff()) {
2657  if (!this->useDP()) {
2660  }
2661  if ( tagged_ ) {
2662  this->setSPlotNtupleDoubleBranchValue("efficiency",sigModel->getEvtEff());
2664  this->setSPlotNtupleDoubleBranchValue("scffraction",sigModel->getEvtScfFraction());
2665  }
2666  } else {
2670  }
2671  }
2672  }
2673  if (this->useDP()) {
2675  if (useSCF_) {
2677  this->setSPlotNtupleDoubleBranchValue("sigTMDPLike",sigDPLike_);
2678  this->setSPlotNtupleDoubleBranchValue("sigSCFDPLike",scfDPLike_);
2679  } else {
2680  this->setSPlotNtupleDoubleBranchValue("sigDPLike",sigDPLike_);
2681  }
2682  if (usingBkgnd_) {
2683  const UInt_t nBkgnds = this->nBkgndClasses();
2684  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
2685  TString name = this->bkgndClassName(iBkgnd);
2686  name += "DPLike";
2687  this->setSPlotNtupleDoubleBranchValue(name,bkgndDPLike_[iBkgnd]);
2688  }
2689  }
2690  } else {
2691  sigTotalLike_ = 1.0;
2692  if (useSCF_) {
2693  scfTotalLike_ = 1.0;
2694  }
2695  if (usingBkgnd_) {
2696  const UInt_t nBkgnds = this->nBkgndClasses();
2697  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
2698  bkgndTotalLike_[iBkgnd] = 1.0;
2699  }
2700  }
2701  }
2702 
2703  // the signal PDF values
2704  if ( useSCF_ ) {
2705  sigTotalLike_ *= this->setSPlotNtupleBranchValues(sigPdfs, "sigTM", iEvt);
2706  scfTotalLike_ *= this->setSPlotNtupleBranchValues(scfPdfs, "sigSCF", iEvt);
2707  } else {
2708  sigTotalLike_ *= this->setSPlotNtupleBranchValues(sigPdfs, "sig", iEvt);
2709  }
2710 
2711  // the background PDF values
2712  if (usingBkgnd_) {
2713  const UInt_t nBkgnds = this->nBkgndClasses();
2714  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
2715  const TString& bkgndClass = this->bkgndClassName(iBkgnd);
2716  LauPdfList& pdfs = (*bkgndPdfs)[iBkgnd];
2717  bkgndTotalLike_[iBkgnd] *= this->setSPlotNtupleBranchValues(&(pdfs), bkgndClass, iEvt);
2718  }
2719  }
2720 
2721  // the total likelihoods
2722  if (useSCF_) {
2723  Double_t scfFrac(0.0);
2724  if ( useSCFHist_ ) {
2725  scfFrac = recoSCFFracs_[iEvt];
2726  } else {
2727  scfFrac = scfFrac_.value();
2728  }
2729  this->setSPlotNtupleDoubleBranchValue("sigSCFFrac",scfFrac);
2730  sigTotalLike_ *= ( 1.0 - scfFrac );
2731  if ( scfMap_ == 0 ) {
2732  scfTotalLike_ *= scfFrac;
2733  }
2734  this->setSPlotNtupleDoubleBranchValue("sigTMTotalLike",sigTotalLike_);
2735  this->setSPlotNtupleDoubleBranchValue("sigSCFTotalLike",scfTotalLike_);
2736  } else {
2737  this->setSPlotNtupleDoubleBranchValue("sigTotalLike",sigTotalLike_);
2738  }
2739  if (usingBkgnd_) {
2740  const UInt_t nBkgnds = this->nBkgndClasses();
2741  for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
2742  TString name = this->bkgndClassName(iBkgnd);
2743  name += "TotalLike";
2745  }
2746  }
2747  // fill the tree
2748  this->fillSPlotNtupleBranches();
2749  }
2750  std::cout << "INFO in LauCPFitModel::storePerEvtLlhds : Finished storing per-event likelihood values." << std::endl;
2751 }
2752 
2753 void LauCPFitModel::embedNegSignal(const TString& fileName, const TString& treeName,
2754  Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment,
2755  Bool_t useReweighting)
2756 {
2757  if (negSignalTree_) {
2758  std::cerr << "ERROR in LauCPFitModel::embedNegSignal : Already embedding signal from a file." << std::endl;
2759  return;
2760  }
2761 
2762  if (!reuseEventsWithinEnsemble && reuseEventsWithinExperiment) {
2763  std::cerr << "WARNING in LauCPFitModel::embedNegSignal : Conflicting options provided, will not reuse events at all." << std::endl;
2764  reuseEventsWithinExperiment = kFALSE;
2765  }
2766 
2767  negSignalTree_ = new LauEmbeddedData(fileName,treeName,reuseEventsWithinExperiment);
2768  Bool_t dataOK = negSignalTree_->findBranches();
2769  if (!dataOK) {
2770  delete negSignalTree_; negSignalTree_ = 0;
2771  std::cerr << "ERROR in LauCPFitModel::embedNegSignal : Problem creating data tree for embedding." << std::endl;
2772  return;
2773  }
2774  reuseSignal_ = reuseEventsWithinEnsemble;
2775  useNegReweighting_ = useReweighting;
2776  if (this->enableEmbedding() == kFALSE) {this->enableEmbedding(kTRUE);}
2777 }
2778 
2779 void LauCPFitModel::embedNegBkgnd(const TString& bkgndClass, const TString& fileName, const TString& treeName,
2780  Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment)
2781 {
2782  if ( ! this->validBkgndClass( bkgndClass ) ) {
2783  std::cerr << "ERROR in LauCPFitModel::embedBkgnd : Invalid background class \"" << bkgndClass << "\"." << std::endl;
2784  std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
2785  return;
2786  }
2787 
2788  UInt_t bkgndID = this->bkgndClassID( bkgndClass );
2789 
2790  if (negBkgndTree_[bkgndID]) {
2791  std::cerr << "ERROR in LauCPFitModel::embedNegBkgnd : Already embedding background from a file." << std::endl;
2792  return;
2793  }
2794 
2795  if (!reuseEventsWithinEnsemble && reuseEventsWithinExperiment) {
2796  std::cerr << "WARNING in LauCPFitModel::embedNegBkgnd : Conflicting options provided, will not reuse events at all." << std::endl;
2797  reuseEventsWithinExperiment = kFALSE;
2798  }
2799 
2800  negBkgndTree_[bkgndID] = new LauEmbeddedData(fileName,treeName,reuseEventsWithinExperiment);
2801  Bool_t dataOK = negBkgndTree_[bkgndID]->findBranches();
2802  if (!dataOK) {
2803  delete negBkgndTree_[bkgndID]; negBkgndTree_[bkgndID] = 0;
2804  std::cerr << "ERROR in LauCPFitModel::embedNegBkgnd : Problem creating data tree for embedding." << std::endl;
2805  return;
2806  }
2807  reuseBkgnd_[bkgndID] = reuseEventsWithinEnsemble;
2808  if (this->enableEmbedding() == kFALSE) {this->enableEmbedding(kTRUE);}
2809 }
2810 
2811 void LauCPFitModel::embedPosSignal(const TString& fileName, const TString& treeName,
2812  Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment,
2813  Bool_t useReweighting)
2814 {
2815  if (posSignalTree_) {
2816  std::cerr << "ERROR in LauCPFitModel::embedPosSignal : Already embedding signal from a file." << std::endl;
2817  return;
2818  }
2819 
2820  if (!reuseEventsWithinEnsemble && reuseEventsWithinExperiment) {
2821  std::cerr << "WARNING in LauCPFitModel::embedPosSignal : Conflicting options provided, will not reuse events at all." << std::endl;
2822  reuseEventsWithinExperiment = kFALSE;
2823  }
2824 
2825  posSignalTree_ = new LauEmbeddedData(fileName,treeName,reuseEventsWithinExperiment);
2826  Bool_t dataOK = posSignalTree_->findBranches();
2827  if (!dataOK) {
2828  delete posSignalTree_; posSignalTree_ = 0;
2829  std::cerr << "ERROR in LauCPFitModel::embedPosSignal : Problem creating data tree for embedding." << std::endl;
2830  return;
2831  }
2832  reuseSignal_ = reuseEventsWithinEnsemble;
2833  usePosReweighting_ = useReweighting;
2834  if (this->enableEmbedding() == kFALSE) {this->enableEmbedding(kTRUE);}
2835 }
2836 
2837 void LauCPFitModel::embedPosBkgnd(const TString& bkgndClass, const TString& fileName, const TString& treeName,
2838  Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment)
2839 {
2840  if ( ! this->validBkgndClass( bkgndClass ) ) {
2841  std::cerr << "ERROR in LauCPFitModel::embedBkgnd : Invalid background class \"" << bkgndClass << "\"." << std::endl;
2842  std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
2843  return;
2844  }
2845 
2846  UInt_t bkgndID = this->bkgndClassID( bkgndClass );
2847 
2848  if (posBkgndTree_[bkgndID]) {
2849  std::cerr << "ERROR in LauCPFitModel::embedPosBkgnd : Already embedding background from a file." << std::endl;
2850  return;
2851  }
2852 
2853  if (!reuseEventsWithinEnsemble && reuseEventsWithinExperiment) {
2854  std::cerr << "WARNING in LauCPFitModel::embedPosBkgnd : Conflicting options provided, will not reuse events at all." << std::endl;
2855  reuseEventsWithinExperiment = kFALSE;
2856  }
2857 
2858  posBkgndTree_[bkgndID] = new LauEmbeddedData(fileName,treeName,reuseEventsWithinExperiment);
2859  Bool_t dataOK = posBkgndTree_[bkgndID]->findBranches();
2860  if (!dataOK) {
2861  delete posBkgndTree_[bkgndID]; posBkgndTree_[bkgndID] = 0;
2862  std::cerr << "ERROR in LauCPFitModel::embedPosBkgnd : Problem creating data tree for embedding." << std::endl;
2863  return;
2864  }
2865  reuseBkgnd_[bkgndID] = reuseEventsWithinEnsemble;
2866  if (this->enableEmbedding() == kFALSE) {this->enableEmbedding(kTRUE);}
2867 }
2868 
2869 void LauCPFitModel::weightEvents( const TString& /*dataFileName*/, const TString& /*dataTreeName*/ )
2870 {
2871  std::cerr << "ERROR in LauCPFitModel::weightEvents : Method not available for this fit model." << std::endl;
2872  return;
2873 }
2874 
Double_t range() const
The range allowed for the parameter.
void generateExtraPdfValues(LauPdfList *extraPdfs, LauEmbeddedData *embeddedData)
Generate from the extra PDFs.
Bool_t generateSignalEvent()
Generate signal event.
virtual void addSPlotNtupleIntegerBranch(const TString &name)
Add a branch to the sPlot tree for storing an integer.
Class for defining the abstract interface for signal Dalitz plot dynamics.
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 Double_t getEvtthPrime() const =0
Retrieve the square Dalitz plot coordinate, theta&#39;, for the current event.
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 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.
LauParameter()
Default constructor.
Definition: LauParameter.cc:30
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.
Bool_t useDP() const
Is the Dalitz plot term in the likelihood.
LauDaughters * getDaughters()
Retrieve the daughters.
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-.
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.
std::vector< Double_t > fakeSCFFracs_
The cached values of the SCF fraction for each bin centre.
std::vector< std::vector< LauParameter > > LauParArray
Type to define an array of parameters.
LauParameter getDPRate() const
Retrieve the overall Dalitz plot rate.
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.
virtual void calcLikelihoodInfo(Double_t m13Sq, Double_t m23Sq)=0
Calculate the likelihood (and all associated information) given values of the Dalitz plot coordinates...
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.
Bool_t squareDP() const
Determine to use or not the square Dalitz plot.
Definition: LauDaughters.hh:72
virtual Bool_t hasResonance(const TString &resName) const
Check whether this model includes a named resonance.
std::vector< LauParameter > getExtraParameters()
Retrieve any extra parameters/quantities (e.g. K-matrix total fit fractions)
LauParameter negDPRate_
The average DP rate for B-.
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.
std::vector< Double_t > bkgndTotalLike_
Total background likelihood(s)
virtual void fillGenNtupleBranches()
Fill the gen tuple branches.
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.
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.
LauParameter getMeanEff() const
Retrieve the mean efficiency across the Dalitz plot.
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.
LauAbsDPDynamics * posSigModel_
The B+ signal Dalitz plot model.
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.
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. ...
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.
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:229
std::vector< LauParameter * > LauParameterPList
List of parameter pointers.
virtual void fillDataTree(const LauFitDataTree &inputFitTree)=0
Obtain data from a fit tree.
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.
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)
virtual void calcExtraInfo(Bool_t init=kFALSE)=0
Calculate the fit fractions, mean efficiency and total DP rate.
Class for defining a CP fit model.
Abstract interface to the fitting and toy MC model.
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 Double_t getEvtScfFraction() const =0
Retrieve the fraction of events that are poorly reconstructed (the self cross feed fraction) in the D...
virtual void addGenNtupleIntegerBranch(const TString &name)
Add a branch to the gen tree for storing an integer.
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.
LauScfMap * scfMap_
The smearing matrix for the SCF DP PDF.
File containing declaration of LauComplex class.
Double_t sigExtraLike_
Signal likelihood from extra PDFs.
virtual TString getConjResName(const TString &resName) const
Retrieve the name of the charge conjugate of a named resonance.
LauFitData getValues(const std::vector< TString > &names) const
Get values of specified branches.
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:33
Class to store the data for embedding in toy experiments.
virtual Double_t getEvtmPrime() const =0
Retrieve the square Dalitz plot coordinate, m&#39;, for the current event.
LauAbsEffModel * getEffModel()
Retrieve the model for the efficiency across the Dalitz plot.
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< 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.
virtual ToyMCStatus checkToyMC(Bool_t printErrorMessages=kTRUE, Bool_t printInfoMessages=kFALSE)=0
Check the status of the toy MC generation.
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.
LauAbsDPDynamics * negSigModel_
The B- signal Dalitz plot 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.
UInt_t getnAmp() const
Retrieve the number of amplitude components.
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.
File containing declaration of LauAbsDPDynamics class.
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 Double_t getEvtm23Sq() const =0
Retrieve the invariant mass squared of the second and third daughters in the current 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.
Double_t getDPNorm() const
Retrieve the normalisation factor for the log-likelihood function.
Class that implements the efficiency description across the signal Dalitz plot.
Definition: LauEffModel.hh:37
TString getNameParent() const
Get name of the parent particle.
UInt_t nUsedEvents() const
Retrieve the number of events that have already been sampled.
virtual void updateCoeffs(const std::vector< LauComplex > &coeffs)
Update the complex coefficients for the resonances.
const LauParArray & getFitFractionsEfficiencyUncorrected() const
Retrieve the fit fractions for the amplitude components.
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.
virtual Double_t getEvtm13Sq() const =0
Retrieve the invariant mass squared of the first and third daughters in the current event...
Bool_t getReweightedEvent(LauAbsDPDynamics *dynamics)
Retrieve an event from the data sample, applying an accept/reject based on the given DP 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
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.
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
const TMatrixD & covarianceMatrix() const
Access the fit covariance matrix.
const LauParameterPList & fitPars() const
Access the fit variables.
LauPdfList posScfPdfs_
The B+ SCF PDFs.
virtual const LauComplex & getEvtDPAmp() const =0
Retrieve the total amplitude of all amplitude components at the current point in the Dalitz plot...
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.
const LauParArray & getFitFractions() const
Retrieve the fit fractions for the amplitude components.
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.
virtual Double_t getEvtEff() const =0
Retrieve the efficiency for the current event.
Double_t getAsymmetry()
Retrieve the asymmetry.
Definition: LauAsymmCalc.hh:42
virtual Double_t getEvtSCFDPLikelihood(UInt_t iEvt)
Calculate the SCF likelihood for the DP for a given event.
Bool_t usingScfModel()
Check whether a self cross feed fraction model is being used.
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
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 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.
virtual void initialise(const std::vector< LauComplex > &coeffs)=0
Initialise the Dalitz plot dynamics.
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 generate()=0
Generate a toy MC signal event.
virtual Bool_t passVeto(const LauKinematics *kinematics) const =0
Determine whether the given DP position is outside the vetoes.