laura is hosted by Hepforge, IPPP Durham
Laura++  3.6.0
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauCPFitModel.cc
Go to the documentation of this file.
1 
2 /*
3 Copyright 2004 University of Warwick
4 
5 Licensed under the Apache License, Version 2.0 (the "License");
6 you may not use this file except in compliance with the License.
7 You may obtain a copy of the License at
8 
9  http://www.apache.org/licenses/LICENSE-2.0
10 
11 Unless required by applicable law or agreed to in writing, software
12 distributed under the License is distributed on an "AS IS" BASIS,
13 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 See the License for the specific language governing permissions and
15 limitations under the License.
16 */
17 
18 /*
19 Laura++ package authors:
20 John Back
21 Paul Harrison
22 Thomas Latham
23 */
24 
29 #include "LauCPFitModel.hh"
30 
31 #include "LauAbsBkgndDPModel.hh"
32 #include "LauAbsCoeffSet.hh"
33 #include "LauAbsPdf.hh"
34 #include "LauAsymmCalc.hh"
35 #include "LauComplex.hh"
36 #include "LauConstants.hh"
37 #include "LauDaughters.hh"
38 #include "LauEffModel.hh"
39 #include "LauEmbeddedData.hh"
40 #include "LauFitNtuple.hh"
41 #include "LauGenNtuple.hh"
42 #include "LauIsobarDynamics.hh"
43 #include "LauKinematics.hh"
44 #include "LauPrint.hh"
45 #include "LauRandom.hh"
46 #include "LauScfMap.hh"
47 
48 #include "TBranch.h"
49 #include "TCanvas.h"
50 #include "TFile.h"
51 #include "TGraph.h"
52 #include "TGraph2D.h"
53 #include "TH2.h"
54 #include "TLeaf.h"
55 #include "TMath.h"
56 #include "TMinuit.h"
57 #include "TRandom.h"
58 #include "TStyle.h"
59 #include "TSystem.h"
60 #include "TTree.h"
61 #include "TVirtualFitter.h"
62 
63 #include <fstream>
64 #include <iomanip>
65 #include <iostream>
66 #include <vector>
67 
68 LauCPFitModel::LauCPFitModel( LauIsobarDynamics* negModel,
69  LauIsobarDynamics* posModel,
70  Bool_t tagged,
71  const TString& tagVarName ) :
72  LauAbsFitModel(),
73  negSigModel_( negModel ),
74  posSigModel_( posModel ),
75  negKinematics_( negModel ? negModel->getKinematics() : 0 ),
76  posKinematics_( posModel ? posModel->getKinematics() : 0 ),
77  usingBkgnd_( kFALSE ),
78  nSigComp_( 0 ),
79  nSigDPPar_( 0 ),
80  nExtraPdfPar_( 0 ),
81  nNormPar_( 0 ),
82  negMeanEff_( "negMeanEff", 0.0, 0.0, 1.0 ),
83  posMeanEff_( "posMeanEff", 0.0, 0.0, 1.0 ),
84  negDPRate_( "negDPRate", 0.0, 0.0, 100.0 ),
85  posDPRate_( "posDPRate", 0.0, 0.0, 100.0 ),
86  signalEvents_( 0 ),
87  signalAsym_( 0 ),
88  forceAsym_( kFALSE ),
89  tagged_( tagged ),
90  tagVarName_( tagVarName ),
91  curEvtCharge_( 0 ),
92  useSCF_( kFALSE ),
93  useSCFHist_( kFALSE ),
94  scfFrac_( "scfFrac", 0.0, 0.0, 1.0 ),
95  scfFracHist_( 0 ),
96  scfMap_( 0 ),
97  compareFitData_( kFALSE ),
98  negParent_( "B-" ),
99  posParent_( "B+" ),
100  negSignalTree_( 0 ),
101  posSignalTree_( 0 ),
102  reuseSignal_( kFALSE ),
103  useNegReweighting_( kFALSE ),
104  usePosReweighting_( kFALSE ),
105  sigDPLike_( 0.0 ),
106  scfDPLike_( 0.0 ),
107  sigExtraLike_( 0.0 ),
108  scfExtraLike_( 0.0 ),
109  sigTotalLike_( 0.0 ),
110  scfTotalLike_( 0.0 )
111 {
112  const LauDaughters* negDaug = negSigModel_->getDaughters();
113  if ( negDaug != 0 ) {
114  negParent_ = negDaug->getNameParent();
115  }
116  const LauDaughters* posDaug = posSigModel_->getDaughters();
117  if ( posDaug != 0 ) {
118  posParent_ = posDaug->getNameParent();
119  }
120 }
121 
122 LauCPFitModel::~LauCPFitModel()
123 {
124  delete negSignalTree_;
125  delete posSignalTree_;
126  for ( LauBkgndEmbDataList::iterator iter = negBkgndTree_.begin(); iter != negBkgndTree_.end();
127  ++iter ) {
128  delete ( *iter );
129  }
130  for ( LauBkgndEmbDataList::iterator iter = posBkgndTree_.begin(); iter != posBkgndTree_.end();
131  ++iter ) {
132  delete ( *iter );
133  }
134  delete scfFracHist_;
135 }
136 
137 void LauCPFitModel::setupBkgndVectors()
138 {
139  UInt_t nBkgnds = this->nBkgndClasses();
140  negBkgndDPModels_.resize( nBkgnds );
141  posBkgndDPModels_.resize( nBkgnds );
142  negBkgndPdfs_.resize( nBkgnds );
143  posBkgndPdfs_.resize( nBkgnds );
144  bkgndEvents_.resize( nBkgnds );
145  bkgndAsym_.resize( nBkgnds );
146  negBkgndTree_.resize( nBkgnds );
147  posBkgndTree_.resize( nBkgnds );
148  reuseBkgnd_.resize( nBkgnds );
149  bkgndDPLike_.resize( nBkgnds );
150  bkgndExtraLike_.resize( nBkgnds );
151  bkgndTotalLike_.resize( nBkgnds );
152 }
153 
154 void LauCPFitModel::setNSigEvents( LauParameter* nSigEvents )
155 {
156  if ( nSigEvents == 0 ) {
157  std::cerr << "ERROR in LauCPFitModel::setNSigEvents : The LauParameter pointer is null."
158  << std::endl;
159  gSystem->Exit( EXIT_FAILURE );
160  }
161  if ( signalEvents_ != 0 ) {
162  std::cerr << "ERROR in LauCPFitModel::setNSigEvents : You are trying to overwrite the signal yield."
163  << std::endl;
164  return;
165  }
166  if ( signalAsym_ != 0 ) {
167  std::cerr << "ERROR in LauCPFitModel::setNSigEvents : You are trying to overwrite the signal asymmetry."
168  << std::endl;
169  return;
170  }
171 
172  signalEvents_ = nSigEvents;
173  TString name = signalEvents_->name();
174  if ( ! name.Contains( "signalEvents" ) &&
175  ! ( name.BeginsWith( "signal" ) && name.EndsWith( "Events" ) ) ) {
176  signalEvents_->name( "signalEvents" );
177  }
178  Double_t value = nSigEvents->value();
179  signalEvents_->range( -2.0 * ( TMath::Abs( value ) + 1.0 ), 2.0 * ( TMath::Abs( value ) + 1.0 ) );
180 
181  signalAsym_ = new LauParameter( "signalAsym", 0.0, -1.0, 1.0, kTRUE );
182 }
183 
184 void LauCPFitModel::setNSigEvents( LauParameter* nSigEvents, LauParameter* sigAsym, Bool_t forceAsym )
185 {
186  if ( nSigEvents == 0 ) {
187  std::cerr << "ERROR in LauCPFitModel::setNSigEvents : The event LauParameter pointer is null."
188  << std::endl;
189  gSystem->Exit( EXIT_FAILURE );
190  }
191  if ( sigAsym == 0 ) {
192  std::cerr << "ERROR in LauCPFitModel::setNSigEvents : The asym LauParameter pointer is null."
193  << std::endl;
194  gSystem->Exit( EXIT_FAILURE );
195  }
196 
197  if ( signalEvents_ != 0 ) {
198  std::cerr << "ERROR in LauCPFitModel::setNSigEvents : You are trying to overwrite the signal yield."
199  << std::endl;
200  return;
201  }
202  if ( signalAsym_ != 0 ) {
203  std::cerr << "ERROR in LauCPFitModel::setNSigEvents : You are trying to overwrite the signal asymmetry."
204  << std::endl;
205  return;
206  }
207 
208  signalEvents_ = nSigEvents;
209  signalEvents_->name( "signalEvents" );
210  Double_t value = nSigEvents->value();
211  signalEvents_->range( -2.0 * ( TMath::Abs( value ) + 1.0 ), 2.0 * ( TMath::Abs( value ) + 1.0 ) );
212 
213  signalAsym_ = sigAsym;
214  signalAsym_->name( "signalAsym" );
215  signalAsym_->range( -1.0, 1.0 );
216 
217  forceAsym_ = forceAsym;
218 }
219 
220 void LauCPFitModel::setNBkgndEvents( LauAbsRValue* nBkgndEvents )
221 {
222  if ( nBkgndEvents == 0 ) {
223  std::cerr << "ERROR in LauCPFitModel::setNBgkndEvents : The background yield LauParameter pointer is null."
224  << std::endl;
225  gSystem->Exit( EXIT_FAILURE );
226  }
227 
228  if ( ! this->validBkgndClass( nBkgndEvents->name() ) ) {
229  std::cerr << "ERROR in LauCPFitModel::setNBkgndEvents : Invalid background class \""
230  << nBkgndEvents->name() << "\"." << std::endl;
231  std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed."
232  << std::endl;
233  gSystem->Exit( EXIT_FAILURE );
234  }
235 
236  UInt_t bkgndID = this->bkgndClassID( nBkgndEvents->name() );
237 
238  if ( bkgndEvents_[bkgndID] != 0 ) {
239  std::cerr << "ERROR in LauCPFitModel::setNBkgndEvents : You are trying to overwrite the background yield."
240  << std::endl;
241  return;
242  }
243 
244  if ( bkgndAsym_[bkgndID] != 0 ) {
245  std::cerr << "ERROR in LauCPFitModel::setNBkgndEvents : You are trying to overwrite the background asymmetry."
246  << std::endl;
247  return;
248  }
249 
250  nBkgndEvents->name( nBkgndEvents->name() + "Events" );
251  if ( nBkgndEvents->isLValue() ) {
252  Double_t value = nBkgndEvents->value();
253  LauParameter* yield = dynamic_cast<LauParameter*>( nBkgndEvents );
254  yield->range( -2.0 * ( TMath::Abs( value ) + 1.0 ), 2.0 * ( TMath::Abs( value ) + 1.0 ) );
255  }
256  bkgndEvents_[bkgndID] = nBkgndEvents;
257 
258  bkgndAsym_[bkgndID] = new LauParameter( nBkgndEvents->name() + "Asym", 0.0, -1.0, 1.0, kTRUE );
259 }
260 
261 void LauCPFitModel::setNBkgndEvents( LauAbsRValue* nBkgndEvents, LauAbsRValue* bkgndAsym )
262 {
263  if ( nBkgndEvents == 0 ) {
264  std::cerr << "ERROR in LauCPFitModel::setNBkgndEvents : The background yield LauParameter pointer is null."
265  << std::endl;
266  gSystem->Exit( EXIT_FAILURE );
267  }
268 
269  if ( bkgndAsym == 0 ) {
270  std::cerr << "ERROR in LauCPFitModel::setNBkgndEvents : The background asym LauParameter pointer is null."
271  << std::endl;
272  gSystem->Exit( EXIT_FAILURE );
273  }
274 
275  if ( ! this->validBkgndClass( nBkgndEvents->name() ) ) {
276  std::cerr << "ERROR in LauCPFitModel::setNBkgndEvents : Invalid background class \""
277  << nBkgndEvents->name() << "\"." << std::endl;
278  std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed."
279  << std::endl;
280  gSystem->Exit( EXIT_FAILURE );
281  }
282 
283  UInt_t bkgndID = this->bkgndClassID( nBkgndEvents->name() );
284 
285  if ( bkgndEvents_[bkgndID] != 0 ) {
286  std::cerr << "ERROR in LauCPFitModel::setNBkgndEvents : You are trying to overwrite the background yield."
287  << std::endl;
288  return;
289  }
290 
291  if ( bkgndAsym_[bkgndID] != 0 ) {
292  std::cerr << "ERROR in LauCPFitModel::setNBkgndEvents : You are trying to overwrite the background asymmetry."
293  << std::endl;
294  return;
295  }
296 
297  nBkgndEvents->name( nBkgndEvents->name() + "Events" );
298  if ( nBkgndEvents->isLValue() ) {
299  Double_t value = nBkgndEvents->value();
300  LauParameter* yield = dynamic_cast<LauParameter*>( nBkgndEvents );
301  yield->range( -2.0 * ( TMath::Abs( value ) + 1.0 ), 2.0 * ( TMath::Abs( value ) + 1.0 ) );
302  }
303  bkgndEvents_[bkgndID] = nBkgndEvents;
304 
305  bkgndAsym->name( nBkgndEvents->name() + "Asym" );
306  if ( bkgndAsym->isLValue() ) {
307  LauParameter* asym = dynamic_cast<LauParameter*>( bkgndAsym );
308  asym->range( -1.0, 1.0 );
309  }
310  bkgndAsym_[bkgndID] = bkgndAsym;
311 }
312 
313 void LauCPFitModel::splitSignalComponent( const TH2* dpHisto,
314  const Bool_t upperHalf,
315  const Bool_t fluctuateBins,
316  LauScfMap* scfMap )
317 {
318  if ( useSCF_ == kTRUE ) {
319  std::cerr << "ERROR in LauCPFitModel::splitSignalComponent : Have already setup SCF."
320  << std::endl;
321  return;
322  }
323 
324  if ( dpHisto == 0 ) {
325  std::cerr << "ERROR in LauCPFitModel::splitSignalComponent : The histogram pointer is null."
326  << std::endl;
327  return;
328  }
329 
330  const LauDaughters* daughters = negSigModel_->getDaughters();
331  scfFracHist_ = new LauEffModel( daughters, 0 );
332  scfFracHist_
333  ->setEffHisto( dpHisto, kTRUE, fluctuateBins, 0.0, 0.0, upperHalf, daughters->squareDP() );
334 
335  scfMap_ = scfMap;
336 
337  useSCF_ = kTRUE;
338  useSCFHist_ = kTRUE;
339 }
340 
341 void LauCPFitModel::splitSignalComponent( const Double_t scfFrac, const Bool_t fixed )
342 {
343  if ( useSCF_ == kTRUE ) {
344  std::cerr << "ERROR in LauCPFitModel::splitSignalComponent : Have already setup SCF."
345  << std::endl;
346  return;
347  }
348 
349  scfFrac_.range( 0.0, 1.0 );
350  scfFrac_.value( scfFrac );
351  scfFrac_.initValue( scfFrac );
352  scfFrac_.genValue( scfFrac );
353  scfFrac_.fixed( fixed );
354 
355  useSCF_ = kTRUE;
356  useSCFHist_ = kFALSE;
357 }
358 
359 void LauCPFitModel::setBkgndDPModels( const TString& bkgndClass,
360  LauAbsBkgndDPModel* negModel,
361  LauAbsBkgndDPModel* posModel )
362 {
363  if ( ( negModel == 0 ) || ( posModel == 0 ) ) {
364  std::cerr << "ERROR in LauCPFitModel::setBkgndDPModels : One or both of the model pointers is null."
365  << std::endl;
366  return;
367  }
368 
369  // check that this background name is valid
370  if ( ! this->validBkgndClass( bkgndClass ) ) {
371  std::cerr << "ERROR in LauCPFitModel::setBkgndDPModel : Invalid background class \""
372  << bkgndClass << "\"." << std::endl;
373  std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed."
374  << std::endl;
375  return;
376  }
377 
378  UInt_t bkgndID = this->bkgndClassID( bkgndClass );
379  negBkgndDPModels_[bkgndID] = negModel;
380  posBkgndDPModels_[bkgndID] = posModel;
381 
382  usingBkgnd_ = kTRUE;
383 }
384 
385 void LauCPFitModel::setSignalPdfs( LauAbsPdf* negPdf, LauAbsPdf* posPdf )
386 {
387  if ( tagged_ ) {
388  if ( negPdf == 0 || posPdf == 0 ) {
389  std::cerr << "ERROR in LauCPFitModel::setSignalPdfs : One or both of the PDF pointers is null."
390  << std::endl;
391  return;
392  }
393  } else {
394  // if we're doing an untagged analysis we will only use the negative PDFs
395  if ( negPdf == 0 ) {
396  std::cerr << "ERROR in LauCPFitModel::setSignalPdfs : The negative PDF pointer is null."
397  << std::endl;
398  return;
399  }
400  if ( posPdf != 0 ) {
401  std::cerr << "WARNING in LauCPFitModel::setSignalPdfs : Doing an untagged fit so will not use the positive PDF."
402  << std::endl;
403  }
404  }
405  negSignalPdfs_.push_back( negPdf );
406  posSignalPdfs_.push_back( posPdf );
407 }
408 
409 void LauCPFitModel::setSCFPdfs( LauAbsPdf* negPdf, LauAbsPdf* posPdf )
410 {
411  if ( tagged_ ) {
412  if ( negPdf == 0 || posPdf == 0 ) {
413  std::cerr << "ERROR in LauCPFitModel::setSCFPdfs : One or both of the PDF pointers is null."
414  << std::endl;
415  return;
416  }
417  } else {
418  // if we're doing an untagged analysis we will only use the negative PDFs
419  if ( negPdf == 0 ) {
420  std::cerr << "ERROR in LauCPFitModel::setSCFPdfs : The negative PDF pointer is null."
421  << std::endl;
422  return;
423  }
424  if ( posPdf != 0 ) {
425  std::cerr << "WARNING in LauCPFitModel::setSCFPdfs : Doing an untagged fit so will not use the positive PDF."
426  << std::endl;
427  }
428  }
429  negScfPdfs_.push_back( negPdf );
430  posScfPdfs_.push_back( posPdf );
431 }
432 
433 void LauCPFitModel::setBkgndPdfs( const TString& bkgndClass, LauAbsPdf* negPdf, LauAbsPdf* posPdf )
434 {
435  if ( tagged_ ) {
436  if ( negPdf == 0 || posPdf == 0 ) {
437  std::cerr << "ERROR in LauCPFitModel::setBkgndPdfs : One or both of the PDF pointers is null."
438  << std::endl;
439  return;
440  }
441  } else {
442  // if we're doing an untagged analysis we will only use the negative PDFs
443  if ( negPdf == 0 ) {
444  std::cerr << "ERROR in LauCPFitModel::setBkgndPdfs : The negative PDF pointer is null."
445  << std::endl;
446  return;
447  }
448  if ( posPdf != 0 ) {
449  std::cerr << "WARNING in LauCPFitModel::setBkgndPdfs : Doing an untagged fit so will not use the positive PDF."
450  << std::endl;
451  }
452  }
453 
454  // check that this background name is valid
455  if ( ! this->validBkgndClass( bkgndClass ) ) {
456  std::cerr << "ERROR in LauCPFitModel::setBkgndPdfs : Invalid background class \""
457  << bkgndClass << "\"." << std::endl;
458  std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed."
459  << std::endl;
460  return;
461  }
462 
463  UInt_t bkgndID = this->bkgndClassID( bkgndClass );
464  negBkgndPdfs_[bkgndID].push_back( negPdf );
465  posBkgndPdfs_[bkgndID].push_back( posPdf );
466 
467  usingBkgnd_ = kTRUE;
468 }
469 
470 void LauCPFitModel::setAmpCoeffSet( LauAbsCoeffSet* coeffSet )
471 {
472  // Resize the coeffPars vector if not already done
473  if ( coeffPars_.empty() ) {
474  const UInt_t nNegAmp = negSigModel_->getnTotAmp();
475  const UInt_t nPosAmp = posSigModel_->getnTotAmp();
476  if ( nNegAmp != nPosAmp ) {
477  std::cerr << "ERROR in LauCPFitModel::setAmpCoeffSet : Unequal number of signal DP components in the negative and positive models: "
478  << nNegAmp << " != " << nPosAmp << std::endl;
479  gSystem->Exit( EXIT_FAILURE );
480  }
481  coeffPars_.resize( nNegAmp );
482  for ( std::vector<LauAbsCoeffSet*>::iterator iter = coeffPars_.begin();
483  iter != coeffPars_.end();
484  ++iter ) {
485  ( *iter ) = 0;
486  }
487  fitFracAsymm_.resize( nNegAmp );
488  acp_.resize( nNegAmp );
489  }
490 
491  // Is there a component called compName in the signal model?
492  TString compName( coeffSet->name() );
493  TString conjName = negSigModel_->getConjResName( compName );
494  const Int_t negIndex = negSigModel_->resonanceIndex( compName );
495  const Int_t posIndex = posSigModel_->resonanceIndex( conjName );
496  if ( negIndex < 0 ) {
497  std::cerr << "ERROR in LauCPFitModel::setAmpCoeffSet : " << negParent_
498  << " signal DP model doesn't contain component \"" << compName << "\"."
499  << std::endl;
500  return;
501  }
502  if ( posIndex < 0 ) {
503  std::cerr << "ERROR in LauCPFitModel::setAmpCoeffSet : " << posParent_
504  << " signal DP model doesn't contain component \"" << conjName << "\"."
505  << std::endl;
506  return;
507  }
508  if ( posIndex != negIndex ) {
509  std::cerr << "ERROR in LauCPFitModel::setAmpCoeffSet : " << negParent_
510  << " signal DP model and " << posParent_
511  << " signal DP model have different indices for components \"" << compName
512  << "\" and \"" << conjName << "\"." << std::endl;
513  return;
514  }
515 
516  // Do we already have it in our list of names?
517  if ( coeffPars_[negIndex] != 0 && coeffPars_[negIndex]->name() == compName ) {
518  std::cerr << "ERROR in LauCPFitModel::setAmpCoeffSet : Have already set coefficients for \""
519  << compName << "\"." << std::endl;
520  return;
521  }
522 
523  coeffSet->index( negIndex );
524  coeffPars_[negIndex] = coeffSet;
525 
526  TString parName = coeffSet->baseName();
527  parName += "FitFracAsym";
528  fitFracAsymm_[negIndex] = LauParameter( parName, 0.0, -1.0, 1.0 );
529 
530  acp_[negIndex] = coeffSet->acp();
531 
532  ++nSigComp_;
533 
534  std::cout << "INFO in LauCPFitModel::setAmpCoeffSet : Added coefficients for component \""
535  << compName << "\" to the fit model." << std::endl;
536  coeffSet->printParValues();
537 }
538 
539 void LauCPFitModel::initialise()
540 {
541  // Initialisation
542  if ( ! this->useDP() && negSignalPdfs_.empty() ) {
543  std::cerr << "ERROR in LauCPFitModel::initialise : Signal model doesn't exist for any variable."
544  << std::endl;
545  gSystem->Exit( EXIT_FAILURE );
546  }
547 
548  if ( this->useDP() ) {
549  // Check that we have all the Dalitz-plot models
550  if ( ( negSigModel_ == 0 ) || ( posSigModel_ == 0 ) ) {
551  std::cerr << "ERROR in LauCPFitModel::initialise : the pointer to one (neg or pos) of the signal DP models is null.\n";
552  std::cerr << " : Removing the Dalitz Plot from the model."
553  << std::endl;
554  this->useDP( kFALSE );
555  }
556  if ( usingBkgnd_ ) {
557  if ( negBkgndDPModels_.empty() || posBkgndDPModels_.empty() ) {
558  std::cerr << "ERROR in LauCPFitModel::initialise : No background DP models found.\n";
559  std::cerr << " : Removing the Dalitz plot from the model."
560  << std::endl;
561  this->useDP( kFALSE );
562  }
563  for ( LauBkgndDPModelList::const_iterator dpmodel_iter = negBkgndDPModels_.begin();
564  dpmodel_iter != negBkgndDPModels_.end();
565  ++dpmodel_iter ) {
566  if ( ( *dpmodel_iter ) == 0 ) {
567  std::cerr << "ERROR in LauCPFitModel::initialise : The pointer to one of the background DP models is null.\n";
568  std::cerr << " : Removing the Dalitz Plot from the model."
569  << std::endl;
570  this->useDP( kFALSE );
571  break;
572  }
573  }
574  for ( LauBkgndDPModelList::const_iterator dpmodel_iter = posBkgndDPModels_.begin();
575  dpmodel_iter != posBkgndDPModels_.end();
576  ++dpmodel_iter ) {
577  if ( ( *dpmodel_iter ) == 0 ) {
578  std::cerr << "ERROR in LauCPFitModel::initialise : The pointer to one of the background DP models is null.\n";
579  std::cerr << " : Removing the Dalitz Plot from the model."
580  << std::endl;
581  this->useDP( kFALSE );
582  break;
583  }
584  }
585  }
586 
587  // Need to check that the number of components we have and that the dynamics has matches up
588  const UInt_t nNegAmp = negSigModel_->getnTotAmp();
589  const UInt_t nPosAmp = posSigModel_->getnTotAmp();
590  if ( nNegAmp != nPosAmp ) {
591  std::cerr << "ERROR in LauCPFitModel::initialise : Unequal number of signal DP components in the negative and positive models: "
592  << nNegAmp << " != " << nPosAmp << std::endl;
593  gSystem->Exit( EXIT_FAILURE );
594  }
595  if ( nNegAmp != nSigComp_ ) {
596  std::cerr << "ERROR in LauCPFitModel::initialise : Number of signal DP components in the model ("
597  << nNegAmp << ") not equal to number of coefficients supplied (" << nSigComp_
598  << ")." << std::endl;
599  gSystem->Exit( EXIT_FAILURE );
600  }
601 
602  if ( ! fixParamFileName_.IsNull() || ! fixParamMap_.empty() ) {
603 
604  // Set coefficients
605 
606  std::vector<LauParameter*> params;
607  for ( auto itr = coeffPars_.begin(); itr != coeffPars_.end(); ++itr ) {
608  std::vector<LauParameter*> p = ( *itr )->getParameters();
609  params.insert( params.end(), p.begin(), p.end() );
610  }
611 
612  this->fixParams( params );
613 
614  // Set resonance parameters (if they exist)
615 
616  negSigModel_->collateResonanceParameters();
617  posSigModel_->collateResonanceParameters();
618 
619  this->fixParams( negSigModel_->getFloatingParameters() );
620  this->fixParams( posSigModel_->getFloatingParameters() );
621  }
622 
623  // From the initial parameter values calculate the coefficients
624  // so they can be passed to the signal model
625  this->updateCoeffs();
626 
627  // If all is well, go ahead and initialise them
628  this->initialiseDPModels();
629  }
630 
631  // Next check that, if a given component is being used we've got the
632  // right number of PDFs for all the variables involved
633  // TODO - should probably check variable names and so on as well
634 
635  UInt_t nsigpdfvars( 0 );
636  for ( LauPdfPList::const_iterator pdf_iter = negSignalPdfs_.begin();
637  pdf_iter != negSignalPdfs_.end();
638  ++pdf_iter ) {
639  std::vector<TString> varNames = ( *pdf_iter )->varNames();
640  for ( std::vector<TString>::const_iterator var_iter = varNames.begin();
641  var_iter != varNames.end();
642  ++var_iter ) {
643  if ( ( *var_iter ) != "m13Sq" && ( *var_iter ) != "m23Sq" ) {
644  ++nsigpdfvars;
645  }
646  }
647  }
648  if ( useSCF_ ) {
649  UInt_t nscfpdfvars( 0 );
650  for ( LauPdfPList::const_iterator pdf_iter = negScfPdfs_.begin();
651  pdf_iter != negScfPdfs_.end();
652  ++pdf_iter ) {
653  std::vector<TString> varNames = ( *pdf_iter )->varNames();
654  for ( std::vector<TString>::const_iterator var_iter = varNames.begin();
655  var_iter != varNames.end();
656  ++var_iter ) {
657  if ( ( *var_iter ) != "m13Sq" && ( *var_iter ) != "m23Sq" ) {
658  ++nscfpdfvars;
659  }
660  }
661  }
662  if ( nscfpdfvars != nsigpdfvars ) {
663  std::cerr << "ERROR in LauCPFitModel::initialise : There are " << nsigpdfvars
664  << " TM signal PDF variables but " << nscfpdfvars
665  << " SCF signal PDF variables." << std::endl;
666  gSystem->Exit( EXIT_FAILURE );
667  }
668  }
669  if ( usingBkgnd_ ) {
670  for ( LauBkgndPdfsList::const_iterator bgclass_iter = negBkgndPdfs_.begin();
671  bgclass_iter != negBkgndPdfs_.end();
672  ++bgclass_iter ) {
673  UInt_t nbkgndpdfvars( 0 );
674  const LauPdfPList& pdfList = ( *bgclass_iter );
675  for ( LauPdfPList::const_iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end();
676  ++pdf_iter ) {
677  std::vector<TString> varNames = ( *pdf_iter )->varNames();
678  for ( std::vector<TString>::const_iterator var_iter = varNames.begin();
679  var_iter != varNames.end();
680  ++var_iter ) {
681  if ( ( *var_iter ) != "m13Sq" && ( *var_iter ) != "m23Sq" ) {
682  ++nbkgndpdfvars;
683  }
684  }
685  }
686  if ( nbkgndpdfvars != nsigpdfvars ) {
687  std::cerr << "ERROR in LauCPFitModel::initialise : There are " << nsigpdfvars
688  << " signal PDF variables but " << nbkgndpdfvars
689  << " bkgnd PDF variables." << std::endl;
690  gSystem->Exit( EXIT_FAILURE );
691  }
692  }
693  }
694 
695  // Clear the vectors of parameter information so we can start from scratch
696  this->clearFitParVectors();
697 
698  // Set the fit parameters for signal and background models
699  this->setSignalDPParameters();
700 
701  // Set the fit parameters for the various extra PDFs
702  this->setExtraPdfParameters();
703 
704  // Set the initial bg and signal events
705  this->setFitNEvents();
706 
707  // Check that we have the expected number of fit variables
708  const LauParameterPList& fitVars = this->fitPars();
709  if ( fitVars.size() != ( nSigDPPar_ + nExtraPdfPar_ + nNormPar_ ) ) {
710  std::cerr << "ERROR in LauCPFitModel::initialise : Number of fit parameters not of expected size. Exiting"
711  << std::endl;
712  gSystem->Exit( EXIT_FAILURE );
713  }
714 
715  this->setExtraNtupleVars();
716 }
717 
718 void LauCPFitModel::recalculateNormalisation()
719 {
720  //std::cout << "INFO in LauCPFitModel::recalculateNormalizationInDPModels : Recalc Norm in DP model" << std::endl;
721  negSigModel_->recalculateNormalisation();
722  posSigModel_->recalculateNormalisation();
723  negSigModel_->modifyDataTree();
724  posSigModel_->modifyDataTree();
725 }
726 
727 void LauCPFitModel::initialiseDPModels()
728 {
729  std::cout << "INFO in LauCPFitModel::initialiseDPModels : Initialising signal DP model"
730  << std::endl;
731  negSigModel_->initialise( negCoeffs_ );
732  posSigModel_->initialise( posCoeffs_ );
733 
734  if ( usingBkgnd_ == kTRUE ) {
735  for ( LauBkgndDPModelList::iterator iter = negBkgndDPModels_.begin();
736  iter != negBkgndDPModels_.end();
737  ++iter ) {
738  ( *iter )->initialise();
739  }
740  for ( LauBkgndDPModelList::iterator iter = posBkgndDPModels_.begin();
741  iter != posBkgndDPModels_.end();
742  ++iter ) {
743  ( *iter )->initialise();
744  }
745  }
746 }
747 
748 void LauCPFitModel::setSignalDPParameters()
749 {
750  // Set the fit parameters for the signal model.
751 
752  nSigDPPar_ = 0;
753  if ( ! this->useDP() ) {
754  return;
755  }
756 
757  std::cout << "INFO in LauCPFitModel::setSignalDPParameters : Setting the initial fit parameters for the signal DP model."
758  << std::endl;
759 
760  // Place isobar coefficient parameters in vector of fit variables
761  for ( UInt_t i = 0; i < nSigComp_; i++ ) {
762  LauParameterPList pars = coeffPars_[i]->getParameters();
763  nSigDPPar_ += this->addFitParameters( pars, kTRUE );
764  }
765 
766  // Obtain the resonance parameters and place them in the vector of fit variables and in a separate vector
767  // Need to make sure that they are unique because some might appear in both DP models
768  LauParameterPList& negSigDPPars = negSigModel_->getFloatingParameters();
769  LauParameterPList& posSigDPPars = posSigModel_->getFloatingParameters();
770  nSigDPPar_ += this->addResonanceParameters( negSigDPPars );
771  nSigDPPar_ += this->addResonanceParameters( posSigDPPars );
772 }
773 
774 void LauCPFitModel::setExtraPdfParameters()
775 {
776  // Include all the parameters of the PDF in the fit
777  // NB all of them are passed to the fit, even though some have been fixed through parameter.fixed(kTRUE)
778  // With the new "cloned parameter" scheme only "original" parameters are passed to the fit.
779  // Their clones are updated automatically when the originals are updated.
780 
781  nExtraPdfPar_ = 0;
782 
783  nExtraPdfPar_ += this->addFitParameters( negSignalPdfs_ );
784  if ( tagged_ ) {
785  nExtraPdfPar_ += this->addFitParameters( posSignalPdfs_ );
786  }
787 
788  if ( useSCF_ == kTRUE ) {
789  nExtraPdfPar_ += this->addFitParameters( negScfPdfs_ );
790  if ( tagged_ ) {
791  nExtraPdfPar_ += this->addFitParameters( posScfPdfs_ );
792  }
793  }
794 
795  if ( usingBkgnd_ == kTRUE ) {
796  for ( LauBkgndPdfsList::iterator iter = negBkgndPdfs_.begin(); iter != negBkgndPdfs_.end();
797  ++iter ) {
798  nExtraPdfPar_ += this->addFitParameters( *iter );
799  }
800  if ( tagged_ ) {
801  for ( LauBkgndPdfsList::iterator iter = posBkgndPdfs_.begin();
802  iter != posBkgndPdfs_.end();
803  ++iter ) {
804  nExtraPdfPar_ += this->addFitParameters( *iter );
805  }
806  }
807  }
808 }
809 
810 void LauCPFitModel::setFitNEvents()
811 {
812  if ( signalEvents_ == 0 ) {
813  std::cerr << "ERROR in LauCPFitModel::setFitNEvents : Signal yield not defined." << std::endl;
814  return;
815  }
816  nNormPar_ = 0;
817 
818  // initialise the total number of events to be the sum of all the hypotheses
819  Double_t nTotEvts = signalEvents_->value();
820  for ( LauBkgndYieldList::const_iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end();
821  ++iter ) {
822  nTotEvts += ( *iter )->value();
823  if ( ( *iter ) == 0 ) {
824  std::cerr << "ERROR in LauCPFitModel::setFitNEvents : Background yield not defined."
825  << std::endl;
826  return;
827  }
828  }
829  this->eventsPerExpt( TMath::FloorNint( nTotEvts ) );
830 
831  // if doing an extended ML fit add the number of signal events into the fit parameters
832  if ( this->doEMLFit() ) {
833  std::cout << "INFO in LauCPFitModel::setFitNEvents : Initialising number of events for signal and background components..."
834  << std::endl;
835  // add the signal fraction to the list of fit parameters
836  nNormPar_ += this->addFitParameters( signalEvents_ );
837  } else {
838  std::cout << "INFO in LauCPFitModel::setFitNEvents : Initialising number of events for background components (and hence signal)..."
839  << std::endl;
840  }
841 
842  // if not using the DP in the model we need an explicit signal asymmetry parameter
843  if ( this->useDP() == kFALSE ) {
844  nNormPar_ += this->addFitParameters( signalAsym_ );
845  }
846 
847  if ( useSCF_ && ! useSCFHist_ ) {
848  nNormPar_ += this->addFitParameters( &scfFrac_ );
849  }
850 
851  if ( usingBkgnd_ == kTRUE ) {
852  nNormPar_ += this->addFitParameters( bkgndEvents_ );
853  nNormPar_ += this->addFitParameters( bkgndAsym_ );
854  }
855 }
856 
857 void LauCPFitModel::setExtraNtupleVars()
858 {
859  // Set-up other parameters derived from the fit results, e.g. fit fractions.
860 
861  if ( this->useDP() != kTRUE ) {
862  return;
863  }
864 
865  // First clear the vectors so we start from scratch
866  this->clearExtraVarVectors();
867 
868  LauParameterList& extraVars = this->extraPars();
869 
870  // Add the positive and negative fit fractions for each signal component
871  negFitFrac_ = negSigModel_->getFitFractions();
872  if ( negFitFrac_.size() != nSigComp_ ) {
873  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "
874  << negFitFrac_.size() << std::endl;
875  gSystem->Exit( EXIT_FAILURE );
876  }
877  for ( UInt_t i( 0 ); i < nSigComp_; ++i ) {
878  if ( negFitFrac_[i].size() != nSigComp_ ) {
879  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "
880  << negFitFrac_[i].size() << std::endl;
881  gSystem->Exit( EXIT_FAILURE );
882  }
883  }
884 
885  posFitFrac_ = posSigModel_->getFitFractions();
886  if ( posFitFrac_.size() != nSigComp_ ) {
887  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "
888  << posFitFrac_.size() << std::endl;
889  gSystem->Exit( EXIT_FAILURE );
890  }
891  for ( UInt_t i( 0 ); i < nSigComp_; ++i ) {
892  if ( posFitFrac_[i].size() != nSigComp_ ) {
893  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "
894  << posFitFrac_[i].size() << std::endl;
895  gSystem->Exit( EXIT_FAILURE );
896  }
897  }
898 
899  // Add the positive and negative fit fractions that have not been corrected for the efficiency for each signal component
900  negFitFracEffUnCorr_ = negSigModel_->getFitFractionsEfficiencyUncorrected();
901  if ( negFitFracEffUnCorr_.size() != nSigComp_ ) {
902  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "
903  << negFitFracEffUnCorr_.size() << std::endl;
904  gSystem->Exit( EXIT_FAILURE );
905  }
906  for ( UInt_t i( 0 ); i < nSigComp_; ++i ) {
907  if ( negFitFracEffUnCorr_[i].size() != nSigComp_ ) {
908  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "
909  << negFitFracEffUnCorr_[i].size() << std::endl;
910  gSystem->Exit( EXIT_FAILURE );
911  }
912  }
913 
914  posFitFracEffUnCorr_ = posSigModel_->getFitFractionsEfficiencyUncorrected();
915  if ( posFitFracEffUnCorr_.size() != nSigComp_ ) {
916  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "
917  << posFitFracEffUnCorr_.size() << std::endl;
918  gSystem->Exit( EXIT_FAILURE );
919  }
920  for ( UInt_t i( 0 ); i < nSigComp_; ++i ) {
921  if ( posFitFracEffUnCorr_[i].size() != nSigComp_ ) {
922  std::cerr << "ERROR in LauCPFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "
923  << posFitFracEffUnCorr_[i].size() << std::endl;
924  gSystem->Exit( EXIT_FAILURE );
925  }
926  }
927 
928  for ( UInt_t i( 0 ); i < nSigComp_; ++i ) {
929  for ( UInt_t j( i ); j < nSigComp_; ++j ) {
930  TString name = negFitFrac_[i][j].name();
931  name.Insert( name.Index( "FitFrac" ), "Neg" );
932  negFitFrac_[i][j].name( name );
933  extraVars.push_back( negFitFrac_[i][j] );
934 
935  name = posFitFrac_[i][j].name();
936  name.Insert( name.Index( "FitFrac" ), "Pos" );
937  posFitFrac_[i][j].name( name );
938  extraVars.push_back( posFitFrac_[i][j] );
939 
940  name = negFitFracEffUnCorr_[i][j].name();
941  name.Insert( name.Index( "FitFrac" ), "Neg" );
942  negFitFracEffUnCorr_[i][j].name( name );
943  extraVars.push_back( negFitFracEffUnCorr_[i][j] );
944 
945  name = posFitFracEffUnCorr_[i][j].name();
946  name.Insert( name.Index( "FitFrac" ), "Pos" );
947  posFitFracEffUnCorr_[i][j].name( name );
948  extraVars.push_back( posFitFracEffUnCorr_[i][j] );
949  }
950  }
951 
952  // Store any extra parameters/quantities from the DP model (e.g. K-matrix total fit fractions)
953  std::vector<LauParameter> negExtraPars = negSigModel_->getExtraParameters();
954  std::vector<LauParameter>::iterator negExtraIter;
955  for ( negExtraIter = negExtraPars.begin(); negExtraIter != negExtraPars.end(); ++negExtraIter ) {
956  LauParameter negExtraParameter = ( *negExtraIter );
957  extraVars.push_back( negExtraParameter );
958  }
959 
960  std::vector<LauParameter> posExtraPars = posSigModel_->getExtraParameters();
961  std::vector<LauParameter>::iterator posExtraIter;
962  for ( posExtraIter = posExtraPars.begin(); posExtraIter != posExtraPars.end(); ++posExtraIter ) {
963  LauParameter posExtraParameter = ( *posExtraIter );
964  extraVars.push_back( posExtraParameter );
965  }
966 
967  // Now add in the DP efficiency value
968  Double_t initMeanEff = negSigModel_->getMeanEff().initValue();
969  negMeanEff_.value( initMeanEff );
970  negMeanEff_.genValue( initMeanEff );
971  negMeanEff_.initValue( initMeanEff );
972  extraVars.push_back( negMeanEff_ );
973 
974  initMeanEff = posSigModel_->getMeanEff().initValue();
975  posMeanEff_.value( initMeanEff );
976  posMeanEff_.genValue( initMeanEff );
977  posMeanEff_.initValue( initMeanEff );
978  extraVars.push_back( posMeanEff_ );
979 
980  // Also add in the DP rates
981  Double_t initDPRate = negSigModel_->getDPRate().initValue();
982  negDPRate_.value( initDPRate );
983  negDPRate_.genValue( initDPRate );
984  negDPRate_.initValue( initDPRate );
985  extraVars.push_back( negDPRate_ );
986 
987  initDPRate = posSigModel_->getDPRate().initValue();
988  posDPRate_.value( initDPRate );
989  posDPRate_.genValue( initDPRate );
990  posDPRate_.initValue( initDPRate );
991  extraVars.push_back( posDPRate_ );
992 
993  // Calculate the CPC and CPV Fit Fractions, ACPs and FitFrac asymmetries
994  this->calcExtraFractions( kTRUE );
995  this->calcAsymmetries( kTRUE );
996 
997  // Add the CP violating and CP conserving fit fractions for each signal component
998  for ( UInt_t i = 0; i < nSigComp_; i++ ) {
999  for ( UInt_t j = i; j < nSigComp_; j++ ) {
1000  extraVars.push_back( CPVFitFrac_[i][j] );
1001  }
1002  }
1003  for ( UInt_t i = 0; i < nSigComp_; i++ ) {
1004  for ( UInt_t j = i; j < nSigComp_; j++ ) {
1005  extraVars.push_back( CPCFitFrac_[i][j] );
1006  }
1007  }
1008 
1009  // Add the Fit Fraction asymmetry for each signal component
1010  for ( UInt_t i = 0; i < nSigComp_; i++ ) {
1011  extraVars.push_back( fitFracAsymm_[i] );
1012  }
1013 
1014  // Add the calculated CP asymmetry for each signal component
1015  for ( UInt_t i = 0; i < nSigComp_; i++ ) {
1016  extraVars.push_back( acp_[i] );
1017  }
1018 }
1019 
1020 void LauCPFitModel::calcExtraFractions( Bool_t initValues )
1021 {
1022  // Calculate the CP-conserving and CP-violating fit fractions
1023 
1024  if ( initValues ) {
1025  // create the structure
1026  CPCFitFrac_.clear();
1027  CPVFitFrac_.clear();
1028  CPCFitFrac_.resize( nSigComp_ );
1029  CPVFitFrac_.resize( nSigComp_ );
1030  for ( UInt_t i( 0 ); i < nSigComp_; ++i ) {
1031  CPCFitFrac_[i].resize( nSigComp_ );
1032  CPVFitFrac_[i].resize( nSigComp_ );
1033 
1034  for ( UInt_t j( i ); j < nSigComp_; ++j ) {
1035  TString name = negFitFrac_[i][j].name();
1036  name.Replace( name.Index( "Neg" ), 3, "CPC" );
1037  CPCFitFrac_[i][j].name( name );
1038  CPCFitFrac_[i][j].valueAndRange( 0.0, -100.0, 100.0 );
1039 
1040  name = negFitFrac_[i][j].name();
1041  name.Replace( name.Index( "Neg" ), 3, "CPV" );
1042  CPVFitFrac_[i][j].name( name );
1043  CPVFitFrac_[i][j].valueAndRange( 0.0, -100.0, 100.0 );
1044  }
1045  }
1046  }
1047 
1048  Double_t denom = negDPRate_.value() + posDPRate_.value();
1049 
1050  for ( UInt_t i( 0 ); i < nSigComp_; ++i ) {
1051  for ( UInt_t j( i ); j < nSigComp_; ++j ) {
1052 
1053  Double_t negTerm = negFitFrac_[i][j].value() * negDPRate_.value();
1054  Double_t posTerm = posFitFrac_[i][j].value() * posDPRate_.value();
1055 
1056  Double_t cpcFitFrac = ( negTerm + posTerm ) / denom;
1057  Double_t cpvFitFrac = ( negTerm - posTerm ) / denom;
1058 
1059  CPCFitFrac_[i][j].value( cpcFitFrac );
1060  CPVFitFrac_[i][j].value( cpvFitFrac );
1061 
1062  if ( initValues ) {
1063  CPCFitFrac_[i][j].genValue( cpcFitFrac );
1064  CPCFitFrac_[i][j].initValue( cpcFitFrac );
1065 
1066  CPVFitFrac_[i][j].genValue( cpvFitFrac );
1067  CPVFitFrac_[i][j].initValue( cpvFitFrac );
1068  }
1069  }
1070  }
1071 }
1072 
1073 void LauCPFitModel::calcAsymmetries( Bool_t initValues )
1074 {
1075  // Calculate the CP asymmetries
1076  // Also calculate the fit fraction asymmetries
1077 
1078  for ( UInt_t i = 0; i < nSigComp_; i++ ) {
1079 
1080  acp_[i] = coeffPars_[i]->acp();
1081 
1082  LauAsymmCalc asymmCalc( negFitFrac_[i][i].value(), posFitFrac_[i][i].value() );
1083  Double_t asym = asymmCalc.getAsymmetry();
1084  fitFracAsymm_[i].value( asym );
1085  if ( initValues ) {
1086  fitFracAsymm_[i].genValue( asym );
1087  fitFracAsymm_[i].initValue( asym );
1088  }
1089  }
1090 }
1091 
1092 void LauCPFitModel::finaliseFitResults( const TString& tablePrefixName )
1093 {
1094  // Retrieve parameters from the fit results for calculations and toy generation
1095  // and eventually store these in output root ntuples/text files
1096 
1097  // Now take the fit parameters and update them as necessary
1098  // i.e. to make mag > 0.0, phase in the right range.
1099  // This function will also calculate any other values, such as the
1100  // fit fractions, using any errors provided by fitParErrors as appropriate.
1101  // Also obtain the pull values: (measured - generated)/(average error)
1102 
1103  if ( this->useDP() == kTRUE ) {
1104  for ( UInt_t i = 0; i < nSigComp_; ++i ) {
1105  // Check whether we have "a/b > 0.0", and phases in the right range
1106  coeffPars_[i]->finaliseValues();
1107  }
1108  }
1109 
1110  // update the pulls on the event fractions and asymmetries
1111  if ( this->doEMLFit() ) {
1112  signalEvents_->updatePull();
1113  }
1114  if ( this->useDP() == kFALSE ) {
1115  signalAsym_->updatePull();
1116  }
1117  if ( useSCF_ && ! useSCFHist_ ) {
1118  scfFrac_.updatePull();
1119  }
1120  if ( usingBkgnd_ == kTRUE ) {
1121  for ( LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end();
1122  ++iter ) {
1123  std::vector<LauParameter*> parameters = ( *iter )->getPars();
1124  for ( LauParameter* parameter : parameters ) {
1125  parameter->updatePull();
1126  }
1127  }
1128  for ( LauBkgndYieldList::iterator iter = bkgndAsym_.begin(); iter != bkgndAsym_.end();
1129  ++iter ) {
1130  std::vector<LauParameter*> parameters = ( *iter )->getPars();
1131  for ( LauParameter* parameter : parameters ) {
1132  parameter->updatePull();
1133  }
1134  }
1135  }
1136 
1137  // Update the pulls on all the extra PDFs' parameters
1138  this->updateFitParameters( negSignalPdfs_ );
1139  this->updateFitParameters( posSignalPdfs_ );
1140  if ( useSCF_ == kTRUE ) {
1141  this->updateFitParameters( negScfPdfs_ );
1142  this->updateFitParameters( posScfPdfs_ );
1143  }
1144  if ( usingBkgnd_ == kTRUE ) {
1145  for ( LauBkgndPdfsList::iterator iter = negBkgndPdfs_.begin(); iter != negBkgndPdfs_.end();
1146  ++iter ) {
1147  this->updateFitParameters( *iter );
1148  }
1149  for ( LauBkgndPdfsList::iterator iter = posBkgndPdfs_.begin(); iter != posBkgndPdfs_.end();
1150  ++iter ) {
1151  this->updateFitParameters( *iter );
1152  }
1153  }
1154 
1155  // Fill the fit results to the ntuple
1156 
1157  // update the coefficients and then calculate the fit fractions and ACP's
1158  if ( this->useDP() == kTRUE ) {
1159  this->updateCoeffs();
1160  negSigModel_->updateCoeffs( negCoeffs_ );
1161  negSigModel_->calcExtraInfo();
1162  posSigModel_->updateCoeffs( posCoeffs_ );
1163  posSigModel_->calcExtraInfo();
1164 
1165  LauParArray negFitFrac = negSigModel_->getFitFractions();
1166  if ( negFitFrac.size() != nSigComp_ ) {
1167  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "
1168  << negFitFrac.size() << std::endl;
1169  gSystem->Exit( EXIT_FAILURE );
1170  }
1171  for ( UInt_t i( 0 ); i < nSigComp_; ++i ) {
1172  if ( negFitFrac[i].size() != nSigComp_ ) {
1173  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "
1174  << negFitFrac[i].size() << std::endl;
1175  gSystem->Exit( EXIT_FAILURE );
1176  }
1177  }
1178 
1179  LauParArray posFitFrac = posSigModel_->getFitFractions();
1180  if ( posFitFrac.size() != nSigComp_ ) {
1181  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "
1182  << posFitFrac.size() << std::endl;
1183  gSystem->Exit( EXIT_FAILURE );
1184  }
1185  for ( UInt_t i( 0 ); i < nSigComp_; ++i ) {
1186  if ( posFitFrac[i].size() != nSigComp_ ) {
1187  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "
1188  << posFitFrac[i].size() << std::endl;
1189  gSystem->Exit( EXIT_FAILURE );
1190  }
1191  }
1192 
1193  LauParArray negFitFracEffUnCorr = negSigModel_->getFitFractionsEfficiencyUncorrected();
1194  if ( negFitFracEffUnCorr.size() != nSigComp_ ) {
1195  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "
1196  << negFitFracEffUnCorr.size() << std::endl;
1197  gSystem->Exit( EXIT_FAILURE );
1198  }
1199  for ( UInt_t i( 0 ); i < nSigComp_; ++i ) {
1200  if ( negFitFracEffUnCorr[i].size() != nSigComp_ ) {
1201  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "
1202  << negFitFracEffUnCorr[i].size() << std::endl;
1203  gSystem->Exit( EXIT_FAILURE );
1204  }
1205  }
1206 
1207  LauParArray posFitFracEffUnCorr = posSigModel_->getFitFractionsEfficiencyUncorrected();
1208  if ( posFitFracEffUnCorr.size() != nSigComp_ ) {
1209  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "
1210  << posFitFracEffUnCorr.size() << std::endl;
1211  gSystem->Exit( EXIT_FAILURE );
1212  }
1213  for ( UInt_t i( 0 ); i < nSigComp_; ++i ) {
1214  if ( posFitFracEffUnCorr[i].size() != nSigComp_ ) {
1215  std::cerr << "ERROR in LauCPFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "
1216  << posFitFracEffUnCorr[i].size() << std::endl;
1217  gSystem->Exit( EXIT_FAILURE );
1218  }
1219  }
1220 
1221  for ( UInt_t i( 0 ); i < nSigComp_; ++i ) {
1222  for ( UInt_t j( i ); j < nSigComp_; ++j ) {
1223  negFitFrac_[i][j].value( negFitFrac[i][j].value() );
1224  posFitFrac_[i][j].value( posFitFrac[i][j].value() );
1225  negFitFracEffUnCorr_[i][j].value( negFitFracEffUnCorr[i][j].value() );
1226  posFitFracEffUnCorr_[i][j].value( posFitFracEffUnCorr[i][j].value() );
1227  }
1228  }
1229 
1230  negMeanEff_.value( negSigModel_->getMeanEff().value() );
1231  posMeanEff_.value( posSigModel_->getMeanEff().value() );
1232  negDPRate_.value( negSigModel_->getDPRate().value() );
1233  posDPRate_.value( posSigModel_->getDPRate().value() );
1234 
1235  this->calcExtraFractions();
1236  this->calcAsymmetries();
1237 
1238  // Then store the final fit parameters, and any extra parameters for
1239  // the signal model (e.g. fit fractions, FF asymmetries, ACPs, mean efficiency and DP rate)
1240 
1241  this->clearExtraVarVectors();
1242  LauParameterList& extraVars = this->extraPars();
1243 
1244  // Add the positive and negative fit fractions for each signal component
1245  for ( UInt_t i( 0 ); i < nSigComp_; ++i ) {
1246  for ( UInt_t j( i ); j < nSigComp_; ++j ) {
1247  extraVars.push_back( negFitFrac_[i][j] );
1248  extraVars.push_back( posFitFrac_[i][j] );
1249  extraVars.push_back( negFitFracEffUnCorr_[i][j] );
1250  extraVars.push_back( posFitFracEffUnCorr_[i][j] );
1251  }
1252  }
1253 
1254  // Store any extra parameters/quantities from the DP model (e.g. K-matrix total fit fractions)
1255  std::vector<LauParameter> negExtraPars = negSigModel_->getExtraParameters();
1256  std::vector<LauParameter>::iterator negExtraIter;
1257  for ( negExtraIter = negExtraPars.begin(); negExtraIter != negExtraPars.end();
1258  ++negExtraIter ) {
1259  LauParameter negExtraParameter = ( *negExtraIter );
1260  extraVars.push_back( negExtraParameter );
1261  }
1262 
1263  std::vector<LauParameter> posExtraPars = posSigModel_->getExtraParameters();
1264  std::vector<LauParameter>::iterator posExtraIter;
1265  for ( posExtraIter = posExtraPars.begin(); posExtraIter != posExtraPars.end();
1266  ++posExtraIter ) {
1267  LauParameter posExtraParameter = ( *posExtraIter );
1268  extraVars.push_back( posExtraParameter );
1269  }
1270 
1271  extraVars.push_back( negMeanEff_ );
1272  extraVars.push_back( posMeanEff_ );
1273  extraVars.push_back( negDPRate_ );
1274  extraVars.push_back( posDPRate_ );
1275 
1276  for ( UInt_t i = 0; i < nSigComp_; i++ ) {
1277  for ( UInt_t j( i ); j < nSigComp_; ++j ) {
1278  extraVars.push_back( CPVFitFrac_[i][j] );
1279  }
1280  }
1281  for ( UInt_t i = 0; i < nSigComp_; i++ ) {
1282  for ( UInt_t j( i ); j < nSigComp_; ++j ) {
1283  extraVars.push_back( CPCFitFrac_[i][j] );
1284  }
1285  }
1286  for ( UInt_t i( 0 ); i < nSigComp_; ++i ) {
1287  extraVars.push_back( fitFracAsymm_[i] );
1288  }
1289  for ( UInt_t i( 0 ); i < nSigComp_; ++i ) {
1290  extraVars.push_back( acp_[i] );
1291  }
1292 
1293  this->printFitFractions( std::cout );
1294  this->printAsymmetries( std::cout );
1295  }
1296 
1297  LauFitNtuple* ntuple = this->fitNtuple();
1298  ntuple->storeParsAndErrors( this->fitPars(), this->multiDimConstrainedPars(), this->extraPars() );
1299 
1300  // find out the correlation matrix for the parameters
1301  ntuple->storeCorrMatrix( this->iExpt(), this->fitStatus(), this->covarianceMatrix() );
1302 
1303  // Fill the data into ntuple
1304  ntuple->updateFitNtuple();
1305 
1306  // Print out the partial fit fractions, phases and the
1307  // averaged efficiency, reweighted by the dynamics (and anything else)
1308  if ( this->writeLatexTable() ) {
1309  TString sigOutFileName( tablePrefixName );
1310  sigOutFileName += "_";
1311  sigOutFileName += this->iExpt();
1312  sigOutFileName += "Expt.tex";
1313  this->writeOutTable( sigOutFileName );
1314  }
1315 }
1316 
1317 void LauCPFitModel::printFitFractions( std::ostream& output )
1318 {
1319  // Print out Fit Fractions, total DP rate and mean efficiency
1320  // First for the B- events
1321  for ( UInt_t i = 0; i < nSigComp_; i++ ) {
1322  const TString compName( coeffPars_[i]->name() );
1323  output << negParent_ << " FitFraction for component " << i << " (" << compName
1324  << ") = " << negFitFrac_[i][i] << std::endl;
1325  }
1326  output << negParent_ << " overall DP rate (integral of matrix element squared) = " << negDPRate_
1327  << std::endl;
1328  output << negParent_ << " average efficiency weighted by whole DP dynamics = " << negMeanEff_
1329  << std::endl;
1330 
1331  // Then for the positive sample
1332  for ( UInt_t i = 0; i < nSigComp_; i++ ) {
1333  const TString compName( coeffPars_[i]->name() );
1334  const TString conjName( negSigModel_->getConjResName( compName ) );
1335  output << posParent_ << " FitFraction for component " << i << " (" << conjName
1336  << ") = " << posFitFrac_[i][i] << std::endl;
1337  }
1338  output << posParent_ << " overall DP rate (integral of matrix element squared) = " << posDPRate_
1339  << std::endl;
1340  output << posParent_ << " average efficiency weighted by whole DP dynamics = " << posMeanEff_
1341  << std::endl;
1342 }
1343 
1344 void LauCPFitModel::printAsymmetries( std::ostream& output )
1345 {
1346  for ( UInt_t i = 0; i < nSigComp_; i++ ) {
1347  const TString compName( coeffPars_[i]->name() );
1348  output << "Fit Fraction asymmetry for component " << i << " (" << compName
1349  << ") = " << fitFracAsymm_[i] << std::endl;
1350  }
1351  for ( UInt_t i = 0; i < nSigComp_; i++ ) {
1352  const TString compName( coeffPars_[i]->name() );
1353  output << "ACP for component " << i << " (" << compName << ") = " << acp_[i] << std::endl;
1354  }
1355 }
1356 
1357 void LauCPFitModel::writeOutTable( const TString& outputFile )
1358 {
1359  // Write out the results of the fit to a tex-readable table
1360  // TODO - need to include the yields in this table
1361  std::ofstream fout( outputFile );
1362  LauPrint print;
1363 
1364  std::cout << "INFO in LauCPFitModel::writeOutTable : Writing out results of the fit to the tex file "
1365  << outputFile << std::endl;
1366 
1367  if ( this->useDP() == kTRUE ) {
1368  // print the fit coefficients in one table
1369  coeffPars_.front()->printTableHeading( fout );
1370  for ( UInt_t i = 0; i < nSigComp_; i++ ) {
1371  coeffPars_[i]->printTableRow( fout );
1372  }
1373  fout << "\\hline" << std::endl;
1374  fout << "\\end{tabular}" << std::endl << std::endl;
1375 
1376  // print the fit fractions and asymmetries in another
1377  fout << "\\begin{tabular}{|l|c|c|c|c|}" << std::endl;
1378  fout << "\\hline" << std::endl;
1379  fout << "Component & " << negParent_ << " Fit Fraction & " << posParent_
1380  << " Fit Fraction & Fit Fraction Asymmetry & ACP \\\\" << std::endl;
1381  fout << "\\hline" << std::endl;
1382 
1383  Double_t negFitFracSum( 0.0 );
1384  Double_t posFitFracSum( 0.0 );
1385 
1386  for ( UInt_t i = 0; i < nSigComp_; i++ ) {
1387  TString resName = coeffPars_[i]->name();
1388  resName = resName.ReplaceAll( "_", "\\_" );
1389 
1390  Double_t negFitFrac = negFitFrac_[i][i].value();
1391  Double_t posFitFrac = posFitFrac_[i][i].value();
1392  negFitFracSum += negFitFrac;
1393  posFitFracSum += posFitFrac;
1394 
1395  Double_t fitFracAsymm = fitFracAsymm_[i].value();
1396 
1397  Double_t acp = acp_[i].value();
1398  Double_t acpErr = acp_[i].error();
1399 
1400  fout << resName << " & $";
1401 
1402  print.printFormat( fout, negFitFrac );
1403  fout << "$ & $";
1404  print.printFormat( fout, posFitFrac );
1405  fout << "$ & $";
1406 
1407  print.printFormat( fout, fitFracAsymm );
1408  fout << "$ & $";
1409 
1410  print.printFormat( fout, acp );
1411  fout << " \\pm ";
1412  print.printFormat( fout, acpErr );
1413  fout << "$ \\\\" << std::endl;
1414  }
1415  fout << "\\hline" << std::endl;
1416 
1417  // Also print out sum of fit fractions
1418  fout << "Fit Fraction Sum & $";
1419  print.printFormat( fout, negFitFracSum );
1420  fout << "$ & $";
1421  print.printFormat( fout, posFitFracSum );
1422  fout << "$ & & \\\\" << std::endl;
1423  fout << "\\hline" << std::endl;
1424 
1425  fout << "DP rate & $";
1426  print.printFormat( fout, negDPRate_.value() );
1427  fout << "$ & $";
1428  print.printFormat( fout, posDPRate_.value() );
1429  fout << "$ & & \\\\" << std::endl;
1430 
1431  fout << "$< \\varepsilon > $ & $";
1432  print.printFormat( fout, negMeanEff_.value() );
1433  fout << "$ & $";
1434  print.printFormat( fout, posMeanEff_.value() );
1435  fout << "$ & & \\\\" << std::endl;
1436  fout << "\\hline" << std::endl;
1437  fout << "\\end{tabular}" << std::endl << std::endl;
1438  }
1439 
1440  if ( ! negSignalPdfs_.empty() ) {
1441  fout << "\\begin{tabular}{|l|c|}" << std::endl;
1442  fout << "\\hline" << std::endl;
1443  if ( useSCF_ == kTRUE ) {
1444  fout << "\\Extra TM Signal PDFs' Parameters: & \\\\" << std::endl;
1445  } else {
1446  fout << "\\Extra Signal PDFs' Parameters: & \\\\" << std::endl;
1447  }
1448  this->printFitParameters( negSignalPdfs_, fout );
1449  if ( tagged_ ) {
1450  this->printFitParameters( posSignalPdfs_, fout );
1451  }
1452  if ( useSCF_ == kTRUE && ! negScfPdfs_.empty() ) {
1453  fout << "\\hline" << std::endl;
1454  fout << "\\Extra SCF Signal PDFs' Parameters: & \\\\" << std::endl;
1455  this->printFitParameters( negScfPdfs_, fout );
1456  if ( tagged_ ) {
1457  this->printFitParameters( posScfPdfs_, fout );
1458  }
1459  }
1460  if ( usingBkgnd_ == kTRUE && ! negBkgndPdfs_.empty() ) {
1461  fout << "\\hline" << std::endl;
1462  fout << "\\Extra Background PDFs' Parameters: & \\\\" << std::endl;
1463  for ( LauBkgndPdfsList::const_iterator iter = negBkgndPdfs_.begin();
1464  iter != negBkgndPdfs_.end();
1465  ++iter ) {
1466  this->printFitParameters( *iter, fout );
1467  }
1468  if ( tagged_ ) {
1469  for ( LauBkgndPdfsList::const_iterator iter = posBkgndPdfs_.begin();
1470  iter != posBkgndPdfs_.end();
1471  ++iter ) {
1472  this->printFitParameters( *iter, fout );
1473  }
1474  }
1475  }
1476  fout << "\\hline \n\\end{tabular}" << std::endl << std::endl;
1477  }
1478 }
1479 
1480 void LauCPFitModel::checkInitFitParams()
1481 {
1482  // Update the number of signal events to be total-sum(background events)
1483  this->updateSigEvents();
1484 
1485  // Check whether we want to have randomised initial fit parameters for the signal model
1486  if ( this->useRandomInitFitPars() == kTRUE ) {
1487  std::cout << "INFO in LauCPFitModel::checkInitFitParams : Setting random parameters for the signal model"
1488  << std::endl;
1489  this->randomiseInitFitPars();
1490  }
1491 }
1492 
1493 void LauCPFitModel::randomiseInitFitPars()
1494 {
1495  // Only randomise those parameters that are not fixed!
1496  std::cout << "INFO in LauCPFitModel::randomiseInitFitPars : Randomising the initial fit magnitudes and phases of the components..."
1497  << std::endl;
1498 
1499  if ( fixParamFileName_.IsNull() && fixParamMap_.empty() ) {
1500  // No params are imported - randomise as normal
1501  for ( UInt_t i = 0; i < nSigComp_; i++ ) {
1502  coeffPars_[i]->randomiseInitValues();
1503  }
1504  } else {
1505  // Only randomise those that are not imported (i.e., not found in allImportedFreeParams_)
1506  // by temporarily fixing all imported parameters, and then freeing those not set to be fixed when imported,
1507  // except those that are previously set to be fixed anyhow.
1508  // Convoluted, but beats changing the behaviour of functions that call checkInitFitParams or the coeffSet
1509  // itself.
1510 
1511  for ( auto p : allImportedFreeParams_ ) {
1512  p->fixed( kTRUE );
1513  }
1514 
1515  for ( UInt_t i = 0; i < nSigComp_; i++ ) {
1516  coeffPars_[i]->randomiseInitValues();
1517  }
1518 
1519  for ( auto p : allImportedFreeParams_ ) {
1520  p->fixed( kFALSE );
1521  }
1522  }
1523 }
1524 
1525 std::pair<LauCPFitModel::LauGenInfo, Bool_t> LauCPFitModel::eventsToGenerate()
1526 {
1527  // Determine the number of events to generate for each hypothesis
1528  // If we're smearing then smear each one individually
1529 
1530  LauGenInfo nEvtsGen;
1531 
1532  // Keep track of whether any yield or asymmetry parameters are blinded
1533  Bool_t blind = kFALSE;
1534 
1535  // Signal
1536  if ( signalEvents_->blind() ) {
1537  blind = kTRUE;
1538  }
1539 
1540  Double_t evtWeight( 1.0 );
1541  Double_t nEvts = signalEvents_->genValue();
1542  if ( nEvts < 0.0 ) {
1543  evtWeight = -1.0;
1544  nEvts = TMath::Abs( nEvts );
1545  }
1546 
1547  Double_t sigAsym( 0.0 );
1548  // need to include this as an alternative in case the DP isn't in the model
1549  if ( ! this->useDP() || forceAsym_ ) {
1550  sigAsym = signalAsym_->genValue();
1551  if ( signalAsym_->blind() ) {
1552  blind = kTRUE;
1553  }
1554  } else {
1555  Double_t negRate = negSigModel_->getDPNorm();
1556  Double_t posRate = posSigModel_->getDPNorm();
1557  if ( negRate + posRate > 1e-30 ) {
1558  sigAsym = ( negRate - posRate ) / ( negRate + posRate );
1559  }
1560  }
1561  Double_t nPosEvts = ( nEvts / 2.0 * ( 1.0 - sigAsym ) );
1562  Double_t nNegEvts = ( nEvts / 2.0 * ( 1.0 + sigAsym ) );
1563 
1564  Int_t nPosEvtsToGen { static_cast<Int_t>( nPosEvts ) };
1565  Int_t nNegEvtsToGen { static_cast<Int_t>( nNegEvts ) };
1566  if ( this->doPoissonSmearing() ) {
1567  nPosEvtsToGen = LauRandom::randomFun()->Poisson( nPosEvts );
1568  nNegEvtsToGen = LauRandom::randomFun()->Poisson( nNegEvts );
1569  }
1570 
1571  nEvtsGen[std::make_pair( "signal", +1 )] = std::make_pair( nPosEvtsToGen, evtWeight );
1572  nEvtsGen[std::make_pair( "signal", -1 )] = std::make_pair( nNegEvtsToGen, evtWeight );
1573 
1574  // backgrounds
1575  const UInt_t nBkgnds = this->nBkgndClasses();
1576  for ( UInt_t bkgndID( 0 ); bkgndID < nBkgnds; ++bkgndID ) {
1577  const LauAbsRValue* evtsPar = bkgndEvents_[bkgndID];
1578  const LauAbsRValue* asymPar = bkgndAsym_[bkgndID];
1579  if ( evtsPar->blind() || asymPar->blind() ) {
1580  blind = kTRUE;
1581  }
1582 
1583  evtWeight = 1.0;
1584  nEvts = evtsPar->genValue();
1585  if ( nEvts < 0 ) {
1586  evtWeight = -1.0;
1587  nEvts = TMath::Abs( nEvts );
1588  }
1589 
1590  const Double_t asym = asymPar->genValue();
1591  nPosEvts = ( nEvts / 2.0 * ( 1.0 - asym ) );
1592  nNegEvts = ( nEvts / 2.0 * ( 1.0 + asym ) );
1593 
1594  nPosEvtsToGen = static_cast<Int_t>( nPosEvts );
1595  nNegEvtsToGen = static_cast<Int_t>( nNegEvts );
1596  if ( this->doPoissonSmearing() ) {
1597  nPosEvtsToGen = LauRandom::randomFun()->Poisson( nPosEvts );
1598  nNegEvtsToGen = LauRandom::randomFun()->Poisson( nNegEvts );
1599  }
1600 
1601  const TString& bkgndClass = this->bkgndClassName( bkgndID );
1602  nEvtsGen[std::make_pair( bkgndClass, +1 )] = std::make_pair( nPosEvtsToGen, evtWeight );
1603  nEvtsGen[std::make_pair( bkgndClass, -1 )] = std::make_pair( nNegEvtsToGen, evtWeight );
1604  }
1605 
1606  // Print out the information on what we're generating, but only if none of the parameters are blind (otherwise we risk unblinding them!)
1607  if ( ! blind ) {
1608  std::cout << "INFO in LauCPFitModel::eventsToGenerate : Generating toy MC with:" << std::endl;
1609  std::cout << " : Signal asymmetry = " << sigAsym
1610  << " and number of signal events = " << signalEvents_->genValue() << std::endl;
1611  for ( UInt_t bkgndID( 0 ); bkgndID < nBkgnds; ++bkgndID ) {
1612  const TString& bkgndClass = this->bkgndClassName( bkgndID );
1613  const LauAbsRValue* evtsPar = bkgndEvents_[bkgndID];
1614  const LauAbsRValue* asymPar = bkgndAsym_[bkgndID];
1615  std::cout << " : " << bkgndClass
1616  << " asymmetry = " << asymPar->genValue() << " and number of " << bkgndClass
1617  << " events = " << evtsPar->genValue() << std::endl;
1618  }
1619  }
1620 
1621  return std::make_pair( nEvtsGen, blind );
1622 }
1623 
1624 Bool_t LauCPFitModel::genExpt()
1625 {
1626  // Routine to generate toy Monte Carlo events according to the various models we have defined.
1627 
1628  // Determine the number of events to generate for each hypothesis
1629  std::pair<LauGenInfo, Bool_t> info = this->eventsToGenerate();
1630  LauGenInfo nEvts = info.first;
1631  const Bool_t blind = info.second;
1632 
1633  Bool_t genOK( kTRUE );
1634  Int_t evtNum( 0 );
1635 
1636  const UInt_t nBkgnds = this->nBkgndClasses();
1637  std::vector<TString> bkgndClassNames( nBkgnds );
1638  std::vector<TString> bkgndClassNamesGen( nBkgnds );
1639  for ( UInt_t iBkgnd( 0 ); iBkgnd < nBkgnds; ++iBkgnd ) {
1640  TString name( this->bkgndClassName( iBkgnd ) );
1641  bkgndClassNames[iBkgnd] = name;
1642  bkgndClassNamesGen[iBkgnd] = "gen" + name;
1643  }
1644 
1645  const Bool_t storeSCFTruthInfo = ( useSCF_ || ( this->enableEmbedding() && negSignalTree_ != 0 &&
1646  negSignalTree_->haveBranch( "mcMatch" ) &&
1647  posSignalTree_ != 0 &&
1648  posSignalTree_->haveBranch( "mcMatch" ) ) );
1649 
1650  // Loop over the hypotheses and generate the requested number of events for each one
1651  for ( LauGenInfo::const_iterator iter = nEvts.begin(); iter != nEvts.end(); ++iter ) {
1652 
1653  const TString& type( iter->first.first );
1654  curEvtCharge_ = iter->first.second;
1655  Double_t evtWeight( iter->second.second );
1656  Int_t nEvtsGen( iter->second.first );
1657 
1658  for ( Int_t iEvt( 0 ); iEvt < nEvtsGen; ++iEvt ) {
1659 
1660  this->setGenNtupleDoubleBranchValue( "evtWeight", evtWeight );
1661  this->setGenNtupleDoubleBranchValue( "efficiency", 1.0 );
1662 
1663  if ( type == "signal" ) {
1664  this->setGenNtupleIntegerBranchValue( "genSig", 1 );
1665  for ( UInt_t iBkgnd( 0 ); iBkgnd < nBkgnds; ++iBkgnd ) {
1666  this->setGenNtupleIntegerBranchValue( bkgndClassNamesGen[iBkgnd], 0 );
1667  }
1668  genOK = this->generateSignalEvent();
1669  if ( curEvtCharge_ > 0 ) {
1670  this->setGenNtupleDoubleBranchValue( "efficiency", posSigModel_->getEvtEff() );
1671  } else {
1672  this->setGenNtupleDoubleBranchValue( "efficiency", negSigModel_->getEvtEff() );
1673  }
1674 
1675  } else {
1676  this->setGenNtupleIntegerBranchValue( "genSig", 0 );
1677  if ( storeSCFTruthInfo ) {
1678  this->setGenNtupleIntegerBranchValue( "genTMSig", 0 );
1679  this->setGenNtupleIntegerBranchValue( "genSCFSig", 0 );
1680  }
1681  UInt_t bkgndID( 0 );
1682  for ( UInt_t iBkgnd( 0 ); iBkgnd < nBkgnds; ++iBkgnd ) {
1683  Int_t gen( 0 );
1684  if ( bkgndClassNames[iBkgnd] == type ) {
1685  gen = 1;
1686  bkgndID = iBkgnd;
1687  }
1688  this->setGenNtupleIntegerBranchValue( bkgndClassNamesGen[iBkgnd], gen );
1689  }
1690  genOK = this->generateBkgndEvent( bkgndID );
1691  }
1692  if ( ! genOK ) {
1693  // If there was a problem with the generation then break out and return.
1694  // The problem model will have adjusted itself so that all should be OK next time.
1695  break;
1696  }
1697  if ( this->useDP() == kTRUE ) {
1698  this->setDPBranchValues();
1699  }
1700 
1701  // Store the event charge
1702  this->setGenNtupleIntegerBranchValue( tagVarName_, curEvtCharge_ );
1703 
1704  // Store the event number (within this experiment)
1705  // and then increment it
1706  this->setGenNtupleIntegerBranchValue( "iEvtWithinExpt", evtNum );
1707  ++evtNum;
1708 
1709  this->fillGenNtupleBranches();
1710  if ( ! blind && ( iEvt % 500 == 0 ) ) {
1711  std::cout << "INFO in LauCPFitModel::genExpt : Generated event number " << iEvt
1712  << " out of " << nEvtsGen << " " << type << " events." << std::endl;
1713  }
1714  }
1715 
1716  if ( ! genOK ) {
1717  break;
1718  }
1719  }
1720 
1721  if ( this->useDP() && genOK ) {
1722 
1723  negSigModel_->checkToyMC( kTRUE, kTRUE );
1724  posSigModel_->checkToyMC( kTRUE, kTRUE );
1725 
1726  // Get the fit fractions if they're to be written into the latex table
1727  if ( this->writeLatexTable() ) {
1728  LauParArray negFitFrac = negSigModel_->getFitFractions();
1729  if ( negFitFrac.size() != nSigComp_ ) {
1730  std::cerr << "ERROR in LauCPFitModel::genExpt : Fit Fraction array of unexpected dimension: "
1731  << negFitFrac.size() << std::endl;
1732  gSystem->Exit( EXIT_FAILURE );
1733  }
1734  for ( UInt_t i( 0 ); i < nSigComp_; ++i ) {
1735  if ( negFitFrac[i].size() != nSigComp_ ) {
1736  std::cerr << "ERROR in LauCPFitModel::genExpt : Fit Fraction array of unexpected dimension: "
1737  << negFitFrac[i].size() << std::endl;
1738  gSystem->Exit( EXIT_FAILURE );
1739  }
1740  }
1741  LauParArray posFitFrac = posSigModel_->getFitFractions();
1742  if ( posFitFrac.size() != nSigComp_ ) {
1743  std::cerr << "ERROR in LauCPFitModel::genExpt : Fit Fraction array of unexpected dimension: "
1744  << posFitFrac.size() << std::endl;
1745  gSystem->Exit( EXIT_FAILURE );
1746  }
1747  for ( UInt_t i( 0 ); i < nSigComp_; ++i ) {
1748  if ( posFitFrac[i].size() != nSigComp_ ) {
1749  std::cerr << "ERROR in LauCPFitModel::genExpt : Fit Fraction array of unexpected dimension: "
1750  << posFitFrac[i].size() << std::endl;
1751  gSystem->Exit( EXIT_FAILURE );
1752  }
1753  }
1754 
1755  for ( UInt_t i( 0 ); i < nSigComp_; ++i ) {
1756  for ( UInt_t j( i ); j < nSigComp_; ++j ) {
1757  negFitFrac_[i][j].value( negFitFrac[i][j].value() );
1758  posFitFrac_[i][j].value( posFitFrac[i][j].value() );
1759  }
1760  }
1761  negMeanEff_.value( negSigModel_->getMeanEff().value() );
1762  posMeanEff_.value( posSigModel_->getMeanEff().value() );
1763  negDPRate_.value( negSigModel_->getDPRate().value() );
1764  posDPRate_.value( posSigModel_->getDPRate().value() );
1765  }
1766  }
1767 
1768  // If we're reusing embedded events or if the generation is being
1769  // reset then clear the lists of used events
1770  if ( reuseSignal_ || ! genOK ) {
1771  if ( negSignalTree_ ) {
1772  negSignalTree_->clearUsedList();
1773  }
1774  if ( posSignalTree_ ) {
1775  posSignalTree_->clearUsedList();
1776  }
1777  }
1778  for ( UInt_t bkgndID( 0 ); bkgndID < nBkgnds; ++bkgndID ) {
1779  LauEmbeddedData* data = negBkgndTree_[bkgndID];
1780  if ( reuseBkgnd_[bkgndID] || ! genOK ) {
1781  if ( data ) {
1782  data->clearUsedList();
1783  }
1784  }
1785  }
1786  for ( UInt_t bkgndID( 0 ); bkgndID < nBkgnds; ++bkgndID ) {
1787  LauEmbeddedData* data = posBkgndTree_[bkgndID];
1788  if ( reuseBkgnd_[bkgndID] || ! genOK ) {
1789  if ( data ) {
1790  data->clearUsedList();
1791  }
1792  }
1793  }
1794 
1795  return genOK;
1796 }
1797 
1798 Bool_t LauCPFitModel::generateSignalEvent()
1799 {
1800  // Generate signal event
1801  Bool_t genOK( kTRUE );
1802  Bool_t genSCF( kFALSE );
1803 
1804  LauIsobarDynamics* model( 0 );
1805  LauKinematics* kinematics( 0 );
1806  LauEmbeddedData* embeddedData( 0 );
1807  LauPdfPList* sigPdfs( 0 );
1808  LauPdfPList* scfPdfs( 0 );
1809 
1810  Bool_t doReweighting( kFALSE );
1811 
1812  if ( curEvtCharge_ < 0 ) {
1813  model = negSigModel_;
1814  kinematics = negKinematics_;
1815  sigPdfs = &negSignalPdfs_;
1816  scfPdfs = &negScfPdfs_;
1817  if ( this->enableEmbedding() ) {
1818  embeddedData = negSignalTree_;
1819  doReweighting = useNegReweighting_;
1820  }
1821  } else {
1822  model = posSigModel_;
1823  kinematics = posKinematics_;
1824  if ( tagged_ ) {
1825  sigPdfs = &posSignalPdfs_;
1826  scfPdfs = &posScfPdfs_;
1827  } else {
1828  sigPdfs = &negSignalPdfs_;
1829  scfPdfs = &negScfPdfs_;
1830  }
1831  if ( this->enableEmbedding() ) {
1832  embeddedData = posSignalTree_;
1833  doReweighting = usePosReweighting_;
1834  }
1835  }
1836 
1837  if ( this->useDP() ) {
1838  if ( embeddedData ) {
1839  if ( doReweighting ) {
1840  // Select a (random) event from the generated data. Then store the
1841  // reconstructed DP co-ords, together with other pdf information,
1842  // as the embedded data.
1843  genOK = embeddedData->getReweightedEvent( model );
1844  } else {
1845  // Just get the information of a (randomly) selected event in the
1846  // embedded data
1847  embeddedData->getEmbeddedEvent( kinematics );
1848  }
1849  genSCF = this->storeSignalMCMatch( embeddedData );
1850  } else {
1851  genOK = model->generate();
1852  if ( genOK && useSCF_ ) {
1853  Double_t frac( 0.0 );
1854  if ( useSCFHist_ ) {
1855  frac = scfFracHist_->calcEfficiency( kinematics );
1856  } else {
1857  frac = scfFrac_.genValue();
1858  }
1859  if ( frac < LauRandom::randomFun()->Rndm() ) {
1860  this->setGenNtupleIntegerBranchValue( "genTMSig", 1 );
1861  this->setGenNtupleIntegerBranchValue( "genSCFSig", 0 );
1862  genSCF = kFALSE;
1863  } else {
1864  this->setGenNtupleIntegerBranchValue( "genTMSig", 0 );
1865  this->setGenNtupleIntegerBranchValue( "genSCFSig", 1 );
1866  genSCF = kTRUE;
1867 
1868  // Optionally smear the DP position
1869  // of the SCF event
1870  if ( scfMap_ != 0 ) {
1871  Double_t xCoord( 0.0 ), yCoord( 0.0 );
1872  if ( kinematics->squareDP() ) {
1873  xCoord = kinematics->getmPrime();
1874  yCoord = kinematics->getThetaPrime();
1875  } else {
1876  xCoord = kinematics->getm13Sq();
1877  yCoord = kinematics->getm23Sq();
1878  }
1879 
1880  // Find the bin number where this event is generated
1881  Int_t binNo = scfMap_->binNumber( xCoord, yCoord );
1882 
1883  // Retrieve the migration histogram
1884  TH2* histo = scfMap_->trueHist( binNo );
1885 
1886  const LauAbsEffModel* effModel = model->getEffModel();
1887  do {
1888  // Get a random point from the histogram
1889  histo->GetRandom2( xCoord, yCoord );
1890 
1891  // Update the kinematics
1892  if ( kinematics->squareDP() ) {
1893  kinematics->updateSqDPKinematics( xCoord, yCoord );
1894  } else {
1895  kinematics->updateKinematics( xCoord, yCoord );
1896  }
1897  } while ( ! effModel->passVeto( kinematics ) );
1898  }
1899  }
1900  }
1901  }
1902  } else {
1903  if ( embeddedData ) {
1904  embeddedData->getEmbeddedEvent( 0 );
1905  genSCF = this->storeSignalMCMatch( embeddedData );
1906  } else if ( useSCF_ ) {
1907  Double_t frac = scfFrac_.genValue();
1908  if ( frac < LauRandom::randomFun()->Rndm() ) {
1909  this->setGenNtupleIntegerBranchValue( "genTMSig", 1 );
1910  this->setGenNtupleIntegerBranchValue( "genSCFSig", 0 );
1911  genSCF = kFALSE;
1912  } else {
1913  this->setGenNtupleIntegerBranchValue( "genTMSig", 0 );
1914  this->setGenNtupleIntegerBranchValue( "genSCFSig", 1 );
1915  genSCF = kTRUE;
1916  }
1917  }
1918  }
1919  if ( genOK ) {
1920  if ( useSCF_ ) {
1921  if ( genSCF ) {
1922  this->generateExtraPdfValues( scfPdfs, embeddedData );
1923  } else {
1924  this->generateExtraPdfValues( sigPdfs, embeddedData );
1925  }
1926  } else {
1927  this->generateExtraPdfValues( sigPdfs, embeddedData );
1928  }
1929  }
1930  // Check for problems with the embedding
1931  if ( embeddedData && ( embeddedData->nEvents() == embeddedData->nUsedEvents() ) ) {
1932  std::cerr << "WARNING in LauCPFitModel::generateSignalEvent : Source of embedded signal events used up, clearing the list of used events."
1933  << std::endl;
1934  embeddedData->clearUsedList();
1935  }
1936 
1937  return genOK;
1938 }
1939 
1940 Bool_t LauCPFitModel::generateBkgndEvent( UInt_t bkgndID )
1941 {
1942  // Generate Bkgnd event
1943  Bool_t genOK( kTRUE );
1944 
1945  LauAbsBkgndDPModel* model( 0 );
1946  LauEmbeddedData* embeddedData( 0 );
1947  LauPdfPList* extraPdfs( 0 );
1948  LauKinematics* kinematics( 0 );
1949 
1950  if ( curEvtCharge_ < 0 ) {
1951  model = negBkgndDPModels_[bkgndID];
1952  if ( this->enableEmbedding() ) {
1953  embeddedData = negBkgndTree_[bkgndID];
1954  }
1955  extraPdfs = &negBkgndPdfs_[bkgndID];
1956  kinematics = negKinematics_;
1957  } else {
1958  model = posBkgndDPModels_[bkgndID];
1959  if ( this->enableEmbedding() ) {
1960  embeddedData = posBkgndTree_[bkgndID];
1961  }
1962  if ( tagged_ ) {
1963  extraPdfs = &posBkgndPdfs_[bkgndID];
1964  } else {
1965  extraPdfs = &negBkgndPdfs_[bkgndID];
1966  }
1967  kinematics = posKinematics_;
1968  }
1969 
1970  if ( this->useDP() ) {
1971  if ( embeddedData ) {
1972  embeddedData->getEmbeddedEvent( kinematics );
1973  } else {
1974  if ( model == 0 ) {
1975  const TString& bkgndClass = this->bkgndClassName( bkgndID );
1976  std::cerr << "ERROR in LauCPFitModel::generateBkgndEvent : Can't find the DP model for background class \""
1977  << bkgndClass << "\"." << std::endl;
1978  gSystem->Exit( EXIT_FAILURE );
1979  }
1980  genOK = model->generate();
1981  }
1982  } else {
1983  if ( embeddedData ) {
1984  embeddedData->getEmbeddedEvent( 0 );
1985  }
1986  }
1987  if ( genOK ) {
1988  this->generateExtraPdfValues( extraPdfs, embeddedData );
1989  }
1990 
1991  // Check for problems with the embedding
1992  if ( embeddedData && ( embeddedData->nEvents() == embeddedData->nUsedEvents() ) ) {
1993  const TString& bkgndClass = this->bkgndClassName( bkgndID );
1994  std::cerr << "WARNING in LauCPFitModel::generateBkgndEvent : Source of embedded "
1995  << bkgndClass << " events used up, clearing the list of used events." << std::endl;
1996  embeddedData->clearUsedList();
1997  }
1998 
1999  return genOK;
2000 }
2001 
2002 void LauCPFitModel::setupGenNtupleBranches()
2003 {
2004  // Setup the required ntuple branches
2005  this->addGenNtupleDoubleBranch( "evtWeight" );
2006  this->addGenNtupleIntegerBranch( "genSig" );
2007  this->addGenNtupleDoubleBranch( "efficiency" );
2008  if ( useSCF_ || ( this->enableEmbedding() && negSignalTree_ != 0 &&
2009  negSignalTree_->haveBranch( "mcMatch" ) && posSignalTree_ != 0 &&
2010  posSignalTree_->haveBranch( "mcMatch" ) ) ) {
2011  this->addGenNtupleIntegerBranch( "genTMSig" );
2012  this->addGenNtupleIntegerBranch( "genSCFSig" );
2013  }
2014  const UInt_t nBkgnds = this->nBkgndClasses();
2015  for ( UInt_t iBkgnd( 0 ); iBkgnd < nBkgnds; ++iBkgnd ) {
2016  TString name( this->bkgndClassName( iBkgnd ) );
2017  name.Prepend( "gen" );
2018  this->addGenNtupleIntegerBranch( name );
2019  }
2020  this->addGenNtupleIntegerBranch( "charge" );
2021  if ( this->useDP() == kTRUE ) {
2022  this->addGenNtupleDoubleBranch( "m12" );
2023  this->addGenNtupleDoubleBranch( "m23" );
2024  this->addGenNtupleDoubleBranch( "m13" );
2025  this->addGenNtupleDoubleBranch( "m12Sq" );
2026  this->addGenNtupleDoubleBranch( "m23Sq" );
2027  this->addGenNtupleDoubleBranch( "m13Sq" );
2028  this->addGenNtupleDoubleBranch( "cosHel12" );
2029  this->addGenNtupleDoubleBranch( "cosHel23" );
2030  this->addGenNtupleDoubleBranch( "cosHel13" );
2031  if ( negKinematics_->squareDP() && posKinematics_->squareDP() ) {
2032  this->addGenNtupleDoubleBranch( "mPrime" );
2033  this->addGenNtupleDoubleBranch( "thPrime" );
2034  }
2035  }
2036  for ( LauPdfPList::const_iterator pdf_iter = negSignalPdfs_.begin();
2037  pdf_iter != negSignalPdfs_.end();
2038  ++pdf_iter ) {
2039  std::vector<TString> varNames = ( *pdf_iter )->varNames();
2040  for ( std::vector<TString>::const_iterator var_iter = varNames.begin();
2041  var_iter != varNames.end();
2042  ++var_iter ) {
2043  if ( ( *var_iter ) != "m13Sq" && ( *var_iter ) != "m23Sq" ) {
2044  this->addGenNtupleDoubleBranch( ( *var_iter ) );
2045  }
2046  }
2047  }
2048 }
2049 
2050 void LauCPFitModel::setDPBranchValues()
2051 {
2052  LauKinematics* kinematics( 0 );
2053  if ( curEvtCharge_ < 0 ) {
2054  kinematics = negKinematics_;
2055  } else {
2056  kinematics = posKinematics_;
2057  }
2058 
2059  // Store all the DP information
2060  this->setGenNtupleDoubleBranchValue( "m12", kinematics->getm12() );
2061  this->setGenNtupleDoubleBranchValue( "m23", kinematics->getm23() );
2062  this->setGenNtupleDoubleBranchValue( "m13", kinematics->getm13() );
2063  this->setGenNtupleDoubleBranchValue( "m12Sq", kinematics->getm12Sq() );
2064  this->setGenNtupleDoubleBranchValue( "m23Sq", kinematics->getm23Sq() );
2065  this->setGenNtupleDoubleBranchValue( "m13Sq", kinematics->getm13Sq() );
2066  this->setGenNtupleDoubleBranchValue( "cosHel12", kinematics->getc12() );
2067  this->setGenNtupleDoubleBranchValue( "cosHel23", kinematics->getc23() );
2068  this->setGenNtupleDoubleBranchValue( "cosHel13", kinematics->getc13() );
2069  if ( kinematics->squareDP() ) {
2070  this->setGenNtupleDoubleBranchValue( "mPrime", kinematics->getmPrime() );
2071  this->setGenNtupleDoubleBranchValue( "thPrime", kinematics->getThetaPrime() );
2072  }
2073 }
2074 
2075 void LauCPFitModel::generateExtraPdfValues( LauPdfPList* extraPdfs, LauEmbeddedData* embeddedData )
2076 {
2077  LauKinematics* kinematics( 0 );
2078  if ( curEvtCharge_ < 0 ) {
2079  kinematics = negKinematics_;
2080  } else {
2081  kinematics = posKinematics_;
2082  }
2083 
2084  if ( ! extraPdfs ) {
2085  std::cerr << "ERROR in LauCPFitModel::generateExtraPdfValues : Null pointer to PDF list."
2086  << std::endl;
2087  gSystem->Exit( EXIT_FAILURE );
2088  }
2089 
2090  if ( extraPdfs->empty() ) {
2091  //std::cerr << "WARNING in LauCPFitModel::generateExtraPdfValues : PDF list is empty." << std::endl;
2092  return;
2093  }
2094 
2095  // Generate from the extra PDFs
2096  for ( LauPdfPList::iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end();
2097  ++pdf_iter ) {
2098  LauFitData genValues;
2099  if ( embeddedData ) {
2100  genValues = embeddedData->getValues( ( *pdf_iter )->varNames() );
2101  } else {
2102  genValues = ( *pdf_iter )->generate( kinematics );
2103  }
2104  for ( LauFitData::const_iterator var_iter = genValues.begin(); var_iter != genValues.end();
2105  ++var_iter ) {
2106  TString varName = var_iter->first;
2107  if ( varName != "m13Sq" && varName != "m23Sq" ) {
2108  Double_t value = var_iter->second;
2109  this->setGenNtupleDoubleBranchValue( varName, value );
2110  }
2111  }
2112  }
2113 }
2114 
2115 Bool_t LauCPFitModel::storeSignalMCMatch( LauEmbeddedData* embeddedData )
2116 {
2117  // Default to TM
2118  Bool_t genSCF( kFALSE );
2119  Int_t match( 1 );
2120 
2121  // Check that we have a valid pointer and that embedded data has
2122  // the mcMatch branch. If so then get the match value.
2123  if ( embeddedData && embeddedData->haveBranch( "mcMatch" ) ) {
2124  match = TMath::Nint( embeddedData->getValue( "mcMatch" ) );
2125  }
2126 
2127  // Set the variables accordingly.
2128  if ( match ) {
2129  this->setGenNtupleIntegerBranchValue( "genTMSig", 1 );
2130  this->setGenNtupleIntegerBranchValue( "genSCFSig", 0 );
2131  genSCF = kFALSE;
2132  } else {
2133  this->setGenNtupleIntegerBranchValue( "genTMSig", 0 );
2134  this->setGenNtupleIntegerBranchValue( "genSCFSig", 1 );
2135  genSCF = kTRUE;
2136  }
2137 
2138  return genSCF;
2139 }
2140 
2141 void LauCPFitModel::propagateParUpdates()
2142 {
2143  // Update the signal parameters and then the total normalisation for the signal likelihood
2144  if ( this->useDP() == kTRUE ) {
2145  this->updateCoeffs();
2146  negSigModel_->updateCoeffs( negCoeffs_ );
2147  posSigModel_->updateCoeffs( posCoeffs_ );
2148  }
2149 
2150  // Update the signal fraction from the background fractions if not doing an extended fit
2151  if ( ! this->doEMLFit() && ! signalEvents_->fixed() ) {
2152  this->updateSigEvents();
2153  }
2154 }
2155 
2156 void LauCPFitModel::updateSigEvents()
2157 {
2158  // The background parameters will have been set from Minuit.
2159  // We need to update the signal events using these.
2160  Double_t nTotEvts = this->eventsPerExpt();
2161 
2162  signalEvents_->range( -2.0 * nTotEvts, 2.0 * nTotEvts );
2163  for ( LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end();
2164  ++iter ) {
2165  LauAbsRValue* nBkgndEvents = ( *iter );
2166  if ( nBkgndEvents->isLValue() ) {
2167  LauParameter* yield = dynamic_cast<LauParameter*>( nBkgndEvents );
2168  yield->range( -2.0 * nTotEvts, 2.0 * nTotEvts );
2169  }
2170  }
2171 
2172  if ( signalEvents_->fixed() ) {
2173  return;
2174  }
2175 
2176  // Subtract background events (if any) from signal.
2177  Double_t signalEvents = nTotEvts;
2178  if ( usingBkgnd_ == kTRUE ) {
2179  for ( LauBkgndYieldList::const_iterator iter = bkgndEvents_.begin();
2180  iter != bkgndEvents_.end();
2181  ++iter ) {
2182  signalEvents -= ( *iter )->value();
2183  }
2184  }
2185  signalEvents_->value( signalEvents );
2186 }
2187 
2188 void LauCPFitModel::cacheInputFitVars()
2189 {
2190  // Fill the internal data trees of the signal and background models.
2191  // Note that we store the events of both charges in both the
2192  // negative and the positive models. It's only later, at the stage
2193  // when the likelihood is being calculated, that we separate them.
2194 
2195  LauFitDataTree* inputFitData = this->fitData();
2196 
2197  // First the Dalitz plot variables (m_ij^2)
2198  if ( this->useDP() == kTRUE ) {
2199 
2200  // need to append SCF smearing bins before caching DP amplitudes
2201  if ( scfMap_ != 0 ) {
2202  this->appendBinCentres( inputFitData );
2203  }
2204 
2205  negSigModel_->fillDataTree( *inputFitData );
2206  posSigModel_->fillDataTree( *inputFitData );
2207 
2208  if ( usingBkgnd_ == kTRUE ) {
2209  for ( LauBkgndDPModelList::iterator iter = negBkgndDPModels_.begin();
2210  iter != negBkgndDPModels_.end();
2211  ++iter ) {
2212  ( *iter )->fillDataTree( *inputFitData );
2213  }
2214  for ( LauBkgndDPModelList::iterator iter = posBkgndDPModels_.begin();
2215  iter != posBkgndDPModels_.end();
2216  ++iter ) {
2217  ( *iter )->fillDataTree( *inputFitData );
2218  }
2219  }
2220  }
2221 
2222  // ...and then the extra PDFs
2223  this->cacheInfo( negSignalPdfs_, *inputFitData );
2224  this->cacheInfo( negScfPdfs_, *inputFitData );
2225  for ( LauBkgndPdfsList::iterator iter = negBkgndPdfs_.begin(); iter != negBkgndPdfs_.end();
2226  ++iter ) {
2227  this->cacheInfo( ( *iter ), *inputFitData );
2228  }
2229  if ( tagged_ ) {
2230  this->cacheInfo( posSignalPdfs_, *inputFitData );
2231  this->cacheInfo( posScfPdfs_, *inputFitData );
2232  for ( LauBkgndPdfsList::iterator iter = posBkgndPdfs_.begin(); iter != posBkgndPdfs_.end();
2233  ++iter ) {
2234  this->cacheInfo( ( *iter ), *inputFitData );
2235  }
2236  }
2237 
2238  // the SCF fractions and jacobians
2239  if ( useSCF_ && useSCFHist_ ) {
2240  if ( ! inputFitData->haveBranch( "m13Sq" ) || ! inputFitData->haveBranch( "m23Sq" ) ) {
2241  std::cerr << "ERROR in LauCPFitModel::cacheInputFitVars : Input data does not contain DP branches and so can't cache the SCF fraction."
2242  << std::endl;
2243  gSystem->Exit( EXIT_FAILURE );
2244  }
2245 
2246  UInt_t nEvents = inputFitData->nEvents();
2247  recoSCFFracs_.clear();
2248  recoSCFFracs_.reserve( nEvents );
2249  if ( negKinematics_->squareDP() ) {
2250  recoJacobians_.clear();
2251  recoJacobians_.reserve( nEvents );
2252  }
2253  for ( UInt_t iEvt = 0; iEvt < nEvents; iEvt++ ) {
2254  const LauFitData& dataValues = inputFitData->getData( iEvt );
2255  LauFitData::const_iterator m13_iter = dataValues.find( "m13Sq" );
2256  LauFitData::const_iterator m23_iter = dataValues.find( "m23Sq" );
2257  negKinematics_->updateKinematics( m13_iter->second, m23_iter->second );
2258  Double_t scfFrac = scfFracHist_->calcEfficiency( negKinematics_ );
2259  recoSCFFracs_.push_back( scfFrac );
2260  if ( negKinematics_->squareDP() ) {
2261  recoJacobians_.push_back( negKinematics_->calcSqDPJacobian() );
2262  }
2263  }
2264  }
2265 
2266  // finally cache the event charge
2267  evtCharges_.clear();
2268  if ( tagged_ ) {
2269  if ( ! inputFitData->haveBranch( tagVarName_ ) ) {
2270  std::cerr << "ERROR in LauCPFitModel::cacheInputFitVars : Input data does not contain branch \""
2271  << tagVarName_ << "\"." << std::endl;
2272  gSystem->Exit( EXIT_FAILURE );
2273  }
2274  UInt_t nEvents = inputFitData->nEvents();
2275  evtCharges_.reserve( nEvents );
2276  for ( UInt_t iEvt = 0; iEvt < nEvents; iEvt++ ) {
2277  const LauFitData& dataValues = inputFitData->getData( iEvt );
2278  LauFitData::const_iterator iter = dataValues.find( tagVarName_ );
2279  curEvtCharge_ = static_cast<Int_t>( iter->second );
2280  evtCharges_.push_back( curEvtCharge_ );
2281  }
2282  }
2283 }
2284 
2285 void LauCPFitModel::appendBinCentres( LauFitDataTree* inputData )
2286 {
2287  // We'll be caching the DP amplitudes and efficiencies of the centres of the true bins.
2288  // To do so, we attach some fake points at the end of inputData, the number of the entry
2289  // minus the total number of events corresponding to the number of the histogram for that
2290  // given true bin in the LauScfMap object. (What this means is that when Laura is provided with
2291  // the LauScfMap object by the user, it's the latter who has to make sure that it contains the
2292  // right number of histograms and in exactly the right order!)
2293 
2294  // Get the x and y co-ordinates of the bin centres
2295  std::vector<Double_t> binCentresXCoords;
2296  std::vector<Double_t> binCentresYCoords;
2297  scfMap_->listBinCentres( binCentresXCoords, binCentresYCoords );
2298 
2299  // The SCF histograms could be in square Dalitz plot histograms.
2300  // The dynamics takes normal Dalitz plot coords, so we might have to convert them back.
2301  Bool_t sqDP = negKinematics_->squareDP();
2302  UInt_t nBins = binCentresXCoords.size();
2303  fakeSCFFracs_.clear();
2304  fakeSCFFracs_.reserve( nBins );
2305  if ( sqDP ) {
2306  fakeJacobians_.clear();
2307  fakeJacobians_.reserve( nBins );
2308  }
2309 
2310  for ( UInt_t iBin = 0; iBin < nBins; ++iBin ) {
2311 
2312  if ( sqDP ) {
2313  negKinematics_->updateSqDPKinematics( binCentresXCoords[iBin], binCentresYCoords[iBin] );
2314  binCentresXCoords[iBin] = negKinematics_->getm13Sq();
2315  binCentresYCoords[iBin] = negKinematics_->getm23Sq();
2316  fakeJacobians_.push_back( negKinematics_->calcSqDPJacobian() );
2317  } else {
2318  negKinematics_->updateKinematics( binCentresXCoords[iBin], binCentresYCoords[iBin] );
2319  }
2320 
2321  fakeSCFFracs_.push_back( scfFracHist_->calcEfficiency( negKinematics_ ) );
2322  }
2323 
2324  // Set up inputFitVars_ object to hold the fake events
2325  inputData->appendFakePoints( binCentresXCoords, binCentresYCoords );
2326 }
2327 
2328 Double_t LauCPFitModel::getTotEvtLikelihood( UInt_t iEvt )
2329 {
2330  // Find out whether we have B- or B+
2331  if ( tagged_ ) {
2332  curEvtCharge_ = evtCharges_[iEvt];
2333 
2334  // check that the charge is either +1 or -1
2335  if ( TMath::Abs( curEvtCharge_ ) != 1 ) {
2336  std::cerr << "ERROR in LauCPFitModel::getTotEvtLikelihood : Charge/tag not accepted value: "
2337  << curEvtCharge_ << std::endl;
2338  if ( curEvtCharge_ > 0 ) {
2339  curEvtCharge_ = +1;
2340  } else {
2341  curEvtCharge_ = -1;
2342  }
2343  std::cerr << " : Making it: " << curEvtCharge_
2344  << "." << std::endl;
2345  }
2346  }
2347 
2348  // Get the DP likelihood for signal and backgrounds
2349  this->getEvtDPLikelihood( iEvt );
2350 
2351  // Get the combined extra PDFs likelihood for signal and backgrounds
2352  this->getEvtExtraLikelihoods( iEvt );
2353 
2354  // If appropriate, combine the TM and SCF likelihoods
2355  Double_t sigLike = sigDPLike_ * sigExtraLike_;
2356  if ( useSCF_ ) {
2357  Double_t scfFrac( 0.0 );
2358  if ( useSCFHist_ ) {
2359  scfFrac = recoSCFFracs_[iEvt];
2360  } else {
2361  scfFrac = scfFrac_.unblindValue();
2362  }
2363  sigLike *= ( 1.0 - scfFrac );
2364  if ( ( scfMap_ != 0 ) && ( this->useDP() == kTRUE ) ) {
2365  // if we're smearing the SCF DP PDF then the SCF frac
2366  // is already included in the SCF DP likelihood
2367  sigLike += ( scfDPLike_ * scfExtraLike_ );
2368  } else {
2369  sigLike += ( scfFrac * scfDPLike_ * scfExtraLike_ );
2370  }
2371  }
2372 
2373  // Get the correct event fractions depending on the charge
2374  // Signal asymmetry is built into the DP model... but when the DP
2375  // isn't in the fit we need an explicit parameter
2376  Double_t signalEvents = signalEvents_->unblindValue() * 0.5;
2377  if ( this->useDP() == kFALSE ) {
2378  signalEvents *= ( 1.0 - curEvtCharge_ * signalAsym_->unblindValue() );
2379  }
2380 
2381  // Construct the total event likelihood
2382  Double_t likelihood( 0.0 );
2383  if ( usingBkgnd_ ) {
2384  likelihood = sigLike * signalEvents;
2385  const UInt_t nBkgnds = this->nBkgndClasses();
2386  for ( UInt_t bkgndID( 0 ); bkgndID < nBkgnds; ++bkgndID ) {
2387  Double_t bkgndEvents = bkgndEvents_[bkgndID]->unblindValue() * 0.5 *
2388  ( 1.0 - curEvtCharge_ * bkgndAsym_[bkgndID]->unblindValue() );
2389  likelihood += bkgndEvents * bkgndDPLike_[bkgndID] * bkgndExtraLike_[bkgndID];
2390  }
2391  } else {
2392  likelihood = sigLike * 0.5;
2393  }
2394  return likelihood;
2395 }
2396 
2397 Double_t LauCPFitModel::getEventSum() const
2398 {
2399  Double_t eventSum( 0.0 );
2400  eventSum += signalEvents_->unblindValue();
2401  if ( usingBkgnd_ ) {
2402  for ( LauBkgndYieldList::const_iterator iter = bkgndEvents_.begin();
2403  iter != bkgndEvents_.end();
2404  ++iter ) {
2405  eventSum += ( *iter )->unblindValue();
2406  }
2407  }
2408  return eventSum;
2409 }
2410 
2411 void LauCPFitModel::getEvtDPLikelihood( UInt_t iEvt )
2412 {
2413  // Function to return the signal and background likelihoods for the
2414  // Dalitz plot for the given event evtNo.
2415 
2416  if ( ! this->useDP() ) {
2417  // There's always going to be a term in the likelihood for the
2418  // signal, so we'd better not zero it.
2419  sigDPLike_ = 1.0;
2420  scfDPLike_ = 1.0;
2421 
2422  const UInt_t nBkgnds = this->nBkgndClasses();
2423  for ( UInt_t bkgndID( 0 ); bkgndID < nBkgnds; ++bkgndID ) {
2424  if ( usingBkgnd_ == kTRUE ) {
2425  bkgndDPLike_[bkgndID] = 1.0;
2426  } else {
2427  bkgndDPLike_[bkgndID] = 0.0;
2428  }
2429  }
2430 
2431  return;
2432  }
2433 
2434  const UInt_t nBkgnds = this->nBkgndClasses();
2435  if ( tagged_ ) {
2436  if ( curEvtCharge_ == +1 ) {
2437  posSigModel_->calcLikelihoodInfo( iEvt );
2438  sigDPLike_ = posSigModel_->getEvtIntensity();
2439 
2440  for ( UInt_t bkgndID( 0 ); bkgndID < nBkgnds; ++bkgndID ) {
2441  if ( usingBkgnd_ == kTRUE ) {
2442  bkgndDPLike_[bkgndID] = posBkgndDPModels_[bkgndID]->getLikelihood( iEvt );
2443  } else {
2444  bkgndDPLike_[bkgndID] = 0.0;
2445  }
2446  }
2447  } else {
2448  negSigModel_->calcLikelihoodInfo( iEvt );
2449  sigDPLike_ = negSigModel_->getEvtIntensity();
2450 
2451  for ( UInt_t bkgndID( 0 ); bkgndID < nBkgnds; ++bkgndID ) {
2452  if ( usingBkgnd_ == kTRUE ) {
2453  bkgndDPLike_[bkgndID] = negBkgndDPModels_[bkgndID]->getLikelihood( iEvt );
2454  } else {
2455  bkgndDPLike_[bkgndID] = 0.0;
2456  }
2457  }
2458  }
2459  } else {
2460  posSigModel_->calcLikelihoodInfo( iEvt );
2461  negSigModel_->calcLikelihoodInfo( iEvt );
2462 
2463  sigDPLike_ = 0.5 * ( posSigModel_->getEvtIntensity() + negSigModel_->getEvtIntensity() );
2464 
2465  for ( UInt_t bkgndID( 0 ); bkgndID < nBkgnds; ++bkgndID ) {
2466  if ( usingBkgnd_ == kTRUE ) {
2467  bkgndDPLike_[bkgndID] = 0.5 * ( posBkgndDPModels_[bkgndID]->getLikelihood( iEvt ) +
2468  negBkgndDPModels_[bkgndID]->getLikelihood( iEvt ) );
2469  } else {
2470  bkgndDPLike_[bkgndID] = 0.0;
2471  }
2472  }
2473  }
2474 
2475  if ( useSCF_ == kTRUE ) {
2476  if ( scfMap_ == 0 ) {
2477  // we're not smearing the SCF DP position
2478  // so the likelihood is the same as the TM
2479  scfDPLike_ = sigDPLike_;
2480  } else {
2481  // calculate the smeared SCF DP likelihood
2482  scfDPLike_ = this->getEvtSCFDPLikelihood( iEvt );
2483  }
2484  }
2485 
2486  // Calculate the signal normalisation
2487  // NB the 2.0 is there so that the 0.5 factor is applied to
2488  // signal and background in the same place otherwise you get
2489  // normalisation problems when you switch off the DP in the fit
2490  Double_t norm = negSigModel_->getDPNorm() + posSigModel_->getDPNorm();
2491  sigDPLike_ *= 2.0 / norm;
2492  scfDPLike_ *= 2.0 / norm;
2493 }
2494 
2495 Double_t LauCPFitModel::getEvtSCFDPLikelihood( UInt_t iEvt )
2496 {
2497  Double_t scfDPLike( 0.0 );
2498 
2499  Double_t recoJacobian( 1.0 );
2500  Double_t xCoord( 0.0 );
2501  Double_t yCoord( 0.0 );
2502 
2503  Bool_t squareDP = negKinematics_->squareDP();
2504  if ( squareDP ) {
2505  xCoord = negSigModel_->getEvtmPrime();
2506  yCoord = negSigModel_->getEvtthPrime();
2507  recoJacobian = recoJacobians_[iEvt];
2508  } else {
2509  xCoord = negSigModel_->getEvtm13Sq();
2510  yCoord = negSigModel_->getEvtm23Sq();
2511  }
2512 
2513  // Find the bin that our reco event falls in
2514  Int_t recoBin = scfMap_->binNumber( xCoord, yCoord );
2515 
2516  // Find out which true Bins contribute to the given reco bin
2517  const std::vector<Int_t>* trueBins = scfMap_->trueBins( recoBin );
2518 
2519  const Int_t nDataEvents = this->eventsPerExpt();
2520 
2521  // Loop over the true bins
2522  for ( std::vector<Int_t>::const_iterator iter = trueBins->begin(); iter != trueBins->end();
2523  ++iter ) {
2524  Int_t trueBin = ( *iter );
2525 
2526  // prob of a true event in the given true bin migrating to the reco bin
2527  Double_t pRecoGivenTrue = scfMap_->prob( recoBin, trueBin );
2528  Double_t pTrue( 0.0 );
2529 
2530  // We've cached the DP amplitudes and the efficiency for the
2531  // true bin centres, just after the data points
2532  if ( tagged_ ) {
2533  LauIsobarDynamics* sigModel( 0 );
2534  if ( curEvtCharge_ < 0 ) {
2535  sigModel = negSigModel_;
2536  } else {
2537  sigModel = posSigModel_;
2538  }
2539 
2540  sigModel->calcLikelihoodInfo( nDataEvents + trueBin );
2541 
2542  pTrue = sigModel->getEvtIntensity();
2543  } else {
2544  posSigModel_->calcLikelihoodInfo( nDataEvents + trueBin );
2545  negSigModel_->calcLikelihoodInfo( nDataEvents + trueBin );
2546 
2547  pTrue = 0.5 * ( posSigModel_->getEvtIntensity() + negSigModel_->getEvtIntensity() );
2548  }
2549 
2550  // Get the cached SCF fraction (and jacobian if we're using the square DP)
2551  Double_t scfFraction = fakeSCFFracs_[trueBin];
2552  Double_t jacobian( 1.0 );
2553  if ( squareDP ) {
2554  jacobian = fakeJacobians_[trueBin];
2555  }
2556 
2557  scfDPLike += pTrue * jacobian * scfFraction * pRecoGivenTrue;
2558  }
2559 
2560  // Divide by the reco jacobian
2561  scfDPLike /= recoJacobian;
2562 
2563  return scfDPLike;
2564 }
2565 
2566 void LauCPFitModel::getEvtExtraLikelihoods( UInt_t iEvt )
2567 {
2568  // Function to return the signal and background likelihoods for the
2569  // extra variables for the given event evtNo.
2570 
2571  sigExtraLike_ = 1.0;
2572 
2573  const UInt_t nBkgnds = this->nBkgndClasses();
2574 
2575  if ( ! tagged_ || curEvtCharge_ < 0 ) {
2576  sigExtraLike_ = this->prodPdfValue( negSignalPdfs_, iEvt );
2577 
2578  if ( useSCF_ ) {
2579  scfExtraLike_ = this->prodPdfValue( negScfPdfs_, iEvt );
2580  }
2581 
2582  for ( UInt_t bkgndID( 0 ); bkgndID < nBkgnds; ++bkgndID ) {
2583  if ( usingBkgnd_ ) {
2584  bkgndExtraLike_[bkgndID] = this->prodPdfValue( negBkgndPdfs_[bkgndID], iEvt );
2585  } else {
2586  bkgndExtraLike_[bkgndID] = 0.0;
2587  }
2588  }
2589  } else {
2590  sigExtraLike_ = this->prodPdfValue( posSignalPdfs_, iEvt );
2591 
2592  if ( useSCF_ ) {
2593  scfExtraLike_ = this->prodPdfValue( posScfPdfs_, iEvt );
2594  }
2595 
2596  for ( UInt_t bkgndID( 0 ); bkgndID < nBkgnds; ++bkgndID ) {
2597  if ( usingBkgnd_ ) {
2598  bkgndExtraLike_[bkgndID] = this->prodPdfValue( posBkgndPdfs_[bkgndID], iEvt );
2599  } else {
2600  bkgndExtraLike_[bkgndID] = 0.0;
2601  }
2602  }
2603  }
2604 }
2605 
2606 void LauCPFitModel::updateCoeffs()
2607 {
2608  negCoeffs_.clear();
2609  posCoeffs_.clear();
2610  negCoeffs_.reserve( nSigComp_ );
2611  posCoeffs_.reserve( nSigComp_ );
2612  for ( UInt_t i = 0; i < nSigComp_; i++ ) {
2613  negCoeffs_.push_back( coeffPars_[i]->antiparticleCoeff() );
2614  posCoeffs_.push_back( coeffPars_[i]->particleCoeff() );
2615  }
2616 }
2617 
2618 void LauCPFitModel::setupSPlotNtupleBranches()
2619 {
2620  // add branches for storing the experiment number and the number of
2621  // the event within the current experiment
2622  this->addSPlotNtupleIntegerBranch( "iExpt" );
2623  this->addSPlotNtupleIntegerBranch( "iEvtWithinExpt" );
2624 
2625  // Store the efficiency of the event (for inclusive BF calculations).
2626  if ( this->storeDPEff() ) {
2627  this->addSPlotNtupleDoubleBranch( "efficiency" );
2628  if ( negSigModel_->usingScfModel() && posSigModel_->usingScfModel() ) {
2629  this->addSPlotNtupleDoubleBranch( "scffraction" );
2630  }
2631  }
2632 
2633  // Store the total event likelihood for each species.
2634  if ( useSCF_ ) {
2635  this->addSPlotNtupleDoubleBranch( "sigTMTotalLike" );
2636  this->addSPlotNtupleDoubleBranch( "sigSCFTotalLike" );
2637  this->addSPlotNtupleDoubleBranch( "sigSCFFrac" );
2638  } else {
2639  this->addSPlotNtupleDoubleBranch( "sigTotalLike" );
2640  }
2641  if ( usingBkgnd_ ) {
2642  const UInt_t nBkgnds = this->nBkgndClasses();
2643  for ( UInt_t iBkgnd( 0 ); iBkgnd < nBkgnds; ++iBkgnd ) {
2644  TString name( this->bkgndClassName( iBkgnd ) );
2645  name += "TotalLike";
2646  this->addSPlotNtupleDoubleBranch( name );
2647  }
2648  }
2649 
2650  // Store the DP likelihoods
2651  if ( this->useDP() ) {
2652  if ( useSCF_ ) {
2653  this->addSPlotNtupleDoubleBranch( "sigTMDPLike" );
2654  this->addSPlotNtupleDoubleBranch( "sigSCFDPLike" );
2655  } else {
2656  this->addSPlotNtupleDoubleBranch( "sigDPLike" );
2657  }
2658  if ( usingBkgnd_ ) {
2659  const UInt_t nBkgnds = this->nBkgndClasses();
2660  for ( UInt_t iBkgnd( 0 ); iBkgnd < nBkgnds; ++iBkgnd ) {
2661  TString name( this->bkgndClassName( iBkgnd ) );
2662  name += "DPLike";
2663  this->addSPlotNtupleDoubleBranch( name );
2664  }
2665  }
2666  }
2667 
2668  // Store the likelihoods for each extra PDF
2669  if ( useSCF_ ) {
2670  this->addSPlotNtupleBranches( &negSignalPdfs_, "sigTM" );
2671  this->addSPlotNtupleBranches( &negScfPdfs_, "sigSCF" );
2672  } else {
2673  this->addSPlotNtupleBranches( &negSignalPdfs_, "sig" );
2674  }
2675  if ( usingBkgnd_ ) {
2676  const UInt_t nBkgnds = this->nBkgndClasses();
2677  for ( UInt_t iBkgnd( 0 ); iBkgnd < nBkgnds; ++iBkgnd ) {
2678  const TString& bkgndClass = this->bkgndClassName( iBkgnd );
2679  const LauPdfPList* pdfList = &( negBkgndPdfs_[iBkgnd] );
2680  this->addSPlotNtupleBranches( pdfList, bkgndClass );
2681  }
2682  }
2683 }
2684 
2685 void LauCPFitModel::addSPlotNtupleBranches( const LauPdfPList* extraPdfs, const TString& prefix )
2686 {
2687  if ( extraPdfs ) {
2688  // Loop through each of the PDFs
2689  for ( LauPdfPList::const_iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end();
2690  ++pdf_iter ) {
2691 
2692  // Count the number of input variables that are not
2693  // DP variables (used in the case where there is DP
2694  // dependence for e.g. MVA)
2695  UInt_t nVars( 0 );
2696  std::vector<TString> varNames = ( *pdf_iter )->varNames();
2697  for ( std::vector<TString>::const_iterator var_iter = varNames.begin();
2698  var_iter != varNames.end();
2699  ++var_iter ) {
2700  if ( ( *var_iter ) != "m13Sq" && ( *var_iter ) != "m23Sq" ) {
2701  ++nVars;
2702  }
2703  }
2704 
2705  if ( nVars == 1 ) {
2706  // If the PDF only has one variable then
2707  // simply add one branch for that variable
2708  TString varName = ( *pdf_iter )->varName();
2709  TString name( prefix );
2710  name += varName;
2711  name += "Like";
2712  this->addSPlotNtupleDoubleBranch( name );
2713  } else if ( nVars == 2 ) {
2714  // If the PDF has two variables then we
2715  // need a branch for them both together and
2716  // branches for each
2717  TString allVars( "" );
2718  for ( std::vector<TString>::const_iterator var_iter = varNames.begin();
2719  var_iter != varNames.end();
2720  ++var_iter ) {
2721  allVars += ( *var_iter );
2722  TString name( prefix );
2723  name += ( *var_iter );
2724  name += "Like";
2725  this->addSPlotNtupleDoubleBranch( name );
2726  }
2727  TString name( prefix );
2728  name += allVars;
2729  name += "Like";
2730  this->addSPlotNtupleDoubleBranch( name );
2731  } else {
2732  std::cerr << "WARNING in LauCPFitModel::addSPlotNtupleBranches : Can't yet deal with 3D PDFs."
2733  << std::endl;
2734  }
2735  }
2736  }
2737 }
2738 
2739 Double_t LauCPFitModel::setSPlotNtupleBranchValues( LauPdfPList* extraPdfs,
2740  const TString& prefix,
2741  UInt_t iEvt )
2742 {
2743  // Store the PDF value for each variable in the list
2744  Double_t totalLike( 1.0 );
2745  Double_t extraLike( 0.0 );
2746  if ( extraPdfs ) {
2747  for ( LauPdfPList::iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end();
2748  ++pdf_iter ) {
2749 
2750  // calculate the likelihood for this event
2751  ( *pdf_iter )->calcLikelihoodInfo( iEvt );
2752  extraLike = ( *pdf_iter )->getLikelihood();
2753  totalLike *= extraLike;
2754 
2755  // Count the number of input variables that are not
2756  // DP variables (used in the case where there is DP
2757  // dependence for e.g. MVA)
2758  UInt_t nVars( 0 );
2759  std::vector<TString> varNames = ( *pdf_iter )->varNames();
2760  for ( std::vector<TString>::const_iterator var_iter = varNames.begin();
2761  var_iter != varNames.end();
2762  ++var_iter ) {
2763  if ( ( *var_iter ) != "m13Sq" && ( *var_iter ) != "m23Sq" ) {
2764  ++nVars;
2765  }
2766  }
2767 
2768  if ( nVars == 1 ) {
2769  // If the PDF only has one variable then
2770  // simply store the value for that variable
2771  TString varName = ( *pdf_iter )->varName();
2772  TString name( prefix );
2773  name += varName;
2774  name += "Like";
2775  this->setSPlotNtupleDoubleBranchValue( name, extraLike );
2776  } else if ( nVars == 2 ) {
2777  // If the PDF has two variables then we
2778  // store the value for them both together
2779  // and for each on their own
2780  TString allVars( "" );
2781  for ( std::vector<TString>::const_iterator var_iter = varNames.begin();
2782  var_iter != varNames.end();
2783  ++var_iter ) {
2784  allVars += ( *var_iter );
2785  TString name( prefix );
2786  name += ( *var_iter );
2787  name += "Like";
2788  Double_t indivLike = ( *pdf_iter )->getLikelihood( ( *var_iter ) );
2789  this->setSPlotNtupleDoubleBranchValue( name, indivLike );
2790  }
2791  TString name( prefix );
2792  name += allVars;
2793  name += "Like";
2794  this->setSPlotNtupleDoubleBranchValue( name, extraLike );
2795  } else {
2796  std::cerr << "WARNING in LauCPFitModel::setSPlotNtupleBranchValues : Can't yet deal with 3D PDFs."
2797  << std::endl;
2798  }
2799  }
2800  }
2801  return totalLike;
2802 }
2803 
2804 LauSPlot::NameSet LauCPFitModel::variableNames() const
2805 {
2806  LauSPlot::NameSet nameSet;
2807  if ( this->useDP() ) {
2808  nameSet.insert( "DP" );
2809  }
2810  // Loop through all the signal PDFs
2811  for ( LauPdfPList::const_iterator pdf_iter = negSignalPdfs_.begin();
2812  pdf_iter != negSignalPdfs_.end();
2813  ++pdf_iter ) {
2814  // Loop over the variables involved in each PDF
2815  std::vector<TString> varNames = ( *pdf_iter )->varNames();
2816  for ( std::vector<TString>::const_iterator var_iter = varNames.begin();
2817  var_iter != varNames.end();
2818  ++var_iter ) {
2819  // If they are not DP coordinates then add them
2820  if ( ( *var_iter ) != "m13Sq" && ( *var_iter ) != "m23Sq" ) {
2821  nameSet.insert( ( *var_iter ) );
2822  }
2823  }
2824  }
2825  return nameSet;
2826 }
2827 
2828 LauSPlot::NumbMap LauCPFitModel::freeSpeciesNames() const
2829 {
2830  LauSPlot::NumbMap numbMap;
2831  if ( ! signalEvents_->fixed() && this->doEMLFit() ) {
2832  numbMap["sig"] = signalEvents_->genValue();
2833  }
2834  if ( usingBkgnd_ ) {
2835  const UInt_t nBkgnds = this->nBkgndClasses();
2836  for ( UInt_t iBkgnd( 0 ); iBkgnd < nBkgnds; ++iBkgnd ) {
2837  const TString& bkgndClass = this->bkgndClassName( iBkgnd );
2838  const LauAbsRValue* par = bkgndEvents_[iBkgnd];
2839  if ( ! par->fixed() ) {
2840  numbMap[bkgndClass] = par->genValue();
2841  if ( ! par->isLValue() ) {
2842  std::cerr << "WARNING in LauCPFitModel::freeSpeciesNames : \"" << par->name()
2843  << "\" is a LauFormulaPar, which implies it is perhaps not entirely free to float in the fit, so the sWeight calculation may not be reliable"
2844  << std::endl;
2845  }
2846  }
2847  }
2848  }
2849  return numbMap;
2850 }
2851 
2852 LauSPlot::NumbMap LauCPFitModel::fixdSpeciesNames() const
2853 {
2854  LauSPlot::NumbMap numbMap;
2855  if ( signalEvents_->fixed() && this->doEMLFit() ) {
2856  numbMap["sig"] = signalEvents_->genValue();
2857  }
2858  if ( usingBkgnd_ ) {
2859  const UInt_t nBkgnds = this->nBkgndClasses();
2860  for ( UInt_t iBkgnd( 0 ); iBkgnd < nBkgnds; ++iBkgnd ) {
2861  const TString& bkgndClass = this->bkgndClassName( iBkgnd );
2862  const LauAbsRValue* par = bkgndEvents_[iBkgnd];
2863  if ( par->fixed() ) {
2864  numbMap[bkgndClass] = par->genValue();
2865  }
2866  }
2867  }
2868  return numbMap;
2869 }
2870 
2871 LauSPlot::TwoDMap LauCPFitModel::twodimPDFs() const
2872 {
2873  // This makes the assumption that the form of the positive and
2874  // negative PDFs are the same, which seems reasonable to me
2875 
2876  LauSPlot::TwoDMap twodimMap;
2877 
2878  for ( LauPdfPList::const_iterator pdf_iter = negSignalPdfs_.begin();
2879  pdf_iter != negSignalPdfs_.end();
2880  ++pdf_iter ) {
2881  // Count the number of input variables that are not DP variables
2882  UInt_t nVars( 0 );
2883  std::vector<TString> varNames = ( *pdf_iter )->varNames();
2884  for ( std::vector<TString>::const_iterator var_iter = varNames.begin();
2885  var_iter != varNames.end();
2886  ++var_iter ) {
2887  if ( ( *var_iter ) != "m13Sq" && ( *var_iter ) != "m23Sq" ) {
2888  ++nVars;
2889  }
2890  }
2891  if ( nVars == 2 ) {
2892  if ( useSCF_ ) {
2893  twodimMap.insert(
2894  std::make_pair( "sigTM", std::make_pair( varNames[0], varNames[1] ) ) );
2895  } else {
2896  twodimMap.insert(
2897  std::make_pair( "sig", std::make_pair( varNames[0], varNames[1] ) ) );
2898  }
2899  }
2900  }
2901 
2902  if ( useSCF_ ) {
2903  for ( LauPdfPList::const_iterator pdf_iter = negScfPdfs_.begin();
2904  pdf_iter != negScfPdfs_.end();
2905  ++pdf_iter ) {
2906  // Count the number of input variables that are not DP variables
2907  UInt_t nVars( 0 );
2908  std::vector<TString> varNames = ( *pdf_iter )->varNames();
2909  for ( std::vector<TString>::const_iterator var_iter = varNames.begin();
2910  var_iter != varNames.end();
2911  ++var_iter ) {
2912  if ( ( *var_iter ) != "m13Sq" && ( *var_iter ) != "m23Sq" ) {
2913  ++nVars;
2914  }
2915  }
2916  if ( nVars == 2 ) {
2917  twodimMap.insert(
2918  std::make_pair( "sigSCF", std::make_pair( varNames[0], varNames[1] ) ) );
2919  }
2920  }
2921  }
2922 
2923  if ( usingBkgnd_ ) {
2924  const UInt_t nBkgnds = this->nBkgndClasses();
2925  for ( UInt_t iBkgnd( 0 ); iBkgnd < nBkgnds; ++iBkgnd ) {
2926  const TString& bkgndClass = this->bkgndClassName( iBkgnd );
2927  const LauPdfPList& pdfList = negBkgndPdfs_[iBkgnd];
2928  for ( LauPdfPList::const_iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end();
2929  ++pdf_iter ) {
2930  // Count the number of input variables that are not DP variables
2931  UInt_t nVars( 0 );
2932  std::vector<TString> varNames = ( *pdf_iter )->varNames();
2933  for ( std::vector<TString>::const_iterator var_iter = varNames.begin();
2934  var_iter != varNames.end();
2935  ++var_iter ) {
2936  if ( ( *var_iter ) != "m13Sq" && ( *var_iter ) != "m23Sq" ) {
2937  ++nVars;
2938  }
2939  }
2940  if ( nVars == 2 ) {
2941  twodimMap.insert(
2942  std::make_pair( bkgndClass, std::make_pair( varNames[0], varNames[1] ) ) );
2943  }
2944  }
2945  }
2946  }
2947 
2948  return twodimMap;
2949 }
2950 
2951 void LauCPFitModel::storePerEvtLlhds()
2952 {
2953  std::cout << "INFO in LauCPFitModel::storePerEvtLlhds : Storing per-event likelihood values..."
2954  << std::endl;
2955 
2956  // if we've not been using the DP model then we need to cache all
2957  // the info here so that we can get the efficiency from it
2958 
2959  LauFitDataTree* inputFitData = this->fitData();
2960 
2961  if ( ! this->useDP() && this->storeDPEff() ) {
2962  negSigModel_->initialise( negCoeffs_ );
2963  posSigModel_->initialise( posCoeffs_ );
2964  negSigModel_->fillDataTree( *inputFitData );
2965  posSigModel_->fillDataTree( *inputFitData );
2966  }
2967 
2968  UInt_t evtsPerExpt( this->eventsPerExpt() );
2969 
2970  LauIsobarDynamics* sigModel( 0 );
2971  LauPdfPList* sigPdfs( 0 );
2972  LauPdfPList* scfPdfs( 0 );
2973  LauBkgndPdfsList* bkgndPdfs( 0 );
2974 
2975  for ( UInt_t iEvt = 0; iEvt < evtsPerExpt; ++iEvt ) {
2976 
2977  this->setSPlotNtupleIntegerBranchValue( "iExpt", this->iExpt() );
2978  this->setSPlotNtupleIntegerBranchValue( "iEvtWithinExpt", iEvt );
2979 
2980  // Find out whether we have B- or B+
2981  if ( tagged_ ) {
2982  const LauFitData& dataValues = inputFitData->getData( iEvt );
2983  LauFitData::const_iterator iter = dataValues.find( "charge" );
2984  curEvtCharge_ = static_cast<Int_t>( iter->second );
2985 
2986  if ( curEvtCharge_ == +1 ) {
2987  sigModel = posSigModel_;
2988  sigPdfs = &posSignalPdfs_;
2989  scfPdfs = &posScfPdfs_;
2990  bkgndPdfs = &posBkgndPdfs_;
2991  } else {
2992  sigModel = negSigModel_;
2993  sigPdfs = &negSignalPdfs_;
2994  scfPdfs = &negScfPdfs_;
2995  bkgndPdfs = &negBkgndPdfs_;
2996  }
2997  } else {
2998  sigPdfs = &negSignalPdfs_;
2999  scfPdfs = &negScfPdfs_;
3000  bkgndPdfs = &negBkgndPdfs_;
3001  }
3002 
3003  // the DP information
3004  this->getEvtDPLikelihood( iEvt );
3005  if ( this->storeDPEff() ) {
3006  if ( ! this->useDP() ) {
3007  posSigModel_->calcLikelihoodInfo( iEvt );
3008  negSigModel_->calcLikelihoodInfo( iEvt );
3009  }
3010  if ( tagged_ ) {
3011  this->setSPlotNtupleDoubleBranchValue( "efficiency", sigModel->getEvtEff() );
3012  if ( negSigModel_->usingScfModel() && posSigModel_->usingScfModel() ) {
3013  this->setSPlotNtupleDoubleBranchValue( "scffraction",
3014  sigModel->getEvtScfFraction() );
3015  }
3016  } else {
3017  this->setSPlotNtupleDoubleBranchValue( "efficiency",
3018  0.5 * ( posSigModel_->getEvtEff() +
3019  negSigModel_->getEvtEff() ) );
3020  if ( negSigModel_->usingScfModel() && posSigModel_->usingScfModel() ) {
3021  this->setSPlotNtupleDoubleBranchValue(
3022  "scffraction",
3023  0.5 * ( posSigModel_->getEvtScfFraction() +
3024  negSigModel_->getEvtScfFraction() ) );
3025  }
3026  }
3027  }
3028  if ( this->useDP() ) {
3029  sigTotalLike_ = sigDPLike_;
3030  if ( useSCF_ ) {
3031  scfTotalLike_ = scfDPLike_;
3032  this->setSPlotNtupleDoubleBranchValue( "sigTMDPLike", sigDPLike_ );
3033  this->setSPlotNtupleDoubleBranchValue( "sigSCFDPLike", scfDPLike_ );
3034  } else {
3035  this->setSPlotNtupleDoubleBranchValue( "sigDPLike", sigDPLike_ );
3036  }
3037  if ( usingBkgnd_ ) {
3038  const UInt_t nBkgnds = this->nBkgndClasses();
3039  for ( UInt_t iBkgnd( 0 ); iBkgnd < nBkgnds; ++iBkgnd ) {
3040  TString name = this->bkgndClassName( iBkgnd );
3041  name += "DPLike";
3042  this->setSPlotNtupleDoubleBranchValue( name, bkgndDPLike_[iBkgnd] );
3043  }
3044  }
3045  } else {
3046  sigTotalLike_ = 1.0;
3047  if ( useSCF_ ) {
3048  scfTotalLike_ = 1.0;
3049  }
3050  if ( usingBkgnd_ ) {
3051  const UInt_t nBkgnds = this->nBkgndClasses();
3052  for ( UInt_t iBkgnd( 0 ); iBkgnd < nBkgnds; ++iBkgnd ) {
3053  bkgndTotalLike_[iBkgnd] = 1.0;
3054  }
3055  }
3056  }
3057 
3058  // the signal PDF values
3059  if ( useSCF_ ) {
3060  sigTotalLike_ *= this->setSPlotNtupleBranchValues( sigPdfs, "sigTM", iEvt );
3061  scfTotalLike_ *= this->setSPlotNtupleBranchValues( scfPdfs, "sigSCF", iEvt );
3062  } else {
3063  sigTotalLike_ *= this->setSPlotNtupleBranchValues( sigPdfs, "sig", iEvt );
3064  }
3065 
3066  // the background PDF values
3067  if ( usingBkgnd_ ) {
3068  const UInt_t nBkgnds = this->nBkgndClasses();
3069  for ( UInt_t iBkgnd( 0 ); iBkgnd < nBkgnds; ++iBkgnd ) {
3070  const TString& bkgndClass = this->bkgndClassName( iBkgnd );
3071  LauPdfPList& pdfs = ( *bkgndPdfs )[iBkgnd];
3072  bkgndTotalLike_[iBkgnd] *=
3073  this->setSPlotNtupleBranchValues( &( pdfs ), bkgndClass, iEvt );
3074  }
3075  }
3076 
3077  // the total likelihoods
3078  if ( useSCF_ ) {
3079  Double_t scfFrac( 0.0 );
3080  if ( useSCFHist_ ) {
3081  scfFrac = recoSCFFracs_[iEvt];
3082  } else {
3083  scfFrac = scfFrac_.unblindValue();
3084  }
3085  this->setSPlotNtupleDoubleBranchValue( "sigSCFFrac", scfFrac );
3086  sigTotalLike_ *= ( 1.0 - scfFrac );
3087  if ( scfMap_ == 0 ) {
3088  scfTotalLike_ *= scfFrac;
3089  }
3090  this->setSPlotNtupleDoubleBranchValue( "sigTMTotalLike", sigTotalLike_ );
3091  this->setSPlotNtupleDoubleBranchValue( "sigSCFTotalLike", scfTotalLike_ );
3092  } else {
3093  this->setSPlotNtupleDoubleBranchValue( "sigTotalLike", sigTotalLike_ );
3094  }
3095  if ( usingBkgnd_ ) {
3096  const UInt_t nBkgnds = this->nBkgndClasses();
3097  for ( UInt_t iBkgnd( 0 ); iBkgnd < nBkgnds; ++iBkgnd ) {
3098  TString name = this->bkgndClassName( iBkgnd );
3099  name += "TotalLike";
3100  this->setSPlotNtupleDoubleBranchValue( name, bkgndTotalLike_[iBkgnd] );
3101  }
3102  }
3103  // fill the tree
3104  this->fillSPlotNtupleBranches();
3105  }
3106  std::cout << "INFO in LauCPFitModel::storePerEvtLlhds : Finished storing per-event likelihood values."
3107  << std::endl;
3108 }
3109 
3110 std::map<TString, LauComplex> LauCPFitModel::getDPAmps( const Double_t m13Sq, const Double_t m23Sq )
3111 {
3112  // Initialise the DP model, if not already done
3113  if ( negCoeffs_.empty() || posCoeffs_.empty() ) {
3114  this->updateCoeffs();
3115  this->initialiseDPModels();
3116  }
3117 
3118  if ( ! posKinematics_->withinDPLimits( m13Sq, m23Sq ) ) {
3119  return {
3120  {"signal_" + negParent_, {}},
3121  {"signal_" + posParent_, {}}
3122  };
3123  }
3124 
3125  negSigModel_->calcLikelihoodInfo( m13Sq, m23Sq );
3126  posSigModel_->calcLikelihoodInfo( m13Sq, m23Sq );
3127 
3128  return {
3129  {"signal_" + negParent_, negSigModel_->getEvtDPAmp()},
3130  {"signal_" + posParent_, posSigModel_->getEvtDPAmp()}
3131  };
3132 }
3133 
3134 std::map<TString, Double_t> LauCPFitModel::getDPLikelihoods( const Double_t m13Sq,
3135  const Double_t m23Sq )
3136 {
3137  // Initialise the DP model, if not already done
3138  if ( negCoeffs_.empty() || posCoeffs_.empty() ) {
3139  this->updateCoeffs();
3140  this->initialiseDPModels();
3141  }
3142 
3143  std::map<TString, Double_t> likelihoods;
3144 
3145  if ( ! posKinematics_->withinDPLimits( m13Sq, m23Sq ) ) {
3146  likelihoods.emplace( "signal_" + negParent_, 0.0 );
3147  likelihoods.emplace( "signal_" + posParent_, 0.0 );
3148  if ( usingBkgnd_ ) {
3149  const UInt_t nBkgnds { this->nBkgndClasses() };
3150  for ( UInt_t bkgndID( 0 ); bkgndID < nBkgnds; ++bkgndID ) {
3151  likelihoods.emplace( this->bkgndClassName( bkgndID ) + "_" + negParent_, 0.0 );
3152  likelihoods.emplace( this->bkgndClassName( bkgndID ) + "_" + posParent_, 0.0 );
3153  }
3154  }
3155  return likelihoods;
3156  }
3157 
3158  negSigModel_->calcLikelihoodInfo( m13Sq, m23Sq );
3159  posSigModel_->calcLikelihoodInfo( m13Sq, m23Sq );
3160 
3161  // See comment in getEvtDPLikelihood for explanation of the 2.0 factor
3162  const Double_t norm { 2.0 / ( negSigModel_->getDPNorm() + posSigModel_->getDPNorm() ) };
3163 
3164  likelihoods.emplace( "signal_" + negParent_, negSigModel_->getEvtIntensity() * norm );
3165  likelihoods.emplace( "signal_" + posParent_, posSigModel_->getEvtIntensity() * norm );
3166 
3167  // TODO - SCF signal
3168  static bool warningIssued { false };
3169  if ( useSCF_ && ! warningIssued ) {
3170  warningIssued = true;
3171  std::cerr << "WARNING in LauCPFitModel::getDPLikelihoods : calculation of SCF likelihoods not currently implemented in this function\n";
3172  std::cerr << " : signal likelihood will just be the truth-matched value";
3173  std::cerr << std::endl;
3174  }
3175 
3176  if ( usingBkgnd_ ) {
3177  const UInt_t nBkgnds { this->nBkgndClasses() };
3178  for ( UInt_t bkgndID( 0 ); bkgndID < nBkgnds; ++bkgndID ) {
3179  likelihoods.emplace( this->bkgndClassName( bkgndID ) + "_" + negParent_,
3180  negBkgndDPModels_[bkgndID]->getLikelihood( m13Sq, m23Sq ) );
3181  likelihoods.emplace( this->bkgndClassName( bkgndID ) + "_" + posParent_,
3182  posBkgndDPModels_[bkgndID]->getLikelihood( m13Sq, m23Sq ) );
3183  }
3184  }
3185 
3186  return likelihoods;
3187 }
3188 
3189 void LauCPFitModel::embedNegSignal( const TString& fileName,
3190  const TString& treeName,
3191  Bool_t reuseEventsWithinEnsemble,
3192  Bool_t reuseEventsWithinExperiment,
3193  Bool_t useReweighting )
3194 {
3195  if ( negSignalTree_ ) {
3196  std::cerr << "ERROR in LauCPFitModel::embedNegSignal : Already embedding signal from a file."
3197  << std::endl;
3198  return;
3199  }
3200 
3201  negSignalTree_ = new LauEmbeddedData( fileName, treeName, reuseEventsWithinExperiment );
3202  Bool_t dataOK = negSignalTree_->findBranches();
3203  if ( ! dataOK ) {
3204  delete negSignalTree_;
3205  negSignalTree_ = 0;
3206  std::cerr << "ERROR in LauCPFitModel::embedNegSignal : Problem creating data tree for embedding."
3207  << std::endl;
3208  return;
3209  }
3210  reuseSignal_ = reuseEventsWithinEnsemble;
3211  useNegReweighting_ = useReweighting;
3212  if ( this->enableEmbedding() == kFALSE ) {
3213  this->enableEmbedding( kTRUE );
3214  }
3215 }
3216 
3217 void LauCPFitModel::embedNegBkgnd( const TString& bkgndClass,
3218  const TString& fileName,
3219  const TString& treeName,
3220  Bool_t reuseEventsWithinEnsemble,
3221  Bool_t reuseEventsWithinExperiment )
3222 {
3223  if ( ! this->validBkgndClass( bkgndClass ) ) {
3224  std::cerr << "ERROR in LauCPFitModel::embedBkgnd : Invalid background class \""
3225  << bkgndClass << "\"." << std::endl;
3226  std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed."
3227  << std::endl;
3228  return;
3229  }
3230 
3231  UInt_t bkgndID = this->bkgndClassID( bkgndClass );
3232 
3233  if ( negBkgndTree_[bkgndID] ) {
3234  std::cerr << "ERROR in LauCPFitModel::embedNegBkgnd : Already embedding background from a file."
3235  << std::endl;
3236  return;
3237  }
3238 
3239  negBkgndTree_[bkgndID] = new LauEmbeddedData( fileName, treeName, reuseEventsWithinExperiment );
3240  Bool_t dataOK = negBkgndTree_[bkgndID]->findBranches();
3241  if ( ! dataOK ) {
3242  delete negBkgndTree_[bkgndID];
3243  negBkgndTree_[bkgndID] = 0;
3244  std::cerr << "ERROR in LauCPFitModel::embedNegBkgnd : Problem creating data tree for embedding."
3245  << std::endl;
3246  return;
3247  }
3248  reuseBkgnd_[bkgndID] = reuseEventsWithinEnsemble;
3249  if ( this->enableEmbedding() == kFALSE ) {
3250  this->enableEmbedding( kTRUE );
3251  }
3252 }
3253 
3254 void LauCPFitModel::embedPosSignal( const TString& fileName,
3255  const TString& treeName,
3256  Bool_t reuseEventsWithinEnsemble,
3257  Bool_t reuseEventsWithinExperiment,
3258  Bool_t useReweighting )
3259 {
3260  if ( posSignalTree_ ) {
3261  std::cerr << "ERROR in LauCPFitModel::embedPosSignal : Already embedding signal from a file."
3262  << std::endl;
3263  return;
3264  }
3265 
3266  posSignalTree_ = new LauEmbeddedData( fileName, treeName, reuseEventsWithinExperiment );
3267  Bool_t dataOK = posSignalTree_->findBranches();
3268  if ( ! dataOK ) {
3269  delete posSignalTree_;
3270  posSignalTree_ = 0;
3271  std::cerr << "ERROR in LauCPFitModel::embedPosSignal : Problem creating data tree for embedding."
3272  << std::endl;
3273  return;
3274  }
3275  reuseSignal_ = reuseEventsWithinEnsemble;
3276  usePosReweighting_ = useReweighting;
3277  if ( this->enableEmbedding() == kFALSE ) {
3278  this->enableEmbedding( kTRUE );
3279  }
3280 }
3281 
3282 void LauCPFitModel::embedPosBkgnd( const TString& bkgndClass,
3283  const TString& fileName,
3284  const TString& treeName,
3285  Bool_t reuseEventsWithinEnsemble,
3286  Bool_t reuseEventsWithinExperiment )
3287 {
3288  if ( ! this->validBkgndClass( bkgndClass ) ) {
3289  std::cerr << "ERROR in LauCPFitModel::embedBkgnd : Invalid background class \""
3290  << bkgndClass << "\"." << std::endl;
3291  std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed."
3292  << std::endl;
3293  return;
3294  }
3295 
3296  UInt_t bkgndID = this->bkgndClassID( bkgndClass );
3297 
3298  if ( posBkgndTree_[bkgndID] ) {
3299  std::cerr << "ERROR in LauCPFitModel::embedPosBkgnd : Already embedding background from a file."
3300  << std::endl;
3301  return;
3302  }
3303 
3304  posBkgndTree_[bkgndID] = new LauEmbeddedData( fileName, treeName, reuseEventsWithinExperiment );
3305  Bool_t dataOK = posBkgndTree_[bkgndID]->findBranches();
3306  if ( ! dataOK ) {
3307  delete posBkgndTree_[bkgndID];
3308  posBkgndTree_[bkgndID] = 0;
3309  std::cerr << "ERROR in LauCPFitModel::embedPosBkgnd : Problem creating data tree for embedding."
3310  << std::endl;
3311  return;
3312  }
3313  reuseBkgnd_[bkgndID] = reuseEventsWithinEnsemble;
3314  if ( this->enableEmbedding() == kFALSE ) {
3315  this->enableEmbedding( kTRUE );
3316  }
3317 }
3318 
3319 void LauCPFitModel::weightEvents( const TString& dataFileName, const TString& dataTreeName )
3320 {
3321  // Routine to provide weights for events that are uniformly distributed
3322  // in the DP (or square DP) so as to reproduce the given DP model
3323 
3324  const Bool_t squareDP { posKinematics_->squareDP() || negKinematics_->squareDP() };
3325  if ( squareDP ) {
3326  std::cout << "INFO in LauCPFitModel::weightEvents : will store DP model weights and the square DP jacobian\n";
3327  std::cout << " : the DP model weights can be used on their own to weight events that were generated flat in phase space\n";
3328  std::cout << " : or they can be multiplied by the jacobian to weight events that were generated flat in the square DP\n";
3329  std::cout << " : or they can be multiplied by max(1.0, jacobian) to weight events that were generated quasi-flat in the square DP"
3330  << std::endl;
3331  } else {
3332  std::cout << "INFO in LauCPFitModel::weightEvents : will store DP model weights suitable for weighting events that were generated flat in phase space"
3333  << std::endl;
3334  }
3335 
3336  // This reads in the given dataFile and creates an input
3337  // fit data tree that stores them for all events and experiments.
3338  const Bool_t dataOK { this->verifyFitData( dataFileName, dataTreeName ) };
3339  if ( ! dataOK ) {
3340  std::cerr << "ERROR in LauCPFitModel::weightEvents : Problem caching the data." << std::endl;
3341  return;
3342  }
3343 
3344  LauFitDataTree* inputFitData { this->fitData() };
3345 
3346  if ( ! inputFitData->haveBranch( "m13Sq_MC" ) || ! inputFitData->haveBranch( "m23Sq_MC" ) ) {
3347  std::cerr << "WARNING in LauCPFitModel::weightEvents : Cannot find MC truth DP coordinate branches in supplied data, aborting."
3348  << std::endl;
3349  return;
3350  }
3351  if ( ! inputFitData->haveBranch( "charge" ) ) {
3352  std::cerr << "WARNING in LauCPFitModel::weightEvents : Cannot find branch specifying event charge in supplied data, aborting."
3353  << std::endl;
3354  return;
3355  }
3356 
3357  // Create the ntuple to hold the DP weights
3358  TString weightsFileName { dataFileName };
3359  const Ssiz_t index { weightsFileName.Last( '.' ) };
3360  weightsFileName.Insert( index, "_DPweights" );
3361 
3362  LauGenNtuple weightsTuple { weightsFileName, dataTreeName };
3363  weightsTuple.addIntegerBranch( "iExpt" );
3364  weightsTuple.addIntegerBranch( "iEvtWithinExpt" );
3365  weightsTuple.addDoubleBranch( "dpModelWeight" );
3366  if ( squareDP ) {
3367  weightsTuple.addDoubleBranch( "sqDPJacobian" );
3368  }
3369 
3370  const UInt_t nExpmt { this->nExpt() };
3371  const UInt_t firstExpmt { this->firstExpt() };
3372  for ( UInt_t iExpmt { firstExpmt }; iExpmt < ( firstExpmt + nExpmt ); ++iExpmt ) {
3373 
3374  inputFitData->readExperimentData( iExpmt );
3375 
3376  const UInt_t nEvents { inputFitData->nEvents() };
3377  if ( nEvents < 1 ) {
3378  std::cerr << "WARNING in LauCPFitModel::weightEvents : Zero events in experiment "
3379  << iExpmt << ", skipping..." << std::endl;
3380  continue;
3381  }
3382 
3383  weightsTuple.setIntegerBranchValue( "iExpt", iExpmt );
3384 
3385  // Calculate and store the weights for the events in this experiment
3386  for ( UInt_t iEvent { 0 }; iEvent < nEvents; ++iEvent ) {
3387 
3388  weightsTuple.setIntegerBranchValue( "iEvtWithinExpt", iEvent );
3389 
3390  const LauFitData& evtData = inputFitData->getData( iEvent );
3391 
3392  const auto m13Sq_MC { evtData.at( "m13Sq_MC" ) };
3393  const auto m23Sq_MC { evtData.at( "m23Sq_MC" ) };
3394  const auto charge { static_cast<Int_t>( evtData.at( "charge" ) ) };
3395 
3396  Double_t dpModelWeight { 0.0 };
3397  Double_t jacobian { 1.0 };
3398 
3399  LauKinematics* kinematics { nullptr };
3400  LauIsobarDynamics* dpModel { nullptr };
3401 
3402  if ( charge > 0 ) {
3403  kinematics = posKinematics_;
3404  dpModel = posSigModel_;
3405  } else {
3406  kinematics = negKinematics_;
3407  dpModel = negSigModel_;
3408  }
3409 
3410  if ( kinematics->withinDPLimits( m13Sq_MC, m23Sq_MC ) ) {
3411 
3412  kinematics->updateKinematics( m13Sq_MC, m23Sq_MC );
3413  dpModelWeight = dpModel->getEventWeight();
3414 
3415  if ( squareDP ) {
3416  jacobian = kinematics->calcSqDPJacobian();
3417  }
3418 
3419  const Double_t norm { 0.5 *
3420  ( negSigModel_->getDPNorm() + posSigModel_->getDPNorm() ) };
3421  dpModelWeight /= norm;
3422  }
3423 
3424  weightsTuple.setDoubleBranchValue( "dpModelWeight", dpModelWeight );
3425  if ( squareDP ) {
3426  weightsTuple.setDoubleBranchValue( "sqDPJacobian", jacobian );
3427  }
3428  weightsTuple.fillBranches();
3429  }
3430  }
3431 
3432  weightsTuple.buildIndex( "iExpt", "iEvtWithinExpt" );
3433  weightsTuple.addFriendTree( dataFileName, dataTreeName );
3434  weightsTuple.writeOutGenResults();
3435 }
3436 
3437 void LauCPFitModel::savePDFPlots( const TString& label )
3438 {
3439  savePDFPlotsWave( label, 0 );
3440  savePDFPlotsWave( label, 1 );
3441  savePDFPlotsWave( label, 2 );
3442 
3443  std::cout << "LauCPFitModel::plot" << std::endl;
3444  // ((LauIsobarDynamics*)negSigModel_)->plot();
3445 
3446  //Double_t minm13 = negSigModel_->getKinematics()->getm13Min();
3447  Double_t minm13 = 0.0;
3448  Double_t maxm13 = negSigModel_->getKinematics()->getm13Max();
3449  //Double_t minm23 = negSigModel_->getKinematics()->getm23Min();
3450  Double_t minm23 = 0.0;
3451  Double_t maxm23 = negSigModel_->getKinematics()->getm23Max();
3452 
3453  Double_t mins13 = minm13 * minm13;
3454  Double_t maxs13 = maxm13 * maxm13;
3455  Double_t mins23 = minm23 * minm23;
3456  Double_t maxs23 = maxm23 * maxm23;
3457 
3458  Double_t s13, s23, posChPdf, negChPdf;
3459 
3460  TString xLabel = "s13";
3461  TString yLabel = "s23";
3462 
3463  if ( negSigModel_->getDaughters()->gotSymmetricalDP() ) {
3464  xLabel = "sHigh";
3465  yLabel = "sLow";
3466  }
3467 
3468  Int_t n13 = 200.00, n23 = 200.00;
3469  Double_t delta13, delta23;
3470  delta13 = ( maxs13 - mins13 ) / n13;
3471  delta23 = ( maxs23 - mins23 ) / n23;
3472  UInt_t nAmp = negSigModel_->getnCohAmp();
3473  for ( UInt_t resID = 0; resID <= nAmp; ++resID ) {
3474  TGraph2D* posDt = new TGraph2D();
3475  TGraph2D* negDt = new TGraph2D();
3476  TGraph2D* acpDt = new TGraph2D();
3477 
3478  TString resName = "TotalAmp";
3479  if ( resID != nAmp ) {
3480  TString tStrResID = Form( "%d", resID );
3481  const LauIsobarDynamics* model = negSigModel_;
3482  const LauAbsResonance* resonance = model->getResonance( resID );
3483  resName = resonance->getResonanceName();
3484  std::cout << "resName = " << resName << std::endl;
3485  }
3486 
3487  resName.ReplaceAll( "(", "" );
3488  resName.ReplaceAll( ")", "" );
3489  resName.ReplaceAll( "*", "Star" );
3490 
3491  posDt->SetName( resName + label );
3492  posDt->SetTitle( resName + " (" + label + ") Positive" );
3493  negDt->SetName( resName + label );
3494  negDt->SetTitle( resName + " (" + label + ") Negative" );
3495  acpDt->SetName( resName + label );
3496  acpDt->SetTitle( resName + " (" + label + ") Asymmetry" );
3497 
3498  Int_t count = 0;
3499  for ( Int_t i = 0; i < n13; i++ ) {
3500  s13 = mins13 + i * delta13;
3501  for ( Int_t j = 0; j < n23; j++ ) {
3502  s23 = mins23 + j * delta23;
3503  if ( negSigModel_->getKinematics()->withinDPLimits2( s23, s13 ) ) {
3504  if ( negSigModel_->getDaughters()->gotSymmetricalDP() && ( s13 > s23 ) )
3505  continue;
3506 
3507  negSigModel_->calcLikelihoodInfo( s13, s23 );
3508  posSigModel_->calcLikelihoodInfo( s13, s23 );
3509 
3510  LauComplex negChAmp = negSigModel_->getEvtDPAmp();
3511  LauComplex posChAmp = posSigModel_->getEvtDPAmp();
3512 
3513  if ( resID != nAmp ) {
3514  negChAmp = negSigModel_->getFullAmplitude( resID );
3515  posChAmp = posSigModel_->getFullAmplitude( resID );
3516  }
3517  negChPdf = negChAmp.abs2();
3518  posChPdf = posChAmp.abs2();
3519  negDt->SetPoint( count, s23, s13, negChPdf ); // s23=sHigh, s13 = sLow
3520  posDt->SetPoint( count, s23, s13, posChPdf ); // s23=sHigh, s13 = sLow
3521  acpDt->SetPoint( count, s23, s13, negChPdf - posChPdf ); // s23=sHigh, s13 = sLow
3522  count++;
3523  }
3524  }
3525  }
3526  gStyle->SetPalette( 1 );
3527  TCanvas* posC = new TCanvas( "c" + resName + label + "Positive",
3528  resName + " (" + label + ") Positive",
3529  0,
3530  0,
3531  600,
3532  400 );
3533  posDt->GetXaxis()->SetTitle( xLabel );
3534  posDt->GetYaxis()->SetTitle( yLabel );
3535  posDt->Draw( "SURF1" );
3536  posDt->GetXaxis()->SetTitle( xLabel );
3537  posDt->GetYaxis()->SetTitle( yLabel );
3538  posC->SaveAs( "plot_2D_" + resName + "_" + label + "Positive.C" );
3539 
3540  TCanvas* negC = new TCanvas( "c" + resName + label + "Negative",
3541  resName + " (" + label + ") Negative",
3542  0,
3543  0,
3544  600,
3545  400 );
3546  negDt->GetXaxis()->SetTitle( xLabel );
3547  negDt->GetYaxis()->SetTitle( yLabel );
3548  negDt->Draw( "SURF1" );
3549  negDt->GetXaxis()->SetTitle( xLabel );
3550  negDt->GetYaxis()->SetTitle( yLabel );
3551  negC->SaveAs( "plot_2D_" + resName + "_" + label + "Negative.C" );
3552 
3553  TCanvas* acpC = new TCanvas( "c" + resName + label + "Asymmetry",
3554  resName + " (" + label + ") Asymmetry",
3555  0,
3556  0,
3557  600,
3558  400 );
3559  acpDt->GetXaxis()->SetTitle( xLabel );
3560  acpDt->GetYaxis()->SetTitle( yLabel );
3561  acpDt->Draw( "SURF1" );
3562  acpDt->GetXaxis()->SetTitle( xLabel );
3563  acpDt->GetYaxis()->SetTitle( yLabel );
3564  acpC->SaveAs( "plot_2D_" + resName + "_" + label + "Asymmetry.C" );
3565  }
3566 }
3567 
3568 void LauCPFitModel::savePDFPlotsWave( const TString& label, const Int_t& spin )
3569 {
3570 
3571  std::cout << "label = " << label << ", spin = " << spin << std::endl;
3572 
3573  TString tStrResID = "S_Wave";
3574  if ( spin == 1 )
3575  tStrResID = "P_Wave";
3576  if ( spin == 2 )
3577  tStrResID = "D_Wave";
3578 
3579  TString xLabel = "s13";
3580  TString yLabel = "s23";
3581 
3582  std::cout << "LauCPFitModel::savePDFPlotsWave: " << tStrResID << std::endl;
3583 
3584  Double_t minm13 = 0.0;
3585  Double_t maxm13 = negSigModel_->getKinematics()->getm13Max();
3586  Double_t minm23 = 0.0;
3587  Double_t maxm23 = negSigModel_->getKinematics()->getm23Max();
3588 
3589  Double_t mins13 = minm13 * minm13;
3590  Double_t maxs13 = maxm13 * maxm13;
3591  Double_t mins23 = minm23 * minm23;
3592  Double_t maxs23 = maxm23 * maxm23;
3593 
3594  Double_t s13, s23, posChPdf, negChPdf;
3595  TGraph2D* posDt = new TGraph2D();
3596  TGraph2D* negDt = new TGraph2D();
3597  TGraph2D* acpDt = new TGraph2D();
3598 
3599  posDt->SetName( tStrResID + label );
3600  posDt->SetTitle( tStrResID + " (" + label + ") Positive" );
3601  negDt->SetName( tStrResID + label );
3602  negDt->SetTitle( tStrResID + " (" + label + ") Negative" );
3603  acpDt->SetName( tStrResID + label );
3604  acpDt->SetTitle( tStrResID + " (" + label + ") Asymmetry" );
3605 
3606  Int_t n13 = 200.00, n23 = 200.00;
3607  Double_t delta13, delta23;
3608  delta13 = ( maxs13 - mins13 ) / n13;
3609  delta23 = ( maxs23 - mins23 ) / n23;
3610  UInt_t nAmp = negSigModel_->getnCohAmp();
3611 
3612  Int_t count = 0;
3613  for ( Int_t i = 0; i < n13; i++ ) {
3614  s13 = mins13 + i * delta13;
3615  for ( Int_t j = 0; j < n23; j++ ) {
3616  s23 = mins23 + j * delta23;
3617  if ( negSigModel_->getKinematics()->withinDPLimits2( s23, s13 ) ) {
3618  if ( negSigModel_->getDaughters()->gotSymmetricalDP() && ( s13 > s23 ) )
3619  continue;
3620 
3621  LauComplex negChAmp( 0, 0 );
3622  LauComplex posChAmp( 0, 0 );
3623  Bool_t noWaveRes = kTRUE;
3624  negSigModel_->calcLikelihoodInfo( s13, s23 );
3625  for ( UInt_t resID = 0; resID < nAmp; ++resID ) {
3626  const LauIsobarDynamics* model = negSigModel_;
3627  const LauAbsResonance* resonance = model->getResonance( resID );
3628  Int_t spin_res = resonance->getSpin();
3629  if ( spin != spin_res )
3630  continue;
3631  noWaveRes = kFALSE;
3632  negChAmp += negSigModel_->getFullAmplitude( resID );
3633  posChAmp += posSigModel_->getFullAmplitude( resID );
3634  }
3635 
3636  if ( noWaveRes )
3637  return;
3638 
3639  negChPdf = negChAmp.abs2();
3640  posChPdf = posChAmp.abs2();
3641 
3642  negDt->SetPoint( count, s23, s13, negChPdf ); // s23=sHigh, s13 = sLow
3643  posDt->SetPoint( count, s23, s13, posChPdf ); // s23=sHigh, s13 = sLow
3644  acpDt->SetPoint( count, s23, s13, negChPdf - posChPdf ); // s23=sHigh, s13 = sLow
3645  count++;
3646  }
3647  }
3648  }
3649  gStyle->SetPalette( 1 );
3650  TCanvas* posC = new TCanvas( "c" + tStrResID + label + "Positive",
3651  tStrResID + " (" + label + ") Positive",
3652  0,
3653  0,
3654  600,
3655  400 );
3656  posDt->GetXaxis()->SetTitle( xLabel );
3657  posDt->GetYaxis()->SetTitle( yLabel );
3658  posDt->Draw( "SURF1" );
3659  posDt->GetXaxis()->SetTitle( xLabel );
3660  posDt->GetYaxis()->SetTitle( yLabel );
3661  posC->SaveAs( "plot_2D_" + tStrResID + "_" + label + "Positive.C" );
3662 
3663  TCanvas* negC = new TCanvas( "c" + tStrResID + label + "Negative",
3664  tStrResID + " (" + label + ") Negative",
3665  0,
3666  0,
3667  600,
3668  400 );
3669  negDt->GetXaxis()->SetTitle( xLabel );
3670  negDt->GetYaxis()->SetTitle( yLabel );
3671  negDt->Draw( "SURF1" );
3672  negDt->GetXaxis()->SetTitle( xLabel );
3673  negDt->GetYaxis()->SetTitle( yLabel );
3674  negC->SaveAs( "plot_2D_" + tStrResID + "_" + label + "Negative.C" );
3675 
3676  TCanvas* acpC = new TCanvas( "c" + tStrResID + label + "Asymmetry",
3677  tStrResID + " (" + label + ") Asymmetry",
3678  0,
3679  0,
3680  600,
3681  400 );
3682  acpDt->GetXaxis()->SetTitle( xLabel );
3683  acpDt->GetYaxis()->SetTitle( yLabel );
3684  acpDt->Draw( "SURF1" );
3685  acpDt->GetXaxis()->SetTitle( xLabel );
3686  acpDt->GetYaxis()->SetTitle( yLabel );
3687  acpC->SaveAs( "plot_2D_" + tStrResID + "_" + label + "Asymmetry.C" );
3688 }
3689 
3690 Double_t LauCPFitModel::getParamFromTree( TTree& tree, const TString& name )
3691 {
3692  TBranch* branch { tree.FindBranch( name ) };
3693  if ( branch ) {
3694  TLeaf* leaf { branch->GetLeaf( name ) };
3695  if ( leaf ) {
3696  tree.GetEntry( 0 );
3697  return leaf->GetValue();
3698  } else {
3699  std::cout << "ERROR in LauCPFitModel::getParamFromTree : Leaf name " + name +
3700  " not found in parameter file!"
3701  << std::endl;
3702  }
3703  } else {
3704  std::cout << "ERROR in LauCPFitModel::getParamFromTree : Branch name " + name +
3705  " not found in parameter file!"
3706  << std::endl;
3707  }
3708 
3709  return -1.1;
3710 }
3711 
3712 void LauCPFitModel::fixParam( LauParameter* param, const Double_t val, const Bool_t fix )
3713 {
3714  std::cout << "INFO in LauCPFitModel::fixParam : Setting " << param->name() << " to " << val
3715  << std::endl;
3716 
3717  param->value( val );
3718  param->genValue( val );
3719  param->initValue( val );
3720 
3721  if ( fix ) {
3722  param->fixed( kTRUE );
3723  } else if ( ! param->fixed() ) {
3724  // Add parameter name to list to indicate that this should not be randomised by randomiseInitFitPars
3725  // (otherwise only those that are fixed are not randomised).
3726  // This is only done to those that are not already fixed (see randomiseInitFitPars).
3727  allImportedFreeParams_.insert( param );
3728  }
3729 }
3730 
3731 void LauCPFitModel::fixParams( std::vector<LauParameter*>& params )
3732 {
3733  const Bool_t fix { fixParams_ };
3734 
3735  // TODO: Allow some parameters to be fixed and some to remain floating (but initialised)
3736 
3737  if ( ! fixParamFileName_.IsNull() ) {
3738 
3739  // Take param values from a file
3740 
3741  TFile* paramFile = TFile::Open( fixParamFileName_, "READ" );
3742  if ( ! paramFile ) {
3743  std::cerr << "ERROR in LauCPFitModel::fixParams : File '" + fixParamFileName_ +
3744  "' could not be opened for reading!"
3745  << std::endl;
3746  return;
3747  }
3748 
3749  TTree* paramTree = dynamic_cast<TTree*>( paramFile->Get( fixParamTreeName_ ) );
3750  if ( ! paramTree ) {
3751  std::cerr << "ERROR in LauCPFitModel::fixParams : Tree '" + fixParamTreeName_ +
3752  "' not found in parameter file!"
3753  << std::endl;
3754  return;
3755  }
3756 
3757  if ( ! fixParamNames_.empty() ) {
3758 
3759  // Fix params from file, according to vector of names
3760 
3761  for ( auto itr = params.begin(); itr != params.end(); ++itr ) {
3762  auto itrName = fixParamNames_.find( ( *itr )->name() );
3763  if ( itrName != fixParamNames_.end() ) {
3764  this->fixParam( *itr, this->getParamFromTree( *paramTree, *itrName ), fix );
3765  }
3766  }
3767 
3768  } else {
3769 
3770  // Fix some (unspecified) parameters from file, prioritising the map (if it exists)
3771 
3772  for ( auto itr = params.begin(); itr != params.end(); ++itr ) {
3773 
3774  const TString& name = ( *itr )->name();
3775 
3776  if ( ! fixParamMap_.empty() ) {
3777  auto nameValItr = fixParamMap_.find( name );
3778  if ( nameValItr != fixParamMap_.end() ) {
3779  this->fixParam( *itr, nameValItr->second, fix );
3780  }
3781  } else {
3782  this->fixParam( *itr, this->getParamFromTree( *paramTree, name ), fix );
3783  }
3784  }
3785 
3786  } // Vector of names?
3787 
3788  } else {
3789 
3790  // Fix param names fom map, leave others floating
3791 
3792  for ( auto itr = params.begin(); itr != params.end(); ++itr ) {
3793  auto nameValItr = this->fixParamMap_.find( ( *itr )->name() );
3794  if ( nameValItr != this->fixParamMap_.end() ) {
3795  this->fixParam( *itr, nameValItr->second, fix );
3796  }
3797  }
3798  }
3799 }
LauCPFitModel::posParent_
TString posParent_
Name of the parent particle.
Definition: LauCPFitModel.hh:686
LauCPFitModel::getTotEvtLikelihood
virtual Double_t getTotEvtLikelihood(UInt_t iEvt)
Get the total likelihood for each event.
Definition: LauCPFitModel.cc:2328
LauAbsCoeffSet::printParValues
virtual void printParValues() const =0
Print the current values of the parameters.
LauEffModel.hh
File containing declaration of LauEffModel class.
LauAbsPdf.hh
File containing declaration of LauAbsPdf class.
LauCPFitModel::getEvtDPLikelihood
virtual void getEvtDPLikelihood(UInt_t iEvt)
Calculate the signal and background likelihoods for the DP for a given event.
Definition: LauCPFitModel.cc:2411
LauCPFitModel::setDPBranchValues
void setDPBranchValues()
Store all of the DP information.
Definition: LauCPFitModel.cc:2050
LauIsobarDynamics::getEventWeight
Double_t getEventWeight()
Calculate the acceptance rate, for events with the current kinematics, when generating events accordi...
Definition: LauIsobarDynamics.cc:2940
LauCPFitModel::bkgndAsym_
LauBkgndYieldList bkgndAsym_
Background asymmetries(s)
Definition: LauCPFitModel.hh:638
LauCPFitModel::scfTotalLike_
Double_t scfTotalLike_
Total SCF likelihood.
Definition: LauCPFitModel.hh:743
LauRandom.hh
File containing LauRandom namespace.
LauAbsFitModel::doEMLFit
Bool_t doEMLFit() const
Determine whether an extended maximum likelihood fit it being performed.
Definition: LauAbsFitModel.hh:117
LauCPFitModel::setExtraNtupleVars
void setExtraNtupleVars()
Set-up other parameters that are derived from the fit results, e.g. fit fractions.
Definition: LauCPFitModel.cc:857
LauIsobarDynamics::resonanceIndex
Int_t resonanceIndex(const TString &resName) const
Retrieve the index for the given resonance.
Definition: LauIsobarDynamics.cc:999
LauIsobarDynamics::getResonance
const LauAbsResonance * getResonance(const UInt_t resIndex) const
Retrieve a resonance by its index.
Definition: LauIsobarDynamics.cc:1043
LauCPFitModel::embedPosSignal
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.
Definition: LauCPFitModel.cc:3254
LauAbsFitModel::addResonanceParameters
UInt_t addResonanceParameters(LauParameter *param)
Add the given parameter to the list of resonance parameters and the list of all fit parameters.
Definition: LauAbsFitModel.cc:994
LauAbsCoeffSet
Class for defining the abstract interface for complex coefficient classes.
Definition: LauAbsCoeffSet.hh:48
LauParameter
Class for defining the fit parameter objects.
Definition: LauParameter.hh:49
LauAbsResonance::getResonanceName
const TString & getResonanceName() const
Get the name of the resonance.
Definition: LauAbsResonance.hh:158
LauScfMap::trueHist
TH2 * trueHist(Int_t trueBin)
Retrieve the migration histogram for trueBin.
Definition: LauScfMap.cc:190
LauParameter::unblindValue
Double_t unblindValue() const
The unblinded value of the parameter.
Definition: LauParameter.hh:178
LauKinematics::getc13
Double_t getc13() const
Get the cosine of the helicity angle theta13.
Definition: LauKinematics.hh:275
LauIsobarDynamics::getFitFractionsEfficiencyUncorrected
const LauParArray & getFitFractionsEfficiencyUncorrected() const
Retrieve the fit fractions for the amplitude components.
Definition: LauIsobarDynamics.hh:507
LauCPFitModel::scfFracHist_
LauEffModel * scfFracHist_
The histogram giving the DP-dependence of the SCF fraction.
Definition: LauCPFitModel.hh:662
LauAbsFitModel::doPoissonSmearing
Bool_t doPoissonSmearing() const
Determine whether Poisson smearing is enabled for the toy MC generation.
Definition: LauAbsFitModel.hh:126
LauKinematics::getc23
Double_t getc23() const
Get the cosine of the helicity angle theta23.
Definition: LauKinematics.hh:267
LauScfMap::listBinCentres
void listBinCentres(std::vector< Double_t > &xCoords, std::vector< Double_t > &yCoords) const
Create lists of the co-ordinates of the centres of true bins.
Definition: LauScfMap.cc:128
LauIsobarDynamics::calcLikelihoodInfo
void calcLikelihoodInfo(const UInt_t iEvt)
Calculate the likelihood (and all associated information) for the given event number.
Definition: LauIsobarDynamics.cc:2725
LauCPFitModel::storePerEvtLlhds
virtual void storePerEvtLlhds()
Store the per event likelihood values.
Definition: LauCPFitModel.cc:2951
LauParameter::value
Double_t value() const
The value of the parameter.
Definition: LauParameter.hh:172
LauIsobarDynamics::getDPNorm
Double_t getDPNorm() const
Retrieve the normalisation factor for the log-likelihood function.
Definition: LauIsobarDynamics.hh:534
LauCPFitModel::cacheInputFitVars
virtual void cacheInputFitVars()
Read in the input fit data variables, e.g. m13Sq and m23Sq.
Definition: LauCPFitModel.cc:2188
LauAbsFitModel::setGenNtupleIntegerBranchValue
virtual void setGenNtupleIntegerBranchValue(const TString &name, Int_t value)
Set the value of an integer branch in the gen tree.
Definition: LauAbsFitModel.cc:382
LauCPFitModel::posBkgndPdfs_
LauBkgndPdfsList posBkgndPdfs_
The B+ background PDFs.
Definition: LauCPFitModel.hh:569
LauCPFitModel::genExpt
virtual Bool_t genExpt()
Toy MC generation and fitting overloaded functions.
Definition: LauCPFitModel.cc:1624
LauEmbeddedData::getValue
Double_t getValue(const TString &name) const
Get the value of a specified branch.
Definition: LauEmbeddedData.cc:193
LauCPFitModel::printAsymmetries
virtual void printAsymmetries(std::ostream &output)
Print the asymmetries.
Definition: LauCPFitModel.cc:1344
LauCPFitModel::setupSPlotNtupleBranches
virtual void setupSPlotNtupleBranches()
Add branches to store experiment number and the event number within the experiment.
Definition: LauCPFitModel.cc:2618
LauKinematics::getm12Sq
Double_t getm12Sq() const
Get the m12 invariant mass square.
Definition: LauKinematics.hh:239
LauIsobarDynamics::getEvtEff
Double_t getEvtEff() const
Retrieve the efficiency for the current event.
Definition: LauIsobarDynamics.hh:355
LauCPFitModel::negKinematics_
LauKinematics * negKinematics_
The B- Dalitz plot kinematics object.
Definition: LauCPFitModel.hh:548
LauCPFitModel::getEvtExtraLikelihoods
virtual void getEvtExtraLikelihoods(UInt_t iEvt)
Determine the signal and background likelihood for the extra variables for a given event.
Definition: LauCPFitModel.cc:2566
LauAbsFitModel::verifyFitData
virtual Bool_t verifyFitData(const TString &dataFileName, const TString &dataTreeName)
Open the input file and verify that all required variables are present.
Definition: LauAbsFitModel.cc:557
LauFitDataTree::getData
const LauFitData & getData(UInt_t iEvt) const
Retrieve the data for a given event.
Definition: LauFitDataTree.cc:376
LauKinematics::getThetaPrime
Double_t getThetaPrime() const
Get theta' value.
Definition: LauKinematics.hh:287
LauFitObject::firstExpt
UInt_t firstExpt() const
Obtain the number of the first experiment.
Definition: LauFitObject.hh:115
LauIsobarDynamics::recalculateNormalisation
void recalculateNormalisation()
recalculate Normalization
Definition: LauIsobarDynamics.cc:232
LauAbsBkgndDPModel
The abstract interface for a background Dalitz plot model.
Definition: LauAbsBkgndDPModel.hh:45
LauKinematics::getmPrime
Double_t getmPrime() const
Get m' value.
Definition: LauKinematics.hh:281
LauDaughters::getNameParent
TString getNameParent() const
Get name of the parent particle.
Definition: LauDaughters.cc:319
LauCPFitModel::writeOutTable
virtual void writeOutTable(const TString &outputFile)
Write the fit results in latex table format.
Definition: LauCPFitModel.cc:1357
LauCPFitModel::posMeanEff_
LauParameter posMeanEff_
The mean efficiency for B+ model.
Definition: LauCPFitModel.hh:617
LauAbsCoeffSet::name
virtual TString name() const
Retrieve the name of the coefficient set.
Definition: LauAbsCoeffSet.hh:134
LauCPFitModel::fakeJacobians_
std::vector< Double_t > fakeJacobians_
The cached values of the sqDP jacobians for each true bin.
Definition: LauCPFitModel.hh:677
LauCPFitModel::negFitFracEffUnCorr_
LauParArray negFitFracEffUnCorr_
Fit B- fractions (uncorrected for the efficiency)
Definition: LauCPFitModel.hh:596
LauCPFitModel::LauCPFitModel
LauCPFitModel(LauIsobarDynamics *negModel, LauIsobarDynamics *posModel, Bool_t tagged=kTRUE, const TString &tagVarName="charge")
Constructor.
Definition: LauCPFitModel.cc:68
LauCPFitModel::setBkgndPdfs
void setBkgndPdfs(const TString &bkgndClass, LauAbsPdf *negPdf, LauAbsPdf *posPdf)
Set the background PDFs.
Definition: LauCPFitModel.cc:433
LauKinematics::withinDPLimits2
Bool_t withinDPLimits2(const Double_t m13Sq, const Double_t m23Sq) const
Check whether a given (m13Sq,m23Sq) point is within the kinematic limits of the Dalitz plot (alternat...
Definition: LauKinematics.cc:477
LauCPFitModel::nNormPar_
UInt_t nNormPar_
Number of normalisation parameters (yields, asymmetries)
Definition: LauCPFitModel.hh:584
LauCPFitModel::useSCFHist_
Bool_t useSCFHist_
Is the SCF fraction DP-dependent.
Definition: LauCPFitModel.hh:656
LauFitDataTree::appendFakePoints
void appendFakePoints(const std::vector< Double_t > &xCoords, const std::vector< Double_t > &yCoords)
Add fake events to the data.
Definition: LauFitDataTree.cc:234
LauAbsFitModel::storeDPEff
Bool_t storeDPEff() const
Determine whether the efficiency information should be stored in the sPlot ntuple.
Definition: LauAbsFitModel.hh:177
LauScfMap.hh
File containing declaration of LauScfMap class.
LauCPFitModel::setBkgndDPModels
void setBkgndDPModels(const TString &bkgndClass, LauAbsBkgndDPModel *negModel, LauAbsBkgndDPModel *posModel)
Set the background DP models.
Definition: LauCPFitModel.cc:359
LauDaughters::squareDP
Bool_t squareDP() const
Determine to use or not the square Dalitz plot.
Definition: LauDaughters.hh:102
LauEffModel::calcEfficiency
Double_t calcEfficiency(const LauKinematics *kinematics) const
Determine the efficiency for a given point in the Dalitz plot.
Definition: LauEffModel.cc:367
LauKinematics::getm13Max
Double_t getm13Max() const
Get the m13 maximum defined as (mParent - m2)
Definition: LauKinematics.hh:371
LauCPFitModel::negSigModel_
LauIsobarDynamics * negSigModel_
The B- signal Dalitz plot model.
Definition: LauCPFitModel.hh:536
LauAbsFitModel::cacheInfo
void cacheInfo(LauPdfPList &pdfList, const LauFitDataTree &theData)
Have all PDFs in the list cache the data.
Definition: LauAbsFitModel.cc:1125
LauCPFitModel::generateExtraPdfValues
void generateExtraPdfValues(LauPdfPList *extraPdfs, LauEmbeddedData *embeddedData)
Generate from the extra PDFs.
Definition: LauCPFitModel.cc:2075
LauCPFitModel::evtCharges_
std::vector< Int_t > evtCharges_
Vector to store event charges.
Definition: LauCPFitModel.hh:650
LauFitNtuple.hh
File containing declaration of LauFitNtuple class.
LauEmbeddedData::clearUsedList
void clearUsedList()
Clear the list of used events.
Definition: LauEmbeddedData.hh:120
LauCPFitModel::calcExtraFractions
void calcExtraFractions(Bool_t initValues=kFALSE)
Calculate the CP-conserving and CP-violating fit fractions.
Definition: LauCPFitModel.cc:1020
LauAbsFitModel::enableEmbedding
Bool_t enableEmbedding() const
Determine whether embedding of events is enabled in the generation.
Definition: LauAbsFitModel.hh:135
LauAbsFitModel::fitData
const LauFitDataTree * fitData() const
Const access the data store.
Definition: LauAbsFitModel.hh:792
LauIsobarDynamics::fillDataTree
void fillDataTree(const LauFitDataTree &fitDataTree)
Fill the internal data structure that caches the resonance dynamics.
Definition: LauIsobarDynamics.cc:2829
LauCPFitModel::scfExtraLike_
Double_t scfExtraLike_
SCF likelihood from extra PDFs.
Definition: LauCPFitModel.hh:734
LauEffModel::setEffHisto
void setEffHisto(const TH2 *effHisto, Bool_t useInterpolation=kTRUE, Bool_t fluctuateBins=kFALSE, Double_t avEff=-1.0, Double_t absError=-1.0, Bool_t useUpperHalfOnly=kFALSE, Bool_t squareDP=kFALSE)
Set the efficiency variation across the phase space using a predetermined 2D histogram.
Definition: LauEffModel.cc:68
LauCPFitModel::negFitFrac_
LauParArray negFitFrac_
The B- fit fractions.
Definition: LauCPFitModel.hh:590
LauKinematics::updateKinematics
void updateKinematics(const Double_t m13Sq, const Double_t m23Sq)
Update all kinematic quantities based on the DP co-ordinates m13Sq and m23Sq.
Definition: LauKinematics.cc:128
LauAbsFitModel::addGenNtupleDoubleBranch
virtual void addGenNtupleDoubleBranch(const TString &name)
Add a branch to the gen tree for storing a double.
Definition: LauAbsFitModel.cc:377
LauCPFitModel::nSigDPPar_
UInt_t nSigDPPar_
Number of signal DP parameters.
Definition: LauCPFitModel.hh:578
LauCPFitModel::setupGenNtupleBranches
void setupGenNtupleBranches()
Setup the required ntuple branches.
Definition: LauCPFitModel.cc:2002
LauCPFitModel::negCoeffs_
std::vector< LauComplex > negCoeffs_
The complex coefficients for B-.
Definition: LauCPFitModel.hh:689
LauAbsRValue::blind
virtual Bool_t blind() const =0
The blinding state.
LauAbsFitModel::fixParamTreeName_
TString fixParamTreeName_
Imported parameters tree name.
Definition: LauAbsFitModel.hh:801
LauCPFitModel::setupBkgndVectors
virtual void setupBkgndVectors()
Define the length of the background vectors.
Definition: LauCPFitModel.cc:137
LauCPFitModel.hh
File containing declaration of LauCPFitModel class.
LauAbsCoeffSet.hh
File containing declaration of LauAbsCoeffSet class.
LauAbsFitModel::allImportedFreeParams_
std::set< LauParameter * > allImportedFreeParams_
The set of parameters that are imported (either from a file or by value) and not.
Definition: LauAbsFitModel.hh:815
LauAbsFitModel::bkgndClassName
const TString & bkgndClassName(UInt_t classID) const
Get the name of a background class from the number.
Definition: LauAbsFitModel.cc:227
LauCPFitModel::negParent_
TString negParent_
Name of the parent particle.
Definition: LauCPFitModel.hh:683
LauCPFitModel::fixParams
void fixParams(std::vector< LauParameter * > &params)
Set a vector of LauParameters according to the specified method.
Definition: LauCPFitModel.cc:3731
LauAbsRValue::value
virtual Double_t value() const =0
Return the value of the parameter.
LauCPFitModel::eventsToGenerate
std::pair< LauGenInfo, Bool_t > eventsToGenerate()
Determine the number of events to generate for each hypothesis.
Definition: LauCPFitModel.cc:1525
LauCPFitModel::fakeSCFFracs_
std::vector< Double_t > fakeSCFFracs_
The cached values of the SCF fraction for each bin centre.
Definition: LauCPFitModel.hh:671
LauCPFitModel::signalAsym_
LauParameter * signalAsym_
Signal asymmetry.
Definition: LauCPFitModel.hh:629
LauEmbeddedData::haveBranch
Bool_t haveBranch(const TString &name) const
Boolean to determine whether branch exists.
Definition: LauEmbeddedData.hh:88
LauIsobarDynamics::getEvtmPrime
Double_t getEvtmPrime() const
Retrieve the square Dalitz plot coordinate, m', for the current event.
Definition: LauIsobarDynamics.hh:343
LauKinematics::withinDPLimits
Bool_t withinDPLimits(const Double_t m13Sq, const Double_t m23Sq) const
Check whether a given (m13Sq,m23Sq) point is within the kinematic limits of the Dalitz plot.
Definition: LauKinematics.cc:439
LauFitObject::multiDimConstrainedPars
const std::set< TString > & multiDimConstrainedPars() const
Const access to the parameter names used in ND constraints.
Definition: LauFitObject.hh:286
LauAbsFitModel::LauParameterList
std::vector< LauParameter > LauParameterList
List of parameters.
Definition: LauAbsFitModel.hh:345
LauCPFitModel::embedNegSignal
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.
Definition: LauCPFitModel.cc:3189
LauCPFitModel::checkInitFitParams
virtual void checkInitFitParams()
Check the initial fit parameters.
Definition: LauCPFitModel.cc:1480
LauAbsFitModel::bkgndClassID
UInt_t bkgndClassID(const TString &className) const
The number assigned to a background class.
Definition: LauAbsFitModel.cc:206
LauIsobarDynamics::getFloatingParameters
std::vector< LauParameter * > & getFloatingParameters()
Retrieve the floating parameters of the resonance models.
Definition: LauIsobarDynamics.hh:576
LauCPFitModel::CPVFitFrac_
LauParArray CPVFitFrac_
The CP violating fit fraction.
Definition: LauCPFitModel.hh:602
LauCPFitModel::usingBkgnd_
Bool_t usingBkgnd_
Background boolean.
Definition: LauCPFitModel.hh:572
LauAbsFitModel::LauParameterPList
std::vector< LauParameter * > LauParameterPList
List of parameter pointers.
Definition: LauAbsFitModel.hh:339
LauFitDataTree::haveBranch
Bool_t haveBranch(const TString &name) const
Check if the named branch is stored.
Definition: LauFitDataTree.cc:293
LauSPlot::NameSet
std::set< TString > NameSet
Type to store names, e.g. of the discriminating/control variables.
Definition: LauSPlot.hh:72
LauCPFitModel::printFitFractions
virtual void printFitFractions(std::ostream &output)
Print the fit fractions, total DP rate and mean efficiency.
Definition: LauCPFitModel.cc:1317
LauAbsResonance::getSpin
Int_t getSpin() const
Get the spin of the resonance.
Definition: LauAbsResonance.hh:176
LauCPFitModel::tagVarName_
const TString tagVarName_
Event charge.
Definition: LauCPFitModel.hh:644
LauFitNtuple::storeCorrMatrix
void storeCorrMatrix(const UInt_t iExpt, const LauAbsFitter::FitStatus &fitStatus, const TMatrixD &covMatrix)
Store the correlation matrix and other fit information.
Definition: LauFitNtuple.cc:76
LauFitDataTree
Class to store the input fit variables.
Definition: LauFitDataTree.hh:50
LauFitObject::fitStatus
const LauAbsFitter::FitStatus & fitStatus() const
Access the fit status information.
Definition: LauFitObject.hh:333
LauAbsFitModel::fixParamFileName_
TString fixParamFileName_
Imported parameters file name.
Definition: LauAbsFitModel.hh:798
LauAbsCoeffSet::baseName
virtual const TString & baseName() const
Retrieve the base name of the coefficient set.
Definition: LauAbsCoeffSet.hh:151
LauKinematics::getm23Max
Double_t getm23Max() const
Get the m23 maximum defined as (mParent - m1)
Definition: LauKinematics.hh:365
LauIsobarDynamics::generate
Bool_t generate()
Generate a toy MC signal event.
Definition: LauIsobarDynamics.cc:2566
LauCPFitModel::LauBkgndPdfsList
std::vector< LauPdfPList > LauBkgndPdfsList
Typedef for a vector of background PDFs.
Definition: LauCPFitModel.hh:274
LauComplex
Class for defining a complex number.
Definition: LauComplex.hh:61
LauCPFitModel::negSignalPdfs_
LauPdfPList negSignalPdfs_
The B- signal PDFs.
Definition: LauCPFitModel.hh:554
LauKinematics::getc12
Double_t getc12() const
Get the cosine of the helicity angle theta12.
Definition: LauKinematics.hh:259
LauIsobarDynamics::getnTotAmp
UInt_t getnTotAmp() const
Retrieve the total number of amplitude components.
Definition: LauIsobarDynamics.hh:516
LauCPFitModel::setNBkgndEvents
virtual void setNBkgndEvents(LauAbsRValue *nBkgndEvents)
Set the background event yield(s)
Definition: LauCPFitModel.cc:220
LauEffModel
Class that implements the efficiency description across the signal Dalitz plot.
Definition: LauEffModel.hh:50
LauCPFitModel::negBkgndPdfs_
LauBkgndPdfsList negBkgndPdfs_
The B- background PDFs.
Definition: LauCPFitModel.hh:566
LauCPFitModel::usePosReweighting_
Bool_t usePosReweighting_
Boolean to use reweighting for B+.
Definition: LauCPFitModel.hh:715
LauCPFitModel::recoSCFFracs_
std::vector< Double_t > recoSCFFracs_
The cached values of the SCF fraction for each event.
Definition: LauCPFitModel.hh:668
LauAsymmCalc.hh
File containing declaration of LauAsymmCalc class.
LauAbsEffModel::passVeto
virtual Bool_t passVeto(const LauKinematics *kinematics) const =0
Determine whether the given DP position is outside the vetoes.
LauCPFitModel::bkgndExtraLike_
std::vector< Double_t > bkgndExtraLike_
Background likelihood value(s) from extra PDFs.
Definition: LauCPFitModel.hh:737
LauCPFitModel::getParamFromTree
Double_t getParamFromTree(TTree &tree, const TString &name)
Retrieve a named parameter from a TTree.
Definition: LauCPFitModel.cc:3690
LauParameter::blind
Bool_t blind() const
The blinding state.
Definition: LauParameter.hh:160
LauAbsFitModel::useRandomInitFitPars
Bool_t useRandomInitFitPars() const
Determine whether the initial values of the fit parameters, in particular the isobar coefficient para...
Definition: LauAbsFitModel.hh:180
LauAbsBkgndDPModel::generate
virtual Bool_t generate()=0
Generate a toy MC event from the model.
LauCPFitModel::posFitFrac_
LauParArray posFitFrac_
The B+ fit fractions.
Definition: LauCPFitModel.hh:593
LauAbsFitModel::addFitParameters
UInt_t addFitParameters(LauParameter *param, const Bool_t addFixed=kFALSE)
Add the given parameter to the list of all fit parameters.
Definition: LauAbsFitModel.cc:918
LauCPFitModel::posDPRate_
LauParameter posDPRate_
The average DP rate for B+.
Definition: LauCPFitModel.hh:623
LauCPFitModel::useNegReweighting_
Bool_t useNegReweighting_
Boolean to use reweighting for B-.
Definition: LauCPFitModel.hh:712
LauCPFitModel::negSignalTree_
LauEmbeddedData * negSignalTree_
The B- signal event tree.
Definition: LauCPFitModel.hh:697
LauCPFitModel::setSignalDPParameters
void setSignalDPParameters()
Set the fit parameters for the DP model.
Definition: LauCPFitModel.cc:748
LauCPFitModel::reuseSignal_
Bool_t reuseSignal_
Boolean to reuse signal events.
Definition: LauCPFitModel.hh:709
LauSPlot::NumbMap
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:75
LauKinematics::squareDP
void squareDP(const Bool_t calcSquareDPCoords)
Enable/disable the calculation of square Dalitz plot co-ordinates.
Definition: LauKinematics.hh:73
LauCPFitModel::~LauCPFitModel
virtual ~LauCPFitModel()
Destructor.
Definition: LauCPFitModel.cc:122
LauIsobarDynamics::getConjResName
TString getConjResName(const TString &resName) const
Retrieve the name of the charge conjugate of a named resonance.
Definition: LauIsobarDynamics.cc:2974
LauKinematics::getm23
Double_t getm23() const
Get the m23 invariant mass.
Definition: LauKinematics.hh:227
LauAbsFitModel::setSPlotNtupleDoubleBranchValue
virtual void setSPlotNtupleDoubleBranchValue(const TString &name, Double_t value)
Set the value of a double branch in the sPlot tree.
Definition: LauAbsFitModel.cc:422
LauFitNtuple
Class to store the results from the fit into an ntuple.
Definition: LauFitNtuple.hh:58
LauComplex.hh
File containing declaration of LauComplex class.
LauFitObject::iExpt
UInt_t iExpt() const
Obtain the number of the current experiment.
Definition: LauFitObject.hh:118
LauAbsFitModel::fillSPlotNtupleBranches
virtual void fillSPlotNtupleBranches()
Fill the sPlot tuple.
Definition: LauAbsFitModel.cc:427
LauCPFitModel::appendBinCentres
void appendBinCentres(LauFitDataTree *inputData)
Append fake data points to the inputData for each bin in the SCF smearing matrix.
Definition: LauCPFitModel.cc:2285
LauDaughters.hh
File containing declaration of LauDaughters class.
LauAbsFitModel::prodPdfValue
Double_t prodPdfValue(LauPdfPList &pdfList, UInt_t iEvt)
Calculate the product of the per-event likelihoods of the PDFs in the list.
Definition: LauAbsFitModel.cc:1132
LauIsobarDynamics::usingScfModel
Bool_t usingScfModel() const
Check whether a self cross feed fraction model is being used.
Definition: LauIsobarDynamics.hh:564
LauIsobarDynamics::getEvtIntensity
Double_t getEvtIntensity() const
Retrieve the total intensity multiplied by the efficiency for the current event.
Definition: LauIsobarDynamics.hh:373
LauIsobarDynamics::checkToyMC
ToyMCStatus checkToyMC(Bool_t printErrorMessages=kTRUE, Bool_t printInfoMessages=kFALSE)
Check the status of the toy MC generation.
Definition: LauIsobarDynamics.cc:2637
LauCPFitModel::setAmpCoeffSet
virtual void setAmpCoeffSet(LauAbsCoeffSet *coeffSet)
Set the DP amplitude coefficients.
Definition: LauCPFitModel.cc:470
LauCPFitModel::recalculateNormalisation
virtual void recalculateNormalisation()
Recalculate Normalization the signal DP models.
Definition: LauCPFitModel.cc:718
LauCPFitModel::acp_
std::vector< LauParameter > acp_
A_CP parameter.
Definition: LauCPFitModel.hh:611
LauCPFitModel::negScfPdfs_
LauPdfPList negScfPdfs_
The B- SCF PDFs.
Definition: LauCPFitModel.hh:560
LauGenNtuple.hh
File containing declaration of LauGenNtuple class.
LauCPFitModel::recoJacobians_
std::vector< Double_t > recoJacobians_
The cached values of the sqDP jacobians for each event.
Definition: LauCPFitModel.hh:674
LauEmbeddedData::nUsedEvents
UInt_t nUsedEvents() const
Retrieve the number of events that have already been sampled.
Definition: LauEmbeddedData.hh:82
LauCPFitModel::embedNegBkgnd
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...
Definition: LauCPFitModel.cc:3217
LauAbsFitModel::addSPlotNtupleIntegerBranch
virtual void addSPlotNtupleIntegerBranch(const TString &name)
Add a branch to the sPlot tree for storing an integer.
Definition: LauAbsFitModel.cc:407
LauCPFitModel::nExtraPdfPar_
UInt_t nExtraPdfPar_
Number of extra PDF parameters.
Definition: LauCPFitModel.hh:581
LauIsobarDynamics::getnCohAmp
UInt_t getnCohAmp() const
Retrieve the number of coherent amplitude components.
Definition: LauIsobarDynamics.hh:522
LauRandom::randomFun
TRandom * randomFun()
Access the singleton random number generator with a particular seed.
Definition: LauRandom.cc:33
LauIsobarDynamics::getExtraParameters
const std::vector< LauParameter > & getExtraParameters() const
Retrieve any extra parameters/quantities (e.g. K-matrix total fit fractions)
Definition: LauIsobarDynamics.hh:570
LauParArray
std::vector< std::vector< LauParameter > > LauParArray
Type to define an array of parameters.
Definition: LauParameter.hh:570
LauCPFitModel::posBkgndDPModels_
LauBkgndDPModelList posBkgndDPModels_
The B+ background Dalitz plot models.
Definition: LauCPFitModel.hh:545
LauAbsPdf
Class for defining the abstract interface for PDF classes.
Definition: LauAbsPdf.hh:54
LauCPFitModel::weightEvents
virtual void weightEvents(const TString &dataFileName, const TString &dataTreeName)
Weight events based on the DP model.
Definition: LauCPFitModel.cc:3319
LauAbsFitModel::setGenNtupleDoubleBranchValue
virtual void setGenNtupleDoubleBranchValue(const TString &name, Double_t value)
Set the value of a double branch in the gen tree.
Definition: LauAbsFitModel.cc:387
LauComplex::abs2
Double_t abs2() const
Obtain the square of the absolute value of the complex number.
Definition: LauComplex.hh:260
LauSimFitTask::fitNtuple
const LauFitNtuple * fitNtuple() const
Const access to the fit ntuple.
Definition: LauSimFitTask.hh:100
LauCPFitModel::negDPRate_
LauParameter negDPRate_
The average DP rate for B-.
Definition: LauCPFitModel.hh:620
LauKinematics::getm12
Double_t getm12() const
Get the m12 invariant mass.
Definition: LauKinematics.hh:221
LauScfMap::binNumber
Int_t binNumber(Double_t xCoord, Double_t yCoord) const
Find the global bin number for the given co-ordinates.
Definition: LauScfMap.cc:158
LauParameter::fixed
Bool_t fixed() const
Check whether the parameter is fixed or floated.
Definition: LauParameter.hh:235
LauParameter::LauParameter
LauParameter()
Default constructor.
Definition: LauParameter.cc:40
LauIsobarDynamics::getDaughters
const LauDaughters * getDaughters() const
Retrieve the daughters.
Definition: LauIsobarDynamics.hh:540
LauCPFitModel::posFitFracEffUnCorr_
LauParArray posFitFracEffUnCorr_
Fit B+ fractions (uncorrected for the efficiency)
Definition: LauCPFitModel.hh:599
LauCPFitModel::twodimPDFs
virtual LauSPlot::TwoDMap twodimPDFs() const
Returns the species and variables for all 2D PDFs in the fit.
Definition: LauCPFitModel.cc:2871
LauCPFitModel::scfMap_
LauScfMap * scfMap_
The smearing matrix for the SCF DP PDF.
Definition: LauCPFitModel.hh:665
LauCPFitModel::embedPosBkgnd
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...
Definition: LauCPFitModel.cc:3282
LauAbsFitModel::writeLatexTable
Bool_t writeLatexTable() const
Determine whether writing out of the latex table is enabled.
Definition: LauAbsFitModel.hh:144
LauKinematics::getm13
Double_t getm13() const
Get the m13 invariant mass.
Definition: LauKinematics.hh:233
LauIsobarDynamics::getEvtm23Sq
Double_t getEvtm23Sq() const
Retrieve the invariant mass squared of the second and third daughters in the current event.
Definition: LauIsobarDynamics.hh:337
LauAbsRValue
Pure abstract base class for defining a parameter containing an R value.
Definition: LauAbsRValue.hh:45
LauEmbeddedData::getReweightedEvent
Bool_t getReweightedEvent(LauIsobarDynamics *dynamics)
Retrieve an event from the data sample, applying an accept/reject based on the given DP model.
Definition: LauEmbeddedData.cc:73
LauCPFitModel::posSignalTree_
LauEmbeddedData * posSignalTree_
The B+ signal event tree.
Definition: LauCPFitModel.hh:700
LauAbsFitModel::extraPars
const LauParameterList & extraPars() const
Const access the extra variables.
Definition: LauAbsFitModel.hh:771
LauCPFitModel::calcAsymmetries
void calcAsymmetries(Bool_t initValues=kFALSE)
Calculate the CP asymmetries.
Definition: LauCPFitModel.cc:1073
LauConstants.hh
File containing LauConstants namespace.
LauCPFitModel::posScfPdfs_
LauPdfPList posScfPdfs_
The B+ SCF PDFs.
Definition: LauCPFitModel.hh:563
LauAbsFitModel::fixParamMap_
std::map< TString, Double_t > fixParamMap_
Map from imported parameter name to value.
Definition: LauAbsFitModel.hh:804
LauCPFitModel::initialise
virtual void initialise()
Initialise the fit.
Definition: LauCPFitModel.cc:539
LauCPFitModel::sigTotalLike_
Double_t sigTotalLike_
Total signal likelihood.
Definition: LauCPFitModel.hh:740
LauCPFitModel::initialiseDPModels
virtual void initialiseDPModels()
Initialise the signal DP models.
Definition: LauCPFitModel.cc:727
LauParameter::name
const TString & name() const
The parameter name.
Definition: LauParameter.hh:154
LauScfMap
Class for representing the 4D smearing matrix for mis-reconstructed signal (self cross feed)
Definition: LauScfMap.hh:45
LauAbsBkgndDPModel.hh
File containing declaration of LauAbsBkgndDPModel class.
LauAbsFitModel::addSPlotNtupleDoubleBranch
virtual void addSPlotNtupleDoubleBranch(const TString &name)
Add a branch to the sPlot tree for storing a double.
Definition: LauAbsFitModel.cc:412
LauCPFitModel::freeSpeciesNames
virtual LauSPlot::NumbMap freeSpeciesNames() const
Returns the names and yields of species that are free in the fit.
Definition: LauCPFitModel.cc:2828
LauIsobarDynamics::initialise
void initialise(const std::vector< LauComplex > &coeffs)
Initialise the Dalitz plot dynamics.
Definition: LauIsobarDynamics.cc:361
LauGenNtuple::addIntegerBranch
void addIntegerBranch(const TString &name)
Add integer branch to tree.
Definition: LauGenNtuple.cc:88
LauCPFitModel::fitFracAsymm_
std::vector< LauParameter > fitFracAsymm_
The fit fraction asymmetries.
Definition: LauCPFitModel.hh:608
LauAbsCoeffSet::acp
virtual LauParameter acp()=0
Calculate the CP asymmetry.
LauAbsFitModel::nBkgndClasses
UInt_t nBkgndClasses() const
Returns the number of background classes.
Definition: LauAbsFitModel.hh:192
LauScfMap::trueBins
const std::vector< Int_t > * trueBins(Int_t recoBin) const
Find which true bins contribute to the given reco bin.
Definition: LauScfMap.cc:169
LauEmbeddedData::nEvents
UInt_t nEvents() const
Retrieve the number of events.
Definition: LauEmbeddedData.hh:76
LauCPFitModel::fixParam
void fixParam(LauParameter *param, const Double_t val, const Bool_t fix)
Set a LauParameter to a given value.
Definition: LauCPFitModel.cc:3712
LauCPFitModel::posKinematics_
LauKinematics * posKinematics_
The B+ Dalitz plot kinematics object.
Definition: LauCPFitModel.hh:551
LauCPFitModel::savePDFPlotsWave
virtual void savePDFPlotsWave(const TString &label, const Int_t &spin)
Save the pdf Plots for the sum of resonances of a given spin.
Definition: LauCPFitModel.cc:3568
LauIsobarDynamics::collateResonanceParameters
void collateResonanceParameters()
Collate the resonance parameters to initialise (or re-initialise) the model.
Definition: LauIsobarDynamics.cc:273
LauCPFitModel::savePDFPlots
virtual void savePDFPlots(const TString &label)
Save the pdf Plots for all the resonances.
Definition: LauCPFitModel.cc:3437
LauAbsFitModel::fitPars
const LauParameterPList & fitPars() const
Const access to the fit variables.
Definition: LauAbsFitModel.hh:765
LauCPFitModel::randomiseInitFitPars
void randomiseInitFitPars()
Randomise the initial fit parameters.
Definition: LauCPFitModel.cc:1493
LauAbsFitModel::clearFitParVectors
void clearFitParVectors()
Clear the vectors containing fit parameters.
Definition: LauAbsFitModel.cc:240
LauCPFitModel::setSignalPdfs
void setSignalPdfs(LauAbsPdf *negPdf, LauAbsPdf *posPdf)
Set the signal PDFs.
Definition: LauCPFitModel.cc:385
LauIsobarDynamics::getEvtDPAmp
const LauComplex & getEvtDPAmp() const
Retrieve the total amplitude for the current event.
Definition: LauIsobarDynamics.hh:325
LauCPFitModel::generateSignalEvent
Bool_t generateSignalEvent()
Generate signal event.
Definition: LauCPFitModel.cc:1798
LauIsobarDynamics::getEvtScfFraction
Double_t getEvtScfFraction() const
Retrieve the fraction of events that are poorly reconstructed (the self cross feed fraction) for the ...
Definition: LauIsobarDynamics.hh:361
LauCPFitModel::addSPlotNtupleBranches
void addSPlotNtupleBranches(const LauPdfPList *extraPdfs, const TString &prefix)
Add sPlot branches for the extra PDFs.
Definition: LauCPFitModel.cc:2685
LauCPFitModel::tagged_
const Bool_t tagged_
IS the analysis tagged?
Definition: LauCPFitModel.hh:641
LauFitNtuple::storeParsAndErrors
void storeParsAndErrors(const std::vector< LauParameter * > &fitVars, const std::set< TString > &constrainedVars, const std::vector< LauParameter > &extraVars)
Store parameters and their errors.
Definition: LauFitNtuple.cc:216
LauCPFitModel::propagateParUpdates
virtual void propagateParUpdates()
Calculate things that depend on the fit parameters after they have been updated by Minuit.
Definition: LauCPFitModel.cc:2141
LauFitData
std::map< TString, Double_t > LauFitData
Type for holding event data.
Definition: LauFitDataTree.hh:45
LauAbsCoeffSet::index
virtual UInt_t index() const
Retrieve the index number of the coefficient set.
Definition: LauAbsCoeffSet.hh:166
LauIsobarDynamics::modifyDataTree
void modifyDataTree()
Recache the amplitude values for those that have changed.
Definition: LauIsobarDynamics.cc:2792
LauCPFitModel::getDPAmps
virtual std::map< TString, LauComplex > getDPAmps(const Double_t m13Sq, const Double_t m23Sq)
Calculate the DP amplitude(s) for a given DP position.
Definition: LauCPFitModel.cc:3110
LauKinematics::updateSqDPKinematics
void updateSqDPKinematics(const Double_t mPrime, const Double_t thetaPrime)
Update all kinematic quantities based on the square DP co-ordinates m' and theta'.
Definition: LauKinematics.cc:147
LauCPFitModel::storeSignalMCMatch
Bool_t storeSignalMCMatch(LauEmbeddedData *embeddedData)
Store the MC truth info on the TM/SCF nature of the embedded signal event.
Definition: LauCPFitModel.cc:2115
LauCPFitModel::coeffPars_
std::vector< LauAbsCoeffSet * > coeffPars_
Magnitudes and Phases.
Definition: LauCPFitModel.hh:587
LauCPFitModel::scfDPLike_
Double_t scfDPLike_
SCF DP likelihood value.
Definition: LauCPFitModel.hh:725
LauCPFitModel::setSPlotNtupleBranchValues
Double_t setSPlotNtupleBranchValues(LauPdfPList *extraPdfs, const TString &prefix, UInt_t iEvt)
Set the branches for the sPlot ntuple with extra PDFs.
Definition: LauCPFitModel.cc:2739
LauAbsFitModel::addGenNtupleIntegerBranch
virtual void addGenNtupleIntegerBranch(const TString &name)
Add a branch to the gen tree for storing an integer.
Definition: LauAbsFitModel.cc:372
LauAbsRValue::fixed
virtual Bool_t fixed() const =0
Check is the parameter is fixed or floated.
LauParameter::range
Double_t range() const
The range allowed for the parameter.
Definition: LauParameter.hh:229
LauParameter::updatePull
void updatePull()
Call to update the bias and pull values.
Definition: LauParameter.cc:611
LauAbsResonance
Abstract class for defining type for resonance amplitude models (Breit-Wigner, Flatte etc....
Definition: LauAbsResonance.hh:51
LauPrint.hh
File containing declaration of LauPrint class and LauOutputLevel enum.
LauCPFitModel::getEventSum
virtual Double_t getEventSum() const
Get the total number of events.
Definition: LauCPFitModel.cc:2397
LauCPFitModel::getEvtSCFDPLikelihood
virtual Double_t getEvtSCFDPLikelihood(UInt_t iEvt)
Calculate the SCF likelihood for the DP for a given event.
Definition: LauCPFitModel.cc:2495
LauCPFitModel::bkgndEvents_
LauBkgndYieldList bkgndEvents_
Background yield(s)
Definition: LauCPFitModel.hh:635
LauPrint
Class to define various output print commands.
Definition: LauPrint.hh:54
LauEmbeddedData::getEmbeddedEvent
void getEmbeddedEvent(LauKinematics *kinematics)
Retrieve an event from the data sample.
Definition: LauEmbeddedData.cc:155
LauCPFitModel::signalEvents_
LauParameter * signalEvents_
Signal yield.
Definition: LauCPFitModel.hh:626
LauAsymmCalc::getAsymmetry
Double_t getAsymmetry()
Retrieve the asymmetry.
Definition: LauAsymmCalc.hh:61
LauCPFitModel::forceAsym_
Bool_t forceAsym_
Option to force an asymmetry.
Definition: LauCPFitModel.hh:632
LauCPFitModel::negBkgndDPModels_
LauBkgndDPModelList negBkgndDPModels_
The B- background Dalitz plot models.
Definition: LauCPFitModel.hh:542
LauAbsFitModel::useDP
Bool_t useDP() const
Is the Dalitz plot term in the likelihood.
Definition: LauAbsFitModel.hh:98
LauCPFitModel::LauGenInfo
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.
Definition: LauCPFitModel.hh:268
LauCPFitModel::setNSigEvents
virtual void setNSigEvents(LauParameter *nSigEvents)
Set the signal event yield.
Definition: LauCPFitModel.cc:154
LauFitObject::covarianceMatrix
const TMatrixD & covarianceMatrix() const
Access the fit covariance matrix.
Definition: LauFitObject.hh:345
LauEmbeddedData::getValues
LauFitData getValues(const std::vector< TString > &names) const
Get values of specified branches.
Definition: LauEmbeddedData.cc:205
LauKinematics
Class for calculating 3-body kinematic quantities.
Definition: LauKinematics.hh:45
LauCPFitModel::generateBkgndEvent
Bool_t generateBkgndEvent(UInt_t bgID)
Generate background event.
Definition: LauCPFitModel.cc:1940
LauIsobarDynamics::getDPRate
const LauParameter & getDPRate() const
Retrieve the overall Dalitz plot rate.
Definition: LauIsobarDynamics.hh:495
LauAbsFitModel::setSPlotNtupleIntegerBranchValue
virtual void setSPlotNtupleIntegerBranchValue(const TString &name, Int_t value)
Set the value of an integer branch in the sPlot tree.
Definition: LauAbsFitModel.cc:417
LauCPFitModel::setExtraPdfParameters
void setExtraPdfParameters()
Set the fit parameters for the extra PDFs.
Definition: LauCPFitModel.cc:774
LauAbsFitModel::clearExtraVarVectors
void clearExtraVarVectors()
Clear the vectors containing extra ntuple variables.
Definition: LauAbsFitModel.cc:250
LauCPFitModel::fixdSpeciesNames
virtual LauSPlot::NumbMap fixdSpeciesNames() const
Returns the names and yields of species that are fixed in the fit.
Definition: LauCPFitModel.cc:2852
LauCPFitModel::updateCoeffs
virtual void updateCoeffs()
Update the coefficients.
Definition: LauCPFitModel.cc:2606
LauCPFitModel::setSCFPdfs
void setSCFPdfs(LauAbsPdf *negPdf, LauAbsPdf *posPdf)
Set the SCF PDF for a given variable.
Definition: LauCPFitModel.cc:409
LauIsobarDynamics::getMeanEff
const LauParameter & getMeanEff() const
Retrieve the mean efficiency across the Dalitz plot.
Definition: LauIsobarDynamics.hh:489
LauCPFitModel::negMeanEff_
LauParameter negMeanEff_
The mean efficiency for B- model.
Definition: LauCPFitModel.hh:614
LauCPFitModel::sigExtraLike_
Double_t sigExtraLike_
Signal likelihood from extra PDFs.
Definition: LauCPFitModel.hh:731
LauKinematics::calcSqDPJacobian
Double_t calcSqDPJacobian(const Double_t mPrime, const Double_t thPrime) const
Calculate the Jacobian for the transformation m23^2, m13^2 -> m', theta' (square DP) at the given poi...
Definition: LauKinematics.cc:176
LauAbsFitModel::fixParams_
Bool_t fixParams_
Whether to fix the loaded parameters (kTRUE) or leave them floating (kFALSE)
Definition: LauAbsFitModel.hh:810
LauParameter::genValue
Double_t genValue() const
The value generated for the parameter.
Definition: LauParameter.hh:205
LauFitObject::nExpt
UInt_t nExpt() const
Obtain the number of experiments.
Definition: LauFitObject.hh:112
LauAsymmCalc
Class for calculating the asymmetry between two variables.
Definition: LauAsymmCalc.hh:38
LauAbsFitModel::printFitParameters
void printFitParameters(const LauPdfPList &pdfList, std::ostream &fout) const
Print the fit parameters for all PDFs in the list.
Definition: LauAbsFitModel.cc:1097
LauAbsEffModel
Pure abstract base class for defining the efficiency description across the signal Dalitz plot.
Definition: LauAbsEffModel.hh:45
LauFitObject::eventsPerExpt
UInt_t eventsPerExpt() const
Obtain the total number of events in the current experiment.
Definition: LauFitObject.hh:109
LauPrint::printFormat
void printFormat(std::ostream &stream, Double_t value) const
Method to choose the printing format to a specified level of precision.
Definition: LauPrint.cc:43
LauCPFitModel::posCoeffs_
std::vector< LauComplex > posCoeffs_
The complex coefficients for B+.
Definition: LauCPFitModel.hh:692
LauCPFitModel::updateSigEvents
void updateSigEvents()
Update the signal events after Minuit sets background parameters.
Definition: LauCPFitModel.cc:2156
LauIsobarDynamics::getEffModel
const LauAbsEffModel * getEffModel() const
Retrieve the model for the efficiency across the Dalitz plot.
Definition: LauIsobarDynamics.hh:558
LauEmbeddedData::findBranches
Bool_t findBranches()
Find and read the branches in data tree.
Definition: LauEmbeddedData.cc:54
LauIsobarDynamics::updateCoeffs
void updateCoeffs(const std::vector< LauComplex > &coeffs)
Update the complex coefficients for the resonances.
Definition: LauIsobarDynamics.cc:2952
LauSPlot::TwoDMap
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:81
LauCPFitModel::bkgndTotalLike_
std::vector< Double_t > bkgndTotalLike_
Total background likelihood(s)
Definition: LauCPFitModel.hh:746
LauAbsRValue::genValue
virtual Double_t genValue() const =0
The value generated for the parameter.
LauEmbeddedData
Class to store the data for embedding in toy experiments.
Definition: LauEmbeddedData.hh:46
LauDaughters
Class that defines the particular 3-body decay under study.
Definition: LauDaughters.hh:47
LauCPFitModel::posBkgndTree_
LauBkgndEmbDataList posBkgndTree_
The B+ background event tree.
Definition: LauCPFitModel.hh:706
LauFitNtuple::updateFitNtuple
void updateFitNtuple()
Update the fit ntuple.
Definition: LauFitNtuple.cc:254
LauKinematics::getm23Sq
Double_t getm23Sq() const
Get the m23 invariant mass square.
Definition: LauKinematics.hh:245
LauIsobarDynamics::getEvtthPrime
Double_t getEvtthPrime() const
Retrieve the square Dalitz plot coordinate, theta', for the current event.
Definition: LauIsobarDynamics.hh:349
LauCPFitModel::scfFrac_
LauParameter scfFrac_
The (global) SCF fraction parameter.
Definition: LauCPFitModel.hh:659
LauIsobarDynamics::getFullAmplitude
LauComplex getFullAmplitude(const Int_t resID) const
Retrieve the Amplitude of resonance resID.
Definition: LauIsobarDynamics.hh:399
LauCPFitModel::curEvtCharge_
Int_t curEvtCharge_
Current event charge.
Definition: LauCPFitModel.hh:647
LauIsobarDynamics
Class for defining signal dynamics using the isobar model.
Definition: LauIsobarDynamics.hh:53
LauCPFitModel::splitSignalComponent
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.
Definition: LauCPFitModel.cc:313
LauAbsFitModel::fillGenNtupleBranches
virtual void fillGenNtupleBranches()
Fill the gen tuple branches.
Definition: LauAbsFitModel.cc:402
LauEmbeddedData.hh
File containing declaration of LauEmbeddedData class.
LauKinematics::getm13Sq
Double_t getm13Sq() const
Get the m13 invariant mass square.
Definition: LauKinematics.hh:251
LauDaughters::gotSymmetricalDP
Bool_t gotSymmetricalDP() const
Is Dalitz plot symmetric, i.e. 2 identical particles.
Definition: LauDaughters.hh:84
LauCPFitModel::setFitNEvents
void setFitNEvents()
Set the initial yields.
Definition: LauCPFitModel.cc:810
LauCPFitModel::posSigModel_
LauIsobarDynamics * posSigModel_
The B+ signal Dalitz plot model.
Definition: LauCPFitModel.hh:539
LauAbsFitModel::updateFitParameters
void updateFitParameters(LauPdfPList &pdfList)
Update the fit parameters for all PDFs in the list.
Definition: LauAbsFitModel.cc:1090
LauCPFitModel::variableNames
virtual LauSPlot::NameSet variableNames() const
Returns the names of all variables in the fit.
Definition: LauCPFitModel.cc:2804
LauAbsRValue::isLValue
virtual Bool_t isLValue() const =0
Is the parameter also an L value or not.
LauCPFitModel::nSigComp_
UInt_t nSigComp_
Number of signal components.
Definition: LauCPFitModel.hh:575
LauCPFitModel::sigDPLike_
Double_t sigDPLike_
Signal DP likelihood value.
Definition: LauCPFitModel.hh:722
LauIsobarDynamics::calcExtraInfo
void calcExtraInfo(const Bool_t init=kFALSE)
Calculate the fit fractions, mean efficiency and total DP rate.
Definition: LauIsobarDynamics.cc:2199
LauCPFitModel::useSCF_
Bool_t useSCF_
Is the signal split into TM and SCF.
Definition: LauCPFitModel.hh:653
LauIsobarDynamics::getKinematics
const LauKinematics * getKinematics() const
Retrieve the Dalitz plot kinematics.
Definition: LauIsobarDynamics.hh:546
LauAbsRValue::name
virtual const TString & name() const =0
Return the name of the parameter.
LauCPFitModel::bkgndDPLike_
std::vector< Double_t > bkgndDPLike_
Background DP likelihood value(s)
Definition: LauCPFitModel.hh:728
LauCPFitModel::reuseBkgnd_
LauBkgndReuseEventsList reuseBkgnd_
Vector of booleans to reuse background events.
Definition: LauCPFitModel.hh:718
LauParameter::initValue
Double_t initValue() const
The initial value of the parameter.
Definition: LauParameter.hh:211
LauCPFitModel::negBkgndTree_
LauBkgndEmbDataList negBkgndTree_
The B- background event tree.
Definition: LauCPFitModel.hh:703
LauFitDataTree::nEvents
UInt_t nEvents() const
Retrieve the number of events.
Definition: LauFitDataTree.hh:119
LauAbsFitModel
Abstract interface to the fitting and toy MC model.
Definition: LauAbsFitModel.hh:88
LauAbsFitModel::fixParamNames_
std::set< TString > fixParamNames_
Imported parameter names.
Definition: LauAbsFitModel.hh:807
LauIsobarDynamics.hh
File containing declaration of LauIsobarDynamics class.
LauAbsFitModel::validBkgndClass
Bool_t validBkgndClass(const TString &className) const
Check if the given background class is in the list.
Definition: LauAbsFitModel.cc:187
LauIsobarDynamics::getEvtm13Sq
Double_t getEvtm13Sq() const
Retrieve the invariant mass squared of the first and third daughters in the current event.
Definition: LauIsobarDynamics.hh:331
LauScfMap::prob
Double_t prob(Int_t recoBin, Int_t trueBin) const
Probability of a true event in the given true bin migrating to the reco bin.
Definition: LauScfMap.cc:179
LauCPFitModel::getDPLikelihoods
virtual std::map< TString, Double_t > getDPLikelihoods(const Double_t m13Sq, const Double_t m23Sq)
Calculate the DP likelihood(s) for a given DP position.
Definition: LauCPFitModel.cc:3134
LauKinematics.hh
File containing declaration of LauKinematics class.
LauGenNtuple
Class to store the results from the toy MC generation into an ntuple.
Definition: LauGenNtuple.hh:45
LauCPFitModel::CPCFitFrac_
LauParArray CPCFitFrac_
The CP conserving fit fraction.
Definition: LauCPFitModel.hh:605
LauCPFitModel::posSignalPdfs_
LauPdfPList posSignalPdfs_
The B+ signal PDFs.
Definition: LauCPFitModel.hh:557
LauAbsFitModel::LauPdfPList
std::vector< LauAbsPdf * > LauPdfPList
List of Pdfs.
Definition: LauAbsFitModel.hh:337
LauCPFitModel::finaliseFitResults
virtual void finaliseFitResults(const TString &tablePrefixName)
Get the fit results and store them.
Definition: LauCPFitModel.cc:1092
LauIsobarDynamics::getFitFractions
const LauParArray & getFitFractions() const
Retrieve the fit fractions for the amplitude components.
Definition: LauIsobarDynamics.hh:501
Generated by   doxygen 1.8.18