laura is hosted by Hepforge, IPPP Durham
Laura++  v1r0
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauIsobarDynamics.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2005 - 2013.
3 // Distributed under the Boost Software License, Version 1.0.
4 // (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5 
6 // Authors:
7 // Thomas Latham
8 // John Back
9 // Paul Harrison
10 
15 #include <iostream>
16 #include <iomanip>
17 #include <fstream>
18 using std::cout;
19 using std::cerr;
20 using std::endl;
21 using std::setprecision;
22 using std::ifstream;
23 using std::ofstream;
24 
25 #include "TFile.h"
26 #include "TRandom.h"
27 #include "TSystem.h"
28 
29 #include "LauAbsResonance.hh"
30 #include "LauBelleNR.hh"
31 #include "LauBelleSymNR.hh"
32 #include "LauConstants.hh"
33 #include "LauDaughters.hh"
34 #include "LauEffModel.hh"
35 #include "LauFitDataTree.hh"
36 #include "LauGounarisSakuraiRes.hh"
37 #include "LauIntegrals.hh"
38 #include "LauIsobarDynamics.hh"
39 #include "LauKinematics.hh"
40 #include "LauKMatrixPropagator.hh"
41 #include "LauKMatrixPropFactory.hh"
42 #include "LauKMatrixProdPole.hh"
43 #include "LauKMatrixProdSVP.hh"
44 #include "LauNRAmplitude.hh"
45 #include "LauPrint.hh"
46 #include "LauRandom.hh"
47 #include "LauRelBreitWignerRes.hh"
48 #include "LauResonanceMaker.hh"
49 
50 ClassImp(LauIsobarDynamics)
51 
52 
53 // for Kpipi: only one scfFraction 2D histogram is needed
54 LauIsobarDynamics::LauIsobarDynamics(LauDaughters* daughters, LauEffModel* effModel, LauEffModel* scfFractionModel) :
55  LauAbsDPDynamics(daughters, effModel, scfFractionModel),
56  symmetricalDP_(kFALSE),
57  integralsDone_(kFALSE),
58  intFileName_("integ.dat"),
59  m13BinWidth_(0.005),
60  m23BinWidth_(0.005),
61  m13Sq_(0.0),
62  m23Sq_(0.0),
63  mPrime_(0.0),
64  thPrime_(0.0),
65  eff_(0.0),
66  scfFraction_(0.0),
67  jacobian_(0.0),
68  ASq_(0.0),
69  evtLike_(0.0),
70  iterationsMax_(100000),
71  nSigGenLoop_(0),
72  aSqMaxSet_(1.25),
73  aSqMaxVar_(0.0),
74  BelleNRAlpha_(0.0),
75  LASSScatteringLength_(0.0),
76  LASSEffectiveRange_(0.0),
77  LASSResonanceMag_(0.0),
78  LASSResonancePhase_(0.0),
79  LASSBackgroundMag_(0.0),
80  LASSBackgroundPhase_(0.0),
81  LASSCutOff_(0.0),
82  changeLASSScatteringLength_(kFALSE),
83  changeLASSEffectiveRange_(kFALSE),
84  changeLASSResonanceMag_(kFALSE),
85  changeLASSResonancePhase_(kFALSE),
86  changeLASSBackgroundMag_(kFALSE),
87  changeLASSBackgroundPhase_(kFALSE),
88  changeLASSCutOff_(kFALSE),
89  FlatteParameterg1_(0.0),
90  FlatteParameterg2_(0.0),
91  changeFlatteParameterg1_(kFALSE),
92  changeFlatteParameterg2_(kFALSE),
93  resBarrierRadius_(4.0),
94  parBarrierRadius_(4.0),
95  barrierType_( LauAbsResonance::BWPrimeBarrier ),
96  flipHelicity_(kTRUE)
97 {
98  // Constructor for the isobar signal model
99  if (daughters != 0) {
100  symmetricalDP_ = daughters->gotSymmetricalDP();
101  typDaug_.push_back(daughters->getTypeDaug1());
102  typDaug_.push_back(daughters->getTypeDaug2());
103  typDaug_.push_back(daughters->getTypeDaug3());
104  }
105  sigResonances_.clear();
106  kMatrixPropagators_.clear();
107  kMatrixPropSet_.clear();
108 }
109 
110 // for Kspipi, we need a scfFraction 2D histogram for each tagging category. They are provided by the map.
111 // Also, we need to know the place that the tagging category of the current event occupies in the data structure inputFitTree
113  LauAbsDPDynamics(daughters, effModel, scfFractionModel),
114  symmetricalDP_(kFALSE),
115  integralsDone_(kFALSE),
116  intFileName_("integ.dat"),
117  m13BinWidth_(0.005),
118  m23BinWidth_(0.005),
119  m13Sq_(0.0),
120  m23Sq_(0.0),
121  mPrime_(0.0),
122  thPrime_(0.0),
123  eff_(0.0),
124  scfFraction_(0.0),
125  jacobian_(0.0),
126  ASq_(0.0),
127  evtLike_(0.0),
128  iterationsMax_(100000),
129  nSigGenLoop_(0),
130  aSqMaxSet_(1.25),
131  aSqMaxVar_(0.0),
132  BelleNRAlpha_(0.0),
133  LASSScatteringLength_(0.0),
134  LASSEffectiveRange_(0.0),
135  LASSResonanceMag_(0.0),
136  LASSResonancePhase_(0.0),
137  LASSBackgroundMag_(0.0),
138  LASSBackgroundPhase_(0.0),
139  LASSCutOff_(0.0),
140  changeLASSScatteringLength_(kFALSE),
141  changeLASSEffectiveRange_(kFALSE),
142  changeLASSResonanceMag_(kFALSE),
143  changeLASSResonancePhase_(kFALSE),
144  changeLASSBackgroundMag_(kFALSE),
145  changeLASSBackgroundPhase_(kFALSE),
146  changeLASSCutOff_(kFALSE),
147  FlatteParameterg1_(0.0),
148  FlatteParameterg2_(0.0),
149  changeFlatteParameterg1_(kFALSE),
150  changeFlatteParameterg2_(kFALSE),
151  resBarrierRadius_(4.0),
152  parBarrierRadius_(4.0),
153  barrierType_( LauAbsResonance::BWPrimeBarrier ),
154  flipHelicity_(kTRUE)
155 {
156  // Constructor for the isobar signal model
157  if (daughters != 0) {
158  symmetricalDP_ = daughters->gotSymmetricalDP();
159  typDaug_.push_back(daughters->getTypeDaug1());
160  typDaug_.push_back(daughters->getTypeDaug2());
161  typDaug_.push_back(daughters->getTypeDaug3());
162  }
163  sigResonances_.clear();
164  kMatrixPropagators_.clear();
165  kMatrixPropSet_.clear();
166 }
167 
169 {
170 }
171 
172 void LauIsobarDynamics::initialise(const std::vector<LauComplex>& coeffs)
173 {
174  // Check whether we have a valid set of integration constants for
175  // the normalisation of the signal likelihood function.
176  this->initialiseVectors();
177 
178  integralsDone_ = kFALSE;
179 
180  // Print summary of what we have so far to screen
181  this->initSummary();
182 
183  if ( nAmp_ == 0 ) {
184  cout << "No contributions to DP model, not performing normalisation integrals." << endl;
185  } else {
186 
187  // We need to calculate the normalisation constants for the
188  // Dalitz plot generation/fitting.
189  cout<<"Starting special run to generate the integrals for normalising the PDF..."<<endl;
190  this->calcDPNormalisation();
191 
192  // Write the integrals to a file (mainly for debugging purposes)
193  this->writeIntegralsFile();
194  }
195 
196 
197  integralsDone_ = kTRUE;
198 
199  cout << setprecision(10);
200 
201  // Calculate and cache the relative normalisations of each resonance _dynamic_ amplitude
202  // (e.g. Breit-Wigner contribution, not from the complex amplitude/phase) w.r.t. the
203  // total DP amplitude. These are stored in fNorm_[i].
204  // The normalisation uses fSqSum[i], which is calculated within the dynamics() function,
205  // which has already been called by the calcDPNomalisation() function above to create the
206  // normalisation integrals.
207  // fSqSum[i] is the event-by-event running total of the dynamical amplitude
208  // squared for a given resonance, i. We require that:
209  // |fNorm_[i]|^2 * |fSqSum[i]|^2 = 1,
210  // i.e. fNorm_[i] normalises each resonance contribution to give the same number of
211  // events in the DP, accounting for the total DP area and the
212  // dynamics of the resonance.
213 
214  for (UInt_t i = 0; i < nAmp_; i++) {
215  fNorm_[i] = 0.0;
216  if (fSqSum_[i] > 0.0) {fNorm_[i] = TMath::Sqrt(1.0/(fSqSum_[i]));}
217  cout<<" fNorm["<<i<<"] = "<<fNorm_[i]<<endl;
218  cout<<" fSqSum["<<i<<"] = "<<fSqSum_[i]<<endl;
219  }
220 
221  for (UInt_t i = 0; i < nAmp_; i++) {
222  for (UInt_t j = 0; j < nAmp_; j++) {
223  cout<<" fifjEffSum["<<i<<"]["<<j<<"] = "<<fifjEffSum_[i][j];
224  }
225  cout<<endl;
226  }
227 
228  for (UInt_t i = 0; i < nAmp_; i++) {
229  for (UInt_t j = 0; j < nAmp_; j++) {
230  cout<<" fifjSum["<<i<<"]["<<j<<"] = "<<fifjSum_[i][j];
231  }
232  cout<<endl;
233  }
234 
235  // Calculate the initial fit fractions (for later comparison with Toy MC, if required)
236  this->updateCoeffs(coeffs);
237  this->calcExtraInfo(kTRUE);
238  for (UInt_t i = 0; i < nAmp_; i++) {
239  for (UInt_t j = i; j < nAmp_; j++) {
240  cout<<"Initial fit fraction for amplitude ("<<i<<","<<j<<") = "<<fitFrac_[i][j].genValue()<<endl;
241  }
242  }
243 
244  cout<<"Initial efficiency = "<<meanDPEff_.initValue()<<endl;
245 
246  cout<<"Initial DPRate = "<<DPRate_.initValue()<<endl;
247 }
248 
250 {
251 
252  UInt_t i(0);
253  TString nameP = daughters_->getNameParent();
254  TString name1 = daughters_->getNameDaug1();
255  TString name2 = daughters_->getNameDaug2();
256  TString name3 = daughters_->getNameDaug3();
257 
258  cout<<"We are going to do a DP with "<<nameP<<" going to "<<name1<<" "<<name2<<" "<<name3<<endl;
259 
260  cout<<"For the following resonance combinations:"<<endl;
261 
262  cout<<"In tracks 2 and 3:"<<endl;
263  for (i = 0; i < nAmp_; i++) {
264  if (resPairAmp_[i] == 1) {
265  cout<<" "<<resTypAmp_[i]<<" to " <<name2<<", "<< name3<<endl;
266  }
267  }
268 
269  cout<<endl;
270 
271  cout<<"In tracks 1 and 3:"<<endl;
272  for (i = 0; i < nAmp_; i++) {
273  if (resPairAmp_[i] == 2) {
274  cout<<" "<<resTypAmp_[i]<<" to "<<name1<<", "<< name3<<endl;
275  }
276  }
277 
278  cout<<endl;
279 
280  cout<<"In tracks 1 and 2:"<<endl;
281  for (i = 0; i < nAmp_; i++) {
282  if (resPairAmp_[i] == 3) {
283  cout<<" "<<resTypAmp_[i]<<" to "<<name1<<", "<< name2<<endl;
284  }
285  }
286 
287  cout<<endl;
288 
289  for (i = 0; i < nAmp_; i++) {
290  if (resPairAmp_[i] == 0) {
291  cout<<"and a non-resonant amplitude."<<endl;
292  }
293  }
294 
295  cout<<endl;
296 }
297 
299 {
300  fSqSum_.clear(); fSqSum_.resize(nAmp_);
301  fifjEffSum_.clear(); fifjEffSum_.resize(nAmp_);
302  fifjSum_.clear(); fifjSum_.resize(nAmp_);
303  ff_.clear(); ff_.resize(nAmp_);
304  fNorm_.clear(); fNorm_.resize(nAmp_);
305  fitFrac_.clear(); fitFrac_.resize(nAmp_);
306 
307  LauComplex null(0.0, 0.0);
308 
309  for (UInt_t i = 0; i < nAmp_; i++) {
310 
311  fSqSum_[i] = 0.0; fNorm_[i] = 0.0;
312  ff_[i] = null;
313  fifjEffSum_[i].resize(nAmp_);
314  fifjSum_[i].resize(nAmp_);
315  fitFrac_[i].resize(nAmp_);
316 
317  for (UInt_t j = 0; j < nAmp_; j++) {
318  fifjEffSum_[i][j] = null;
319  fifjSum_[i][j] = null;
320  fitFrac_[i][j].valueAndRange(0.0, -100.0, 100.0);
321  }
322  }
323 
324  UInt_t nKMatrixProps = kMatrixPropagators_.size();
325  extraParameters_.clear();
326  extraParameters_.resize( nKMatrixProps );
327 
328  for ( UInt_t i = 0; i < nKMatrixProps; ++i ) {
329  extraParameters_[i].valueAndRange(0.0, -100.0, 100.0);
330  }
331 }
332 
333 
335 {
336  // Routine to end integration calculation for loglike normalisation.
337  // This writes out the normalisation integral output into the file named
338  // outputFileName.
339  cout<<"Writing integral output to integrals file "<<intFileName_.Data()<<endl;
340 
341  UInt_t i(0), j(0);
342  ofstream getChar(intFileName_.Data());
343 
344  getChar << setprecision(10);
345 
346  // Write out daughter types (pi, pi0, K, K0s?)
347  for (i = 0; i < 3; i++) {
348  getChar << typDaug_[i] << " ";
349  }
350 
351  // Write out number of resonances in the Dalitz plot model
352  getChar << nAmp_ << endl;
353 
354  // Write out the resonances
355  for (i = 0; i < nAmp_; i++) {
356  getChar << resTypAmp_[i] << " ";
357  }
358 
359  getChar << endl;
360 
361  // Write out the resonance model types (BW, RelBW etc...)
362  for (i = 0; i < nAmp_; i++) {
363  LauAbsResonance* theResonance = sigResonances_[i];
364  Int_t resModelInt = theResonance->getResonanceModel();
365  getChar << resModelInt << " ";
366  }
367 
368  getChar << endl;
369 
370  // Write out the track pairings for each resonance. This is specified
371  // by the resPairAmpInt integer in the addResonance function.
372  for (i = 0; i < nAmp_; i++) {
373  getChar << resPairAmp_[i] << " ";
374  }
375 
376  getChar << endl;
377 
378  // Write out the fSqSum = |ff|^2, where ff = resAmp()
379  for (i = 0; i < nAmp_; i++) {
380  getChar << fSqSum_[i] << " ";
381  }
382 
383  getChar << endl;
384 
385  // Write out the f_i*f_j_conj*eff values = resAmp_i*resAmp_j_conj*eff.
386  // Note that only the top half of the i*j "matrix" is required, as it
387  // is symmetric w.r.t i, j.
388  for (i = 0; i < nAmp_; i++) {
389  for (j = i; j < nAmp_; j++) {
390  getChar << fifjEffSum_[i][j] << " ";
391  }
392  }
393 
394  getChar << endl;
395 
396  // Similar to fifjEffSum, but without the efficiency term included.
397  for (i = 0; i < nAmp_; i++) {
398  for (j = i; j < nAmp_; j++) {
399  getChar << fifjSum_[i][j] << " ";
400  }
401  }
402 
403  getChar << endl;
404 
405 }
406 
407 void LauIsobarDynamics::addResonance(const TString& resName, Int_t resPairAmpInt, const TString& resType, Double_t newMass, Double_t newWidth, Int_t newSpin)
408 {
409  // Function to add a resonance in a Dalitz plot (for signal). To add a resonance
410  // in the continuum background, please use addBgResonance().
411  // No check is made w.r.t flavour
412  // and charge conservation rules, and so the user is responsible for checking the
413  // internal consistency of their function statements with these laws. For example,
414  // the program will not prevent the user from asking for a rho resonance
415  // in a K-pi pair or a K* resonance in a pi-pi pair.
416  // However, to assist the user, a summary of the resonant structure requested
417  // by the user is printed before the program runs. It is important to check this
418  // information when you first define your Dalitz plot model before doing
419  // any fitting/generating.
420  // Arguments are: resonance name, integer to specify the resonance track pairing
421  // (1 => m_23, 2 => m_13, 3 => m_12), i.e. the bachelor track number.
422  // Also required are the initial values of the magnitude and phase of the amplitude
423  // for this resonance, and whether they are to be fixed in any fit.
424  // Supported resonances for the resName (1st) argument are (case-sensitive):
425  // Rho, Kst, f0_980, phi1020, f2_1270, Kst1410, Kst1430
426  // K2*1430, Rho1450, Kst1680, Chic0, NonReson, f2_1525, f0_1500,
427  // Rho1700, f0_1370, (alternative) NRModel, a2_1320, Rho3_1690, a0_980
428  // a0_1450, phi1680, a6_2450, phi3_1850.
429  // The final argument resType specifies whether the resonance is a Breit-Wigner (BW)
430  // Relativistic Breit-Wigner (RelBW) or Flatte distribution (Flatte), for example.
431 
432  LauAbsResonance *theResonance =
433  resonanceMaker_->getResonance(resName, resPairAmpInt, resType);
434 
435  if (theResonance == 0) {
436  cout<<"Error. Couldn't find the resonance "<<resName<<endl;
437  cout<<"Sorry. Resonance type is not recognised. Not added."<<endl;
438  return;
439  }
440 
441  // Change resonance and lineshape parameters as required.
442  if (newMass > 0.0 || newWidth > 0.0 || newSpin > -1) {
443  theResonance->changeResonance(newMass, newWidth, newSpin);
444  }
445 
446  // implement the helicity flip here
447  if (flipHelicity_ && daughters_->getCharge(resPairAmpInt) == 0) {
448  if ( daughters_->getChargeParent() == 0 && daughters_->getTypeParent() > 0 ) {
449  theResonance->flipHelicity(kTRUE);
450  }
451  }
452 
453  TString resTypeName(resType);
454 
455  // Reset the Blatt-Weisskopf barrier factors as appropriate
456  if (!resTypeName.CompareTo("RelBW")) {
457  LauRelBreitWignerRes* theRBW = dynamic_cast<LauRelBreitWignerRes*>(theResonance);
459  } else if (!resTypeName.CompareTo("GS")) {
460  LauGounarisSakuraiRes* theGS = dynamic_cast<LauGounarisSakuraiRes*>(theResonance);
462  }
463 
464  if (changeLASSScatteringLength_ == kTRUE && !resTypeName.CompareTo("LASS") ) {
465  theResonance->setResonanceParameter(LASSScatteringLength_, "a");
466  }
467  if (changeLASSEffectiveRange_ == kTRUE && !resTypeName.CompareTo("LASS") ) {
468  theResonance->setResonanceParameter(LASSEffectiveRange_, "r");
469  }
470  if (changeLASSResonanceMag_ == kTRUE && !resTypeName.CompareTo("LASS") ) {
471  theResonance->setResonanceParameter(LASSResonanceMag_, "R");
472  }
473  if (changeLASSResonancePhase_ == kTRUE && !resTypeName.CompareTo("LASS") ) {
474  theResonance->setResonanceParameter(LASSResonancePhase_, "phiR");
475  }
476  if (changeLASSBackgroundMag_ == kTRUE && !resTypeName.CompareTo("LASS") ) {
477  theResonance->setResonanceParameter(LASSBackgroundMag_, "B");
478  }
479  if (changeLASSBackgroundPhase_ == kTRUE && !resTypeName.CompareTo("LASS") ) {
480  theResonance->setResonanceParameter(LASSBackgroundPhase_, "phiB");
481  }
482  if (changeLASSCutOff_ == kTRUE && !resTypeName.CompareTo("LASS") ) {
483  theResonance->setResonanceParameter(LASSCutOff_, "cutOff");
484  }
485 
486  if (changeFlatteParameterg1_ == kTRUE && !resTypeName.CompareTo("Flatte") ) {
487  theResonance->setResonanceParameter(FlatteParameterg1_, "g1");
488  }
489  if (changeFlatteParameterg2_ == kTRUE && !resTypeName.CompareTo("Flatte") ) {
490  theResonance->setResonanceParameter(FlatteParameterg2_, "g2");
491  }
492 
493  TString resonanceName = theResonance->getResonanceName();
494 
495  if ( resonanceName.BeginsWith("BelleNR", TString::kExact) ) {
496  LauBelleNR* belleNR = dynamic_cast<LauBelleNR*>(theResonance);
497  if (belleNR != 0) {
498  belleNR->setAlpha(BelleNRAlpha_);
499  } else {
500  cout<<"Belle non-resonant object is null"<<endl;
501  }
502  } else if ( resonanceName.BeginsWith("BelleSymNR", TString::kExact) ) {
503  LauBelleSymNR* belleNR = dynamic_cast<LauBelleSymNR*>(theResonance);
504  if (belleNR != 0) {
505  belleNR->initialise(symmetricalDP_, BelleNRAlpha_, resTypeName);
506  } else {
507  cout<<"Belle non-resonant object is null"<<endl;
508  }
509  }
510 
511  // Initialise the resonance model
512  theResonance->initialise();
513 
514  // Set the resonance name and what track is the bachelor
515  resTypAmp_.push_back(resonanceName);
516  resIntAmp_.push_back(resonanceMaker_->resTypeInt(resonanceName));
517 
518  // Always force the non-resonant amplitude pair to have resPairAmp = 0
519  // in case the user chooses the wrong number.
520  if ( (resonanceName.BeginsWith("NonReson", TString::kExact) == kTRUE) ||
521  (resonanceName.BeginsWith("BelleSymNR", TString::kExact) == kTRUE) ||
522  (resonanceName.BeginsWith("NRModel", TString::kExact) == kTRUE)) {
523  cout<<"Setting resPairAmp to 0 for "<<resonanceName<<" contribution."<<endl;
524  resPairAmp_.push_back(0);
525  } else {
526  resPairAmp_.push_back(resPairAmpInt);
527  }
528 
529  // Increment the number of resonance amplitudes we have so far
530  nAmp_++;
531 
532  // Finally, add the resonance object to the internal array
533  sigResonances_.push_back(theResonance);
534 
535  cout<<"Successfully added resonance. Total number of resonances so far = "<<nAmp_<<endl;
536 
537 }
538 
539 void LauIsobarDynamics::defineKMatrixPropagator(const TString& propName, const TString& paramFileName, Int_t resPairAmpInt,
540  Int_t nChannels, Int_t nPoles, Int_t rowIndex)
541 {
542  // Define the K-matrix propagator. The resPairAmpInt integer specifies which mass combination should be used
543  // for the invariant mass-squared variable "s". The pole masses and coupling constants are defined in the
544  // paramFileName parameter file. The number of channels and poles are defined by the nChannels and nPoles integers, respectively.
545  // The integer rowIndex specifies which row of the propagator should be used when
546  // summing over all amplitude channels: S-wave will be the first row, so rowIndex = 1.
547 
548  if (rowIndex < 1) {
549  cerr<<"Error in defineKMatrixPropagator: rowIndex must be > 0 but is equal to "<<rowIndex<<endl;
550  return;
551  }
552 
553  TString propagatorName(propName), parameterFile(paramFileName);
554 
555  LauKMatrixPropagator* thePropagator = LauKMatrixPropFactory::getInstance()->getPropagator(propagatorName, parameterFile,
556  resPairAmpInt, nChannels,
557  nPoles, rowIndex);
558  kMatrixPropagators_[propagatorName] = thePropagator;
559 
560 }
561 
562 void LauIsobarDynamics::addKMatrixProdPole(const TString& poleName, const TString& propName, Int_t poleIndex)
563 {
564 
565  // Add a K-matrix production pole term, using the K-matrix propagator given by the propName.
566  // Here, poleIndex is the integer specifying the pole number.
567 
568  // First, find the K-matrix propagator.
569  KMPropMap::iterator mapIter = kMatrixPropagators_.find(propName);
570  if (mapIter != kMatrixPropagators_.end()) {
571 
572  LauKMatrixPropagator* thePropagator = mapIter->second;
573 
574  // Make sure the pole index is valid
575  Int_t nPoles = thePropagator->getNPoles();
576  if (poleIndex < 1 || poleIndex > nPoles) {
577  cerr<<"Error in LauIsobarDynamics::addKMatrixProdPole. The pole index "<<poleIndex
578  <<" is not between 1 and "<<nPoles<<". Not adding production pole "<<poleName
579  <<" for K-matrix propagator "<<propName<<endl;
580  return;
581  }
582 
583  // Now add the K-matrix production pole amplitude to the vector of LauAbsResonance pointers.
584  Int_t resPairAmpInt = thePropagator->getResPairAmpInt();
585  LauAbsResonance* prodPole = new LauKMatrixProdPole(poleName, poleIndex, resPairAmpInt, thePropagator, daughters_);
586 
587  resTypAmp_.push_back(poleName);
588  resIntAmp_.push_back(0);
589  resPairAmp_.push_back(resPairAmpInt);
590 
591  nAmp_++;
592  sigResonances_.push_back(prodPole);
593 
594  // Also store the propName-poleName pair for calculating total fit fractions later on
595  // (avoiding the need to use dynamic casts to check which resonances are of the K-matrix type)
596  kMatrixPropSet_[poleName] = propName;
597 
598  cout<<"Successfully added K-matrix production pole term. Total number of resonances so far = "<<nAmp_<<endl;
599 
600  } else {
601 
602  cerr<<"Error in LauIsobarDynamics::addKMatrixProdPole. The propagator of the name "<<propName
603  <<" could not be found for the production pole "<<poleName<<endl;
604 
605  }
606 
607 }
608 
609 
610 void LauIsobarDynamics::addKMatrixProdSVP(const TString& SVPName, const TString& propName, Int_t channelIndex)
611 {
612 
613  // Add a K-matrix production "slowly-varying part" (SVP) term, using the K-matrix propagator
614  // given by the propName. Here, channelIndex is the integer specifying the channel number.
615 
616  // First, find the K-matrix propagator.
617  KMPropMap::iterator mapIter = kMatrixPropagators_.find(propName);
618  if (mapIter != kMatrixPropagators_.end()) {
619 
620  LauKMatrixPropagator* thePropagator = mapIter->second;
621 
622  // Make sure the channel index is valid
623  Int_t nChannels = thePropagator->getNChannels();
624  if (channelIndex < 1 || channelIndex > nChannels) {
625  cerr<<"Error in LauIsobarDynamics::addKMatrixProdSVP. The channel index "<<channelIndex
626  <<" is not between 1 and "<<nChannels<<". Not adding production slowly-varying part "<<SVPName
627  <<" for K-matrix propagator "<<propName<<endl;
628  return;
629  }
630 
631  // Now add the K-matrix production SVP amplitude to the vector of LauAbsResonance pointers.
632  Int_t resPairAmpInt = thePropagator->getResPairAmpInt();
633  LauAbsResonance* prodSVP = new LauKMatrixProdSVP(SVPName, channelIndex, resPairAmpInt, thePropagator, daughters_);
634 
635  resTypAmp_.push_back(SVPName);
636  resIntAmp_.push_back(0);
637  resPairAmp_.push_back(resPairAmpInt);
638 
639  nAmp_++;
640  sigResonances_.push_back(prodSVP);
641 
642  // Also store the SVPName-propName pair for calculating total fit fractions later on
643  // (avoiding the need to use dynamic casts to check which resonances are of the K-matrix type)
644  kMatrixPropSet_[SVPName] = propName;
645 
646  cout<<"Successfully added K-matrix production slowly-varying (SVP) term. Total number of resonances so far = "
647  <<nAmp_<<endl;
648 
649  } else {
650 
651  cerr<<"Error in LauIsobarDynamics::addKMatrixProdSVP. The propagator of the name "<<propName
652  <<" could not be found for the production slowly-varying part "<<SVPName<<endl;
653 
654  }
655 
656 }
657 
659 {
660  TString testName(name);
661  testName.ToLower();
662  LauAbsResonance* theResonance(0);
663 
664  for (std::vector<LauAbsResonance*>::iterator iter=sigResonances_.begin(); iter!=sigResonances_.end(); ++iter) {
665  theResonance = (*iter);
666  if (theResonance != 0) {
667  TString resString = theResonance->getResonanceName();
668  resString.ToLower();
669  if (resString.BeginsWith(testName, TString::kExact)) {
670  return theResonance;
671  }
672  }
673  }
674 
675  cout<<"ERROR in LauIsobarDynamics::findResonance : Couldn't find resonance \""<<name<<"\" in the model."<<endl;
676  return 0;
677 }
678 
680 {
681  TString testName(name);
682  testName.ToLower();
683  const LauAbsResonance* theResonance(0);
684 
685  for (std::vector<LauAbsResonance*>::const_iterator iter=sigResonances_.begin(); iter!=sigResonances_.end(); ++iter) {
686  theResonance = (*iter);
687  if (theResonance != 0) {
688  TString resString = theResonance->getResonanceName();
689  resString.ToLower();
690  if (resString.BeginsWith(testName, TString::kExact)) {
691  return theResonance;
692  }
693  }
694  }
695 
696  cout<<"ERROR in LauIsobarDynamics::findResonance : Couldn't find resonance \""<<name<<"\" in the model."<<endl;
697  return 0;
698 }
699 
700 void LauIsobarDynamics::removeCharge(TString& string) const
701 {
702  Ssiz_t index = string.Index("+");
703  if (index != -1) {
704  string.Remove(index,1);
705  }
706  index = string.Index("-");
707  if (index != -1) {
708  string.Remove(index,1);
709  }
710 }
711 
712 void LauIsobarDynamics::changeResonance(const TString& resName, Double_t newMass, Double_t newWidth, Int_t newSpin)
713 {
714  // Change the mass, width or spin of a resonance.
715 
716  if (newMass > 0.0 || newWidth > 0.0 || newSpin > -1) {
717  cout<<"ERROR in LauIsobarDynamics::changeResonance : mass, width and spin parameters all out of range."<<endl;
718  return;
719  }
720 
721  LauAbsResonance* theRes = this->findResonance(resName);
722  if (theRes != 0) {
723  theRes->changeResonance(newMass, newWidth, newSpin);
724  }
725 
726 }
727 
729 {
730  // Use Gauss-Legendre quadrature integration
731 
732  // Get the rectangle that encloses the DP
733  Double_t minm13 = kinematics_->getm13Min();
734  Double_t maxm13 = kinematics_->getm13Max();
735  Double_t minm23 = kinematics_->getm23Min();
736  Double_t maxm23 = kinematics_->getm23Max();
737  Double_t minm12 = kinematics_->getm12Min();
738  Double_t maxm12 = kinematics_->getm12Max();
739 
740  // Find out whether we have narrow resonances in the DP (defined here as width < 20 MeV).
741  std::multimap<Double_t,Double_t> m13NarrowRes;
742  std::multimap<Double_t,Double_t> m23NarrowRes;
743  std::multimap<Double_t,Double_t> m12NarrowRes;
744  for ( std::vector<LauAbsResonance*>::const_iterator iter = sigResonances_.begin(); iter != sigResonances_.end(); ++iter ) {
745  Double_t width = (*iter)->getWidth();
746  if ( width > 0.020 || width == 0.0 ) { continue; }
747  Double_t mass = (*iter)->getMass();
748  Int_t pair = (*iter)->getPairInt();
749  TString name = (*iter)->getResonanceName();
750  cout<<"Found narrow resonance: "<<name<<", mass = "<<mass<<", width = "<<width<<", pair int = "<<pair<<endl;
751  if ( pair == 1 ) {
752  if ( mass < minm23 || mass > maxm23 ){ continue; }
753  m23NarrowRes.insert( std::make_pair(width,mass) );
754  } else if ( pair == 2 ) {
755  if ( mass < minm13 || mass > maxm13 ){ continue; }
756  m13NarrowRes.insert( std::make_pair(width,mass) );
757  } else if ( pair == 3 ) {
758  if ( mass < minm12 || mass > maxm12 ){ continue; }
759  m12NarrowRes.insert( std::make_pair(width,mass) );
760  } else {
761  cerr<<"WARNING in LauIsobarDynamics::calcDPNormalisation : strange pair integer, "<<pair<<", for resonance \""<<(*iter)->getResonanceName()<<endl;
762  }
763  }
764 
765  // Find the narrowest resonance in each mass pairing
766  const Bool_t e12 = m12NarrowRes.empty();
767  const Bool_t e13 = m13NarrowRes.empty();
768  const Bool_t e23 = m23NarrowRes.empty();
769  //const UInt_t s12 = e12 ? 0 : m12NarrowRes.size();
770  const UInt_t s13 = e13 ? 0 : m13NarrowRes.size();
771  const UInt_t s23 = e23 ? 0 : m23NarrowRes.size();
772  const Double_t w12 = e12 ? DBL_MAX : m12NarrowRes.begin()->first / 100.0;
773  const Double_t w13 = e13 ? DBL_MAX : m13NarrowRes.begin()->first / 100.0;
774  const Double_t w23 = e23 ? DBL_MAX : m23NarrowRes.begin()->first / 100.0;
775 
776  // Start off with default bin width (5 MeV)
777  Double_t m13BinWidth = m13BinWidth_;
778  Double_t m23BinWidth = m23BinWidth_;
779 
780  // Depending on how many narrow resonances we have and where they
781  // are we adopt different approaches
782  if ( e12 && e13 && e23 ) {
783  // If we have no narrow resonances just integrate the whole
784  // DP with the standard bin widths
785  cout<<"No narrow resonances found, integrating over whole Dalitz plot..."<<endl;
786  this->calcDPPartialIntegral(minm13, maxm13, minm23, maxm23, m13BinWidth, m23BinWidth);
787  } else if ( ! e12 ) {
788  // If we have a narrow resonance on the diagonal then we'll have to
789  // just use a narrow bin width over the whole DP (1/10 of the width
790  // of the narrowest resonance in any mass pair)
791  m13BinWidth = m23BinWidth = 10.0*TMath::Min( w12, TMath::Min( w13, w23 ) );
792  cout<<"One or more narrow resonances found in m12, integrating over whole Dalitz plot with bin width of "<<m13BinWidth<<" GeV/c2..."<<endl;
793  this->calcDPPartialIntegral(minm13, maxm13, minm23, maxm23, m13BinWidth, m23BinWidth);
794  } else if ( s13==1 && e23 ) {
795  // We have a single narrow resonance in m13
796  // Divide the plot into 3 regions: the resonance band and
797  // the two areas either side.
798  Double_t mass = m13NarrowRes.begin()->second;
799  Double_t width = m13NarrowRes.begin()->first;
800  Double_t resMin = mass - 5.0*width;
801  Double_t resMax = mass + 5.0*width;
802  // if the resonance is close to threshold just go from
803  // threshold to resMax, otherwise treat threshold to resMin
804  // as a separate region
805  if ( resMin < (minm13+50.0*m13BinWidth_) ) {
806  cout<<"One narrow resonance found in m13, close to threshold, dividing Dalitz plot into two regions..."<<endl;
807  this->calcDPPartialIntegral(resMax, maxm13, minm23, maxm23, m13BinWidth, m23BinWidth);
808  m13BinWidth = w13;
809  this->calcDPPartialIntegral(minm13, resMax, minm23, maxm23, m13BinWidth, m23BinWidth);
810  } else {
811  cout<<"One narrow resonance found in m13, dividing Dalitz plot into three regions..."<<endl;
812  this->calcDPPartialIntegral(minm13, resMin, minm23, maxm23, m13BinWidth, m23BinWidth);
813  this->calcDPPartialIntegral(resMax, maxm13, minm23, maxm23, m13BinWidth, m23BinWidth);
814  m13BinWidth = w13;
815  this->calcDPPartialIntegral(resMin, resMax, minm23, maxm23, m13BinWidth, m23BinWidth);
816  }
817  } else if ( s13==2 && e23 ) {
818  // We have a two narrow resonances in m13
819  // Divide the plot into 5 regions: the resonance bands,
820  // the two areas either side and the area in between.
821  std::multimap<Double_t,Double_t> massordered;
822  for ( std::map<Double_t,Double_t>::const_iterator iter = m13NarrowRes.begin(); iter != m13NarrowRes.end(); ++iter ) {
823  massordered.insert( std::make_pair( iter->second, iter->first ) );
824  }
825  Double_t res1Mass = massordered.begin()->first;
826  Double_t res1Width = massordered.begin()->second;
827  Double_t res1Min = res1Mass - 5.0*res1Width;
828  Double_t res1Max = res1Mass + 5.0*res1Width;
829  Double_t res2Mass = (++(massordered.begin()))->first;
830  Double_t res2Width = (++(massordered.begin()))->second;
831  Double_t res2Min = res2Mass - 5.0*res2Width;
832  Double_t res2Max = res2Mass + 5.0*res2Width;
833  // if the resonance is close to threshold just go from
834  // threshold to resMax, otherwise treat threshold to resMin
835  // as a separate region
836  if ( res1Min < (minm13+50.0*m13BinWidth_) ) {
837  cout<<"Two narrow resonances found in m13, one close to threshold, dividing Dalitz plot into four regions..."<<endl;
838  this->calcDPPartialIntegral(res1Max, res2Min, minm23, maxm23, m13BinWidth, m23BinWidth);
839  this->calcDPPartialIntegral(res2Max, maxm13, minm23, maxm23, m13BinWidth, m23BinWidth);
840  m13BinWidth = res1Width/100.0;
841  this->calcDPPartialIntegral(minm13, res1Max, minm23, maxm23, m13BinWidth, m23BinWidth);
842  m13BinWidth = res2Width/100.0;
843  this->calcDPPartialIntegral(res2Min, res2Max, minm23, maxm23, m13BinWidth, m23BinWidth);
844  } else {
845  cout<<"Two narrow resonances found in m13, dividing Dalitz plot into five regions..."<<endl;
846  this->calcDPPartialIntegral(minm13, res1Min, minm23, maxm23, m13BinWidth, m23BinWidth);
847  this->calcDPPartialIntegral(res1Max, res2Min, minm23, maxm23, m13BinWidth, m23BinWidth);
848  this->calcDPPartialIntegral(res2Max, maxm13, minm23, maxm23, m13BinWidth, m23BinWidth);
849  m13BinWidth = res1Width/100.0;
850  this->calcDPPartialIntegral(res1Min, res1Max, minm23, maxm23, m13BinWidth, m23BinWidth);
851  m13BinWidth = res2Width/100.0;
852  this->calcDPPartialIntegral(res2Min, res2Max, minm23, maxm23, m13BinWidth, m23BinWidth);
853  }
854  } else if ( s23==1 && e13 ) {
855  // We have a single narrow resonance in m23
856  // Divide the plot into 3 regions: the resonance band and
857  // the two areas either side.
858  Double_t mass = m23NarrowRes.begin()->second;
859  Double_t width = m23NarrowRes.begin()->first;
860  Double_t resMin = mass - 5.0*width;
861  Double_t resMax = mass + 5.0*width;
862  // if the resonance is close to threshold just go from
863  // threshold to resMax, otherwise treat threshold to resMin
864  // as a separate region
865  if ( resMin < (minm23+50.0*m23BinWidth_) ) {
866  cout<<"One narrow resonance found in m23, close to threshold, dividing Dalitz plot into two regions..."<<endl;
867  this->calcDPPartialIntegral(minm13, maxm13, resMax, maxm23, m13BinWidth, m23BinWidth);
868  m23BinWidth = w23;
869  this->calcDPPartialIntegral(minm13, maxm13, minm23, resMax, m13BinWidth, m23BinWidth);
870  } else {
871  cout<<"One narrow resonance found in m23, dividing Dalitz plot into three regions..."<<endl;
872  this->calcDPPartialIntegral(minm13, maxm13, minm23, resMin, m13BinWidth, m23BinWidth);
873  this->calcDPPartialIntegral(minm13, maxm13, resMax, maxm23, m13BinWidth, m23BinWidth);
874  m23BinWidth = w23;
875  this->calcDPPartialIntegral(minm13, maxm13, resMin, resMax, m13BinWidth, m23BinWidth);
876  }
877  } else if ( s23==2 && e13 ) {
878  // We have a two narrow resonances in m23
879  // Divide the plot into 5 regions: the resonance bands,
880  // the two areas either side and the area in between.
881  std::multimap<Double_t,Double_t> massordered;
882  for ( std::map<Double_t,Double_t>::const_iterator iter = m23NarrowRes.begin(); iter != m23NarrowRes.end(); ++iter ) {
883  massordered.insert( std::make_pair( iter->second, iter->first ) );
884  }
885  Double_t res1Mass = massordered.begin()->first;
886  Double_t res1Width = massordered.begin()->second;
887  Double_t res1Min = res1Mass - 5.0*res1Width;
888  Double_t res1Max = res1Mass + 5.0*res1Width;
889  Double_t res2Mass = (++(massordered.begin()))->first;
890  Double_t res2Width = (++(massordered.begin()))->second;
891  Double_t res2Min = res2Mass - 5.0*res2Width;
892  Double_t res2Max = res2Mass + 5.0*res2Width;
893  // if the resonance is close to threshold just go from
894  // threshold to resMax, otherwise treat threshold to resMin
895  // as a separate region
896  if ( res1Min < (minm23+50.0*m23BinWidth_) ) {
897  cout<<"Two narrow resonances found in m23, one close to threshold, dividing Dalitz plot into four regions..."<<endl;
898  this->calcDPPartialIntegral(minm13, maxm13, res1Max, res2Min, m13BinWidth, m23BinWidth);
899  this->calcDPPartialIntegral(minm13, maxm13, res2Max, maxm23, m13BinWidth, m23BinWidth);
900  m23BinWidth = res1Width/100.0;
901  this->calcDPPartialIntegral(minm13, maxm13, minm23, res1Max, m13BinWidth, m23BinWidth);
902  m23BinWidth = res2Width/100.0;
903  this->calcDPPartialIntegral(minm13, maxm13, res2Min, res2Max, m13BinWidth, m23BinWidth);
904  } else {
905  cout<<"Two narrow resonances found in m23, dividing Dalitz plot into five regions..."<<endl;
906  this->calcDPPartialIntegral(minm13, maxm13, minm23, res1Min, m13BinWidth, m23BinWidth);
907  this->calcDPPartialIntegral(minm13, maxm13, res1Max, res2Min, m13BinWidth, m23BinWidth);
908  this->calcDPPartialIntegral(minm13, maxm13, res2Max, maxm23, m13BinWidth, m23BinWidth);
909  m23BinWidth = res1Width/100.0;
910  this->calcDPPartialIntegral(minm13, maxm13, res1Min, res1Max, m13BinWidth, m23BinWidth);
911  m23BinWidth = res2Width/100.0;
912  this->calcDPPartialIntegral(minm13, maxm13, res2Min, res2Max, m13BinWidth, m23BinWidth);
913  }
914  } else if ( s13==1 && s23==1 ) {
915  // We have a single narrow resonance in both m13 and m23
916  // Divide the plot into 9 regions: the point where the
917  // resonance bands cross, the four other parts of the bands
918  // and the four remaining areas of the DP.
919  Double_t mass23 = m23NarrowRes.begin()->second;
920  Double_t width23 = m23NarrowRes.begin()->first;
921  Double_t resMin23 = mass23 - 5.0*width23;
922  Double_t resMax23 = mass23 + 5.0*width23;
923  Double_t mass13 = m13NarrowRes.begin()->second;
924  Double_t width13 = m13NarrowRes.begin()->first;
925  Double_t resMin13 = mass13 - 5.0*width13;
926  Double_t resMax13 = mass13 + 5.0*width13;
927  // if either resonance is close to threshold just go from
928  // threshold to resMax, otherwise treat threshold to resMin
929  // as a separate region
930  if ( resMin13 < (minm13+50.0*m13BinWidth_) && resMin23 < (minm23+50.0*m23BinWidth_) ) {
931  cout<<"One narrow resonance found in m13 and one in m23, both close to threshold, dividing Dalitz plot into four regions..."<<endl;
932  m13BinWidth = m13BinWidth_;
933  m23BinWidth = m23BinWidth_;
934  this->calcDPPartialIntegral(resMax13, maxm13, resMax23, maxm23, m13BinWidth, m23BinWidth);
935  m13BinWidth = m13BinWidth_;
936  m23BinWidth = w23;
937  this->calcDPPartialIntegral(resMax13, maxm13, minm23, resMax23, m13BinWidth, m23BinWidth);
938  m13BinWidth = w13;
939  m23BinWidth = m23BinWidth_;
940  this->calcDPPartialIntegral(minm13, resMax13, resMax23, maxm23, m13BinWidth, m23BinWidth);
941  m13BinWidth = w13;
942  m23BinWidth = w23;
943  this->calcDPPartialIntegral(minm13, resMax13, minm23, resMax23, m13BinWidth, m23BinWidth);
944  } else if ( resMin13 < (minm13+50.0*m13BinWidth_) ) {
945  cout<<"One narrow resonance found in m13, close to threshold, and one in m23, not close to threshold, dividing Dalitz plot into six regions..."<<endl;
946  this->calcDPPartialIntegral(resMax13, maxm13, minm23, resMin23, m13BinWidth, m23BinWidth);
947  this->calcDPPartialIntegral(resMax13, maxm13, resMax23, maxm23, m13BinWidth, m23BinWidth);
948  m13BinWidth = m13BinWidth_;
949  m23BinWidth = w23;
950  this->calcDPPartialIntegral(resMax13, maxm13, resMin23, resMax23, m13BinWidth, m23BinWidth);
951  m13BinWidth = w13;
952  m23BinWidth = m23BinWidth_;
953  this->calcDPPartialIntegral(minm13, resMax13, minm23, resMin23, m13BinWidth, m23BinWidth);
954  this->calcDPPartialIntegral(minm13, resMax13, resMax23, maxm23, m13BinWidth, m23BinWidth);
955  m13BinWidth = w13;
956  m23BinWidth = w23;
957  this->calcDPPartialIntegral(minm13, resMax13, resMin23, resMax23, m13BinWidth, m23BinWidth);
958  } else if ( resMin23 < (minm23+50.0*m23BinWidth_) ) {
959  cout<<"One narrow resonance found in m23, close to threshold, and one in m13, not close to threshold, dividing Dalitz plot into six regions..."<<endl;
960  this->calcDPPartialIntegral(minm13, resMin13, resMax23, maxm23, m13BinWidth, m23BinWidth);
961  this->calcDPPartialIntegral(resMax13, maxm13, resMax23, maxm23, m13BinWidth, m23BinWidth);
962  m13BinWidth = m13BinWidth_;
963  m23BinWidth = w23;
964  this->calcDPPartialIntegral(minm13, resMin13, minm23, resMax23, m13BinWidth, m23BinWidth);
965  this->calcDPPartialIntegral(resMax13, maxm13, minm23, resMax23, m13BinWidth, m23BinWidth);
966  m13BinWidth = w13;
967  m23BinWidth = m23BinWidth_;
968  this->calcDPPartialIntegral(resMin13, resMax13, resMax23, maxm23, m13BinWidth, m23BinWidth);
969  m13BinWidth = w13;
970  m23BinWidth = w23;
971  this->calcDPPartialIntegral(resMin13, resMax13, minm23, resMax23, m13BinWidth, m23BinWidth);
972  } else {
973  cout<<"One narrow resonance found in both m13 and m23, neither close to threshold, dividing Dalitz plot into nine regions..."<<endl;
974  this->calcDPPartialIntegral(minm13, resMin13, minm23, resMin23, m13BinWidth, m23BinWidth);
975  this->calcDPPartialIntegral(minm13, resMin13, resMax23, maxm23, m13BinWidth, m23BinWidth);
976  this->calcDPPartialIntegral(resMax13, maxm13, minm23, resMin23, m13BinWidth, m23BinWidth);
977  this->calcDPPartialIntegral(resMax13, maxm13, resMax23, maxm23, m13BinWidth, m23BinWidth);
978  m13BinWidth = m13BinWidth_;
979  m23BinWidth = w23;
980  this->calcDPPartialIntegral(minm13, resMin13, resMin23, resMax23, m13BinWidth, m23BinWidth);
981  this->calcDPPartialIntegral(resMax13, maxm13, resMin23, resMax23, m13BinWidth, m23BinWidth);
982  m13BinWidth = w13;
983  m23BinWidth = m23BinWidth_;
984  this->calcDPPartialIntegral(resMin13, resMax13, minm23, resMin23, m13BinWidth, m23BinWidth);
985  this->calcDPPartialIntegral(resMin13, resMax13, resMax23, maxm23, m13BinWidth, m23BinWidth);
986  m13BinWidth = w13;
987  m23BinWidth = w23;
988  this->calcDPPartialIntegral(resMin13, resMax13, resMin23, resMax23, m13BinWidth, m23BinWidth);
989  }
990  } else if ( e23 && s13>1 ) {
991  // We have multiple narrow resonances in m13 only.
992  // Divide the plot into 2 regions: threshold to the most
993  // massive of the narrow resonances, and the rest
994  cout<<"Multiple narrow resonances found in m13, dividing Dalitz plot into two regions..."<<endl;
995  Double_t mass = 0.0;
996  Double_t width = 0.0;
997  for ( std::map<Double_t,Double_t>::const_iterator iter = m13NarrowRes.begin(); iter != m13NarrowRes.end(); ++iter ) {
998  if ( mass < iter->second ) {
999  mass = iter->second;
1000  width = iter->first;
1001  }
1002  }
1003  Double_t resMax = mass + 5.0*width;
1004  this->calcDPPartialIntegral(resMax, maxm13, minm23, maxm23, m13BinWidth, m23BinWidth);
1005  m13BinWidth = w13;
1006  this->calcDPPartialIntegral(minm13, resMax, minm23, maxm23, m13BinWidth, m23BinWidth);
1007  } else if ( e13 && s23>1 ) {
1008  // We have multiple narrow resonances in m23 only.
1009  // Divide the plot into 2 regions: threshold to the most
1010  // massive of the narrow resonances, and the rest
1011  cout<<"Multiple narrow resonances found in m23, dividing Dalitz plot into two regions..."<<endl;
1012  Double_t mass = 0.0;
1013  Double_t width = 0.0;
1014  for ( std::map<Double_t,Double_t>::const_iterator iter = m23NarrowRes.begin(); iter != m23NarrowRes.end(); ++iter ) {
1015  if ( mass < iter->second ) {
1016  mass = iter->second;
1017  width = iter->first;
1018  }
1019  }
1020  Double_t resMax = mass + 5.0*width;
1021  this->calcDPPartialIntegral(minm13, maxm13, resMax, maxm23, m13BinWidth, m23BinWidth);
1022  m23BinWidth = w23;
1023  this->calcDPPartialIntegral(minm13, maxm13, minm23, resMax, m13BinWidth, m23BinWidth);
1024  } else if ( s13==1 && s23>1 ) {
1025  // We've got a single narrow resonance in m13 and multiple
1026  // narrow resonances in m23. Divide the plot into 6 regions.
1027  Double_t mass23 = 0.0;
1028  Double_t width23 = 0.0;
1029  for ( std::map<Double_t,Double_t>::const_iterator iter = m23NarrowRes.begin(); iter != m23NarrowRes.end(); ++iter ) {
1030  if ( mass23 < iter->second ) {
1031  mass23 = iter->second;
1032  width23 = iter->first;
1033  }
1034  }
1035  Double_t resMax23 = mass23 + 5.0*width23;
1036  Double_t mass13 = m13NarrowRes.begin()->second;
1037  Double_t width13 = m13NarrowRes.begin()->first;
1038  Double_t resMin13 = mass13 - 5.0*width13;
1039  Double_t resMax13 = mass13 + 5.0*width13;
1040  // if the m13 resonance is close to threshold just go from
1041  // threshold to resMax, otherwise treat threshold to resMin
1042  // as a separate region
1043  if ( resMin13 < (minm13+50.0*m13BinWidth_) ) {
1044  cout<<"Multiple narrow resonances found in m23 and one in m13, close to threshold, dividing Dalitz plot into four regions..."<<endl;
1045  m13BinWidth = m13BinWidth_;
1046  m23BinWidth = m23BinWidth_;
1047  this->calcDPPartialIntegral(resMax13, maxm13, resMax23, maxm23, m13BinWidth, m23BinWidth);
1048  m13BinWidth = m13BinWidth_;
1049  m23BinWidth = w23;
1050  this->calcDPPartialIntegral(resMax13, maxm13, minm23, resMax23, m13BinWidth, m23BinWidth);
1051  m13BinWidth = w13;
1052  m23BinWidth = m23BinWidth_;
1053  this->calcDPPartialIntegral(minm13, resMax13, resMax23, maxm23, m13BinWidth, m23BinWidth);
1054  m13BinWidth = w13;
1055  m23BinWidth = w23;
1056  this->calcDPPartialIntegral(minm13, resMax13, minm23, resMax23, m13BinWidth, m23BinWidth);
1057  } else {
1058  cout<<"Multiple narrow resonances found in m23 and one in m13, not close to threshold, dividing Dalitz plot into six regions..."<<endl;
1059  m13BinWidth = m13BinWidth_;
1060  m23BinWidth = m23BinWidth_;
1061  this->calcDPPartialIntegral(minm13, resMin13, resMax23, maxm23, m13BinWidth, m23BinWidth);
1062  this->calcDPPartialIntegral(resMax13, maxm13, resMax23, maxm23, m13BinWidth, m23BinWidth);
1063  m13BinWidth = m13BinWidth_;
1064  m23BinWidth = w23;
1065  this->calcDPPartialIntegral(minm13, resMin13, minm23, resMax23, m13BinWidth, m23BinWidth);
1066  this->calcDPPartialIntegral(resMax13, maxm13, minm23, resMax23, m13BinWidth, m23BinWidth);
1067  m13BinWidth = w13;
1068  m23BinWidth = m23BinWidth_;
1069  this->calcDPPartialIntegral(resMin13, resMax13, resMax23, maxm23, m13BinWidth, m23BinWidth);
1070  m13BinWidth = w13;
1071  m23BinWidth = w23;
1072  this->calcDPPartialIntegral(resMin13, resMax13, minm23, resMax23, m13BinWidth, m23BinWidth);
1073  }
1074  } else if ( s13>1 && s23==1 ) {
1075  // We've got a single narrow resonance in m23 and multiple
1076  // narrow resonances in m13. Divide the plot into 6 regions.
1077  Double_t mass13 = 0.0;
1078  Double_t width13 = 0.0;
1079  for ( std::map<Double_t,Double_t>::const_iterator iter = m13NarrowRes.begin(); iter != m13NarrowRes.end(); ++iter ) {
1080  if ( mass13 < iter->second ) {
1081  mass13 = iter->second;
1082  width13 = iter->first;
1083  }
1084  }
1085  Double_t resMax13 = mass13 + 5.0*width13;
1086  Double_t mass23 = m23NarrowRes.begin()->second;
1087  Double_t width23 = m23NarrowRes.begin()->first;
1088  Double_t resMin23 = mass23 - 5.0*width23;
1089  Double_t resMax23 = mass23 + 5.0*width23;
1090  // if the m23 resonance is close to threshold just go from
1091  // threshold to resMax, otherwise treat threshold to resMin
1092  // as a separate region
1093  if ( resMin23 < (minm23+50.0*m23BinWidth_) ) {
1094  cout<<"Multiple narrow resonances found in m13 and one in m23, close to threshold, dividing Dalitz plot into four regions..."<<endl;
1095  m13BinWidth = m13BinWidth_;
1096  m23BinWidth = m23BinWidth_;
1097  this->calcDPPartialIntegral(resMax13, maxm13, resMax23, maxm23, m13BinWidth, m23BinWidth);
1098  m13BinWidth = m13BinWidth_;
1099  m23BinWidth = w23;
1100  this->calcDPPartialIntegral(resMax13, maxm13, minm23, resMax23, m13BinWidth, m23BinWidth);
1101  m13BinWidth = w13;
1102  m23BinWidth = m23BinWidth_;
1103  this->calcDPPartialIntegral(minm13, resMax13, resMax23, maxm23, m13BinWidth, m23BinWidth);
1104  m13BinWidth = w13;
1105  m23BinWidth = w23;
1106  this->calcDPPartialIntegral(minm13, resMax13, minm23, resMax23, m13BinWidth, m23BinWidth);
1107  } else {
1108  cout<<"Multiple narrow resonances found in m13 and one in m23, not close to threshold, dividing Dalitz plot into six regions..."<<endl;
1109  this->calcDPPartialIntegral(resMax13, maxm13, minm23, resMin23, m13BinWidth, m23BinWidth);
1110  this->calcDPPartialIntegral(resMax13, maxm13, resMax23, maxm23, m13BinWidth, m23BinWidth);
1111  m13BinWidth = m13BinWidth_;
1112  m23BinWidth = w23;
1113  this->calcDPPartialIntegral(resMax13, maxm13, resMin23, resMax23, m13BinWidth, m23BinWidth);
1114  m13BinWidth = w13;
1115  m23BinWidth = m23BinWidth_;
1116  this->calcDPPartialIntegral(minm13, resMax13, minm23, resMin23, m13BinWidth, m23BinWidth);
1117  this->calcDPPartialIntegral(minm13, resMax13, resMax23, maxm23, m13BinWidth, m23BinWidth);
1118  m13BinWidth = w13;
1119  m23BinWidth = w23;
1120  this->calcDPPartialIntegral(minm13, resMax13, resMin23, resMax23, m13BinWidth, m23BinWidth);
1121  }
1122  } else {
1123  // We've got multiple narrow resonances in both m13 and m23.
1124  // Divide the plot into 4 regions.
1125  cout<<"Multiple narrow resonances found in both m13 and m23, dividing Dalitz plot into four regions..."<<endl;
1126  Double_t mass13 = 0.0;
1127  Double_t width13 = 0.0;
1128  for ( std::map<Double_t,Double_t>::const_iterator iter = m13NarrowRes.begin(); iter != m13NarrowRes.end(); ++iter ) {
1129  if ( mass13 < iter->second ) {
1130  mass13 = iter->second;
1131  width13 = iter->first;
1132  }
1133  }
1134  Double_t resMax13 = mass13 + 5.0*width13;
1135 
1136  Double_t mass23 = 0.0;
1137  Double_t width23 = 0.0;
1138  for ( std::map<Double_t,Double_t>::const_iterator iter = m23NarrowRes.begin(); iter != m23NarrowRes.end(); ++iter ) {
1139  if ( mass23 < iter->second ) {
1140  mass23 = iter->second;
1141  width23 = iter->first;
1142  }
1143  }
1144  Double_t resMax23 = mass23 + 5.0*width23;
1145 
1146  m13BinWidth = m13BinWidth_;
1147  m23BinWidth = m23BinWidth_;
1148  this->calcDPPartialIntegral(resMax13, maxm13, resMax23, maxm23, m13BinWidth, m23BinWidth);
1149  m13BinWidth = m13BinWidth_;
1150  m23BinWidth = w23;
1151  this->calcDPPartialIntegral(resMax13, maxm13, minm23, resMax23, m13BinWidth, m23BinWidth);
1152  m13BinWidth = w13;
1153  m23BinWidth = m23BinWidth_;
1154  this->calcDPPartialIntegral(minm13, resMax13, resMax23, maxm23, m13BinWidth, m23BinWidth);
1155  m13BinWidth = w13;
1156  m23BinWidth = w23;
1157  this->calcDPPartialIntegral(minm13, resMax13, minm23, resMax23, m13BinWidth, m23BinWidth);
1158  }
1159 }
1160 
1161 void LauIsobarDynamics::setIntegralBinWidths(Double_t m13BinWidth, Double_t m23BinWidth)
1162 {
1163  // Specify whether we're going to use Gauss-Legendre integration to calculate the normalisation
1164  // integrals, and the bin widths we require for the m13 and m23 axes. Note that the integration
1165  // is done over m13, m23 space, with the appropriate Jacobian applied, and not m13^2, m23^2 space.
1166  // The default bin widths in m13 and m23 space are 5 MeV.
1167 
1168  m13BinWidth_ = m13BinWidth;
1169  m23BinWidth_ = m23BinWidth;
1170 }
1171 
1172 void LauIsobarDynamics::calcDPPartialIntegral(Double_t minm13, Double_t maxm13, Double_t minm23, Double_t maxm23,
1173  Double_t m13BinWidth, Double_t m23BinWidth)
1174 {
1175  // Calculate the total DP area, as well as finding the normalisation terms for
1176  // the signal resonances
1177 
1178  Int_t i(0), j(0);
1179  Double_t precision(1e-6);
1180 
1181  Double_t meanm13 = 0.5*(minm13 + maxm13);
1182  Double_t rangem13 = maxm13 - minm13;
1183  Double_t halfRangem13 = 0.5*rangem13;
1184 
1185  Double_t meanm23 = 0.5*(minm23 + maxm23);
1186  Double_t rangem23 = maxm23 - minm23;
1187  Double_t halfRangem23 = 0.5*rangem23;
1188 
1189  Double_t intFactor = halfRangem13*halfRangem23;
1190 
1191  // Choose smallest of mass ranges to set number of bins, given specified bin width
1192  Int_t nm13Points = static_cast<Int_t>((rangem13/m13BinWidth));
1193  Int_t nm23Points = static_cast<Int_t>((rangem23/m23BinWidth));
1194 
1195  // Avoid integral if we have no points in either x or y space
1196  if (nm13Points == 0 || nm23Points == 0) {return;}
1197 
1198  cout<<"nm13Points = "<<nm13Points<<", nm23Points = "<<nm23Points<<endl;
1199  cout<<"m13BinWidth = "<<m13BinWidth<<", m23BinWidth = "<<m23BinWidth<<endl;
1200  cout<<"Integrating over m13 = "<<minm13<<" to "<<maxm13<<", m23 = "<<minm23<<" to "<<maxm23<<endl;
1201 
1202  LauIntegrals DPIntegrals(precision);
1203  std::vector<Double_t> m13Weights, m23Weights;
1204  std::vector<Double_t> m13Abscissas, m23Abscissas;
1205 
1206  DPIntegrals.calcGaussLegendreWeights(nm13Points, m13Abscissas, m13Weights);
1207  DPIntegrals.calcGaussLegendreWeights(nm23Points, m23Abscissas, m23Weights);
1208 
1209  Int_t nm13Weights = static_cast<Int_t>(m13Weights.size());
1210  Int_t nm23Weights = static_cast<Int_t>(m23Weights.size());
1211 
1212  //cout<<"nm13Weights = "<<nm13Weights<<", nm23Weights = "<<nm23Weights<<endl;
1213  // Print out abscissas and weights for the integration
1214  Double_t totm13Weight(0.0), totm23Weight(0.0);
1215  for (i = 0; i < nm13Weights; i++) {
1216  totm13Weight += m13Weights[i];
1217  }
1218  for (i = 0; i < nm23Weights; i++) {
1219  totm23Weight += m23Weights[i];
1220  }
1221  cout<<"totm13Weight = "<<totm13Weight<<", totm23Weight = "<<totm23Weight<<endl;
1222 
1223  std::vector<Double_t> m13(nm13Weights), m23(nm23Weights);
1224  std::vector<Double_t> m13Sq(nm13Weights), m23Sq(nm23Weights);
1225 
1226  // Use same number of abscissas for x and y co-ordinates
1227  Int_t m = (nm13Weights + 1)/2;
1228  for (i = 0; i < m; i++) {
1229 
1230  Int_t ii = nm13Weights - 1 - i; // symmetric i index
1231 
1232  Double_t dm13 = halfRangem13*m13Abscissas[i];
1233  Double_t m13Val = meanm13 - dm13;
1234  m13[i] = m13Val;
1235  m13Sq[i] = m13Val*m13Val;
1236 
1237  m13Val = meanm13 + dm13;
1238  m13[ii] = m13Val;
1239  m13Sq[ii] = m13Val*m13Val;
1240 
1241  }
1242 
1243  m = (nm23Weights +1)/2;
1244  for (i = 0; i < m; i++) {
1245 
1246  Int_t ii = nm23Weights - 1 - i; // symmetric i index
1247 
1248  Double_t dm23 = halfRangem23*m23Abscissas[i];
1249  Double_t m23Val = meanm23 - dm23;
1250  m23[i] = m23Val;
1251  m23Sq[i] = m23Val*m23Val;
1252  m23Val = meanm23 + dm23;
1253  m23[ii] = m23Val;
1254  m23Sq[ii] = m23Val*m23Val;
1255  }
1256 
1257  // Now compute the integral
1258  Double_t dpArea(0.0);
1259  for (i = 0; i < nm13Weights; i++) {
1260 
1261  for (j = 0; j < nm23Weights; j++) {
1262 
1263  Double_t weight = m13Weights[i]*m23Weights[j];
1264  Double_t Jacobian = 4.0*m13[i]*m23[j];
1265  weight *= (Jacobian*intFactor);
1266 
1267  // Calculate the integral contributions for each resonance.
1268  // Only resonances within the DP area contribute.
1269  // This also calculates the total DP area as a check.
1270  Bool_t withinDP = kinematics_->withinDPLimits(m13Sq[i], m23Sq[j]);
1271  if (withinDP == kTRUE) {
1272 
1273  kinematics_->updateKinematics(m13Sq[i], m23Sq[j]);
1274  this->dynamics(kFALSE, weight);
1275  // Increment total DP area
1276  dpArea += weight;
1277 
1278  }
1279 
1280  } // j weights loop
1281  } // i weights loop
1282 
1283  // Print out DP area to check whether we have a sensible output
1284  cout<<"dpArea = "<<dpArea<<endl;
1285 
1286 }
1287 
1288 void LauIsobarDynamics::dynamics(Bool_t cacheResData, Double_t weight, Bool_t useEff)
1289 {
1290  // Routine that calculates the Dalitz plot amplitude, incorporating
1291  // resonance dynamics (using resAmp()) and any interference between them.
1292  // Used by the fit() and sigGen() functions.
1293 
1294  UInt_t i(0), j(0);
1295 
1296  // Reset the total amplitude to zero
1297  totAmp_.zero();
1298 
1299  // Loop over the number of resonance amplitudes defined in the model
1300  // Have we already calculated this for this event (during fit?)
1301  // Or do we have a resonance that has varying pole mass/width/other factors?
1302  if (cacheResData == kFALSE) {
1303 
1304  for (i = 0; i < nAmp_; i++) {
1305 
1306  // Calculate the dynamics for this resonance, using the resAmp function.
1307  ff_[i] = resAmp(i);
1308 
1309  //cout<<"ff_["<<i<<"] = "<<ff_[i]<<endl;
1310 
1311  // If we have a symmetrical Dalitz plot, flip the m_13^2 and m_23^2
1312  // variables, recalculate the dynamics, and average both contributions.
1313  // Although, the factor of 0.5 cancels out with the fact that
1314  // the resonance appears on both sides of the Dalitz plot.
1315  if (symmetricalDP_ == kTRUE) { // was tieSg_ == 12
1317  ff_[i] += resAmp(i);
1318  // Flip the m_13^2 and m_23^2 variables back to their original values
1320  }
1321 
1322  } // Loop over amplitudes
1323 
1324  // If we haven't cached the data, then we need to find out the efficiency.
1325  eff_ = this->retrieveEfficiency();
1326 
1327  } // Already cached data?
1328 
1329  if (integralsDone_ == kTRUE) {
1330 
1331  // Loop over all signal amplitudes
1332  LauComplex ATerm;
1333  for (i = 0; i < nAmp_; i++) {
1334 
1335  // Get the partial complex amplitude - (mag, phase)*(resonance dynamics)
1336  ATerm = Amp_[i]*ff_[i];
1337  // Scale this contribution by its relative normalisation w.r.t. the whole dynamics
1338  ATerm.rescale(fNorm_[i]);
1339 
1340  // Add this partial amplitude to the sum
1341  //cout<<"For i = "<<i<<", ATerm = "<<ATerm<<", Amp = "<<Amp_[i]<<", ff = "<<ff_[i]<<endl;
1342  totAmp_ += ATerm;
1343 
1344  } // Loop over amplitudes
1345 
1346  // |Sum of partial amplitudes|^2
1347  ASq_ = totAmp_.abs2();
1348 
1349  // Apply the efficiency correction for this event.
1350  // Multiply the amplitude squared sum by the DP efficiency
1351  if ( useEff ) {
1352  ASq_ *= eff_;
1353  }
1354 
1355  } else { // integrals not done
1356 
1357  // Find the efficiency correction to be applied for this event.
1358  eff_ = this->retrieveEfficiency();
1359 
1360  Double_t effWeight = eff_*weight;
1361 
1362  // Need this for integrals for normalisation of likelihood function.
1363  LauComplex fifjEffSumTerm;
1364  LauComplex fifjSumTerm;
1365  for (i = 0; i < nAmp_; i++) {
1366 
1367  // Add the dynamical amplitude squared for this resonance.
1368  Double_t fSqVal = ff_[i].abs2();
1369  fSqSum_[i] += fSqVal*weight;
1370 
1371  for (j = i; j < nAmp_; j++) {
1372 
1373  fifjEffSumTerm = fifjSumTerm = ff_[i]*ff_[j].conj();
1374 
1375  fifjEffSumTerm.rescale(effWeight);
1376  fifjEffSum_[i][j] += fifjEffSumTerm;
1377 
1378  fifjSumTerm.rescale(weight);
1379  fifjSum_[i][j] += fifjSumTerm;
1380  }
1381  }
1382  }
1383 }
1384 
1386 {
1387  // Routine to calculate the resonance dynamics (amplitude)
1388  // using the appropriate Breit-Wigner/Form Factors.
1389  // Called by the dynamics() function.
1390 
1391  // Get the signal resonance from the stored vector
1392  LauAbsResonance* sigResonance = sigResonances_[index];
1393 
1394  if (sigResonance == 0) {
1395  cout<<"Error in resAmp. Couldn't retrieve resonance with index = "<<index<<endl;
1396  return LauComplex(0.0, 0.0);
1397  }
1398 
1399  // Get the integer index of the resonance.
1400  Int_t resInt = resIntAmp_[index];
1401 
1402  LauComplex resAmplitude(0.0, 0.0);
1403 
1404  // Fortran code has this set to 0 -> set NR index to -1
1405  if (resInt == -1) {
1406 
1407  // Non-resonant
1408  resAmplitude = LauComplex(0.3, 0.0);
1409 
1410  } else if (resInt < -1 || resInt >= static_cast<Int_t>(this->getnDefinedResonances())) {
1411 
1412  cout<<"Error in resAmp. Probably bad resonance name."<<endl;
1413  resAmplitude = LauComplex(0.0, 0.0);
1414 
1415  } else {
1416 
1417  resAmplitude = sigResonance->amplitude(kinematics_);
1418 
1419  }
1420 
1421  return resAmplitude;
1422 }
1423 
1424 void LauIsobarDynamics::setFFTerm(UInt_t index, Double_t realPart, Double_t imagPart)
1425 {
1426  // Function to set the internal ff = resAmp() term.
1427  if ( index >= nAmp_ ) {
1428  cout<<"Warning in setFFTerm. index = "<<index<<" is not within the range 0 and "
1429  <<nAmp_-1<<". Bailing out."<<endl;
1430  return;
1431  }
1432 
1433  ff_[index].setRealImagPart( realPart, imagPart );
1434 
1435 }
1436 
1438 {
1439  // This method calculates the fit fractions, mean efficiency and total DP rate
1440 
1441  Double_t fifjEffTot(0.0), fifjTot(0.0);
1442  UInt_t i, j;
1443  for (i = 0; i < nAmp_; i++) {
1444 
1445  // Calculate the diagonal terms
1446  TString name = "A"; name += i; name += "Sq_FitFrac";
1447  fitFrac_[i][i].name(name);
1448 
1449  Double_t fifjSumReal = fifjSum_[i][i].re();
1450  Double_t sumTerm = Amp_[i].abs2()*fifjSumReal*fNorm_[i]*fNorm_[i];
1451  fifjTot += sumTerm;
1452 
1453  Double_t fifjEffSumReal = fifjEffSum_[i][i].re();
1454  Double_t sumEffTerm = Amp_[i].abs2()*fifjEffSumReal*fNorm_[i]*fNorm_[i];
1455  fifjEffTot += sumEffTerm;
1456 
1457  fitFrac_[i][i] = sumTerm;
1458  }
1459 
1460  for (i = 0; i < nAmp_; i++) {
1461  for (j = i+1; j < nAmp_; j++) {
1462  // Calculate the cross-terms
1463  TString name = "A"; name += i; name += "A"; name += j; name += "_FitFrac";
1464  fitFrac_[i][j].name(name);
1465 
1466  LauComplex AmpjConj = Amp_[j].conj();
1467  LauComplex AmpTerm = Amp_[i]*AmpjConj;
1468 
1469  Double_t crossTerm = 2.0*(AmpTerm*fifjSum_[i][j]).re()*fNorm_[i]*fNorm_[j];
1470  fifjTot += crossTerm;
1471 
1472  Double_t crossEffTerm = 2.0*(AmpTerm*fifjEffSum_[i][j]).re()*fNorm_[i]*fNorm_[j];
1473  fifjEffTot += crossEffTerm;
1474 
1475  fitFrac_[i][j] = crossTerm;
1476  }
1477  }
1478 
1479  if (TMath::Abs(fifjTot) > 1e-10) {
1480  meanDPEff_ = fifjEffTot/fifjTot;
1481  if (init) {
1484  }
1485  }
1486  DPRate_ = fifjTot;
1487  if (init) {
1490  }
1491 
1492  // Now divide the fitFraction sums by the overall integral
1493  for (i = 0; i < nAmp_; i++) {
1494  for (j = i; j < nAmp_; j++) {
1495  // Get the actual fractions by dividing by the total DP rate
1496  fitFrac_[i][j] /= fifjTot;
1497  if (init) {
1498  fitFrac_[i][j].genValue( fitFrac_[i][j].value() );
1499  fitFrac_[i][j].initValue( fitFrac_[i][j].value() );
1500  }
1501  }
1502  }
1503 
1504  // Work out total fit fraction over all K-matrix components (for each propagator)
1505  KMPropMap::iterator mapIter;
1506  Int_t propInt(0);
1507 
1508  for (mapIter = kMatrixPropagators_.begin(); mapIter != kMatrixPropagators_.end(); ++mapIter) {
1509 
1510  LauKMatrixPropagator* thePropagator = mapIter->second;
1511 
1512  TString propName = thePropagator->getName();
1513 
1514  // Now loop over all resonances and find those which are K-matrix components for this propagator
1515  Double_t kMatrixTotFitFrac(0.0);
1516 
1517  for (i = 0; i < nAmp_; i++) {
1518 
1519  Bool_t gotKMRes1 = this->gotKMatrixMatch(i, propName);
1520  if (gotKMRes1 == kFALSE) {continue;}
1521 
1522  Double_t fifjSumReal = fifjSum_[i][i].re();
1523  Double_t sumTerm = Amp_[i].abs2()*fifjSumReal*fNorm_[i]*fNorm_[i];
1524 
1525  //Double_t fifjEffSumReal = fifjEffSum_[i][i].re();
1526  //Double_t sumEffTerm = Amp_[i].abs2()*fifjEffSumReal*fNorm_[i]*fNorm_[i];
1527 
1528  kMatrixTotFitFrac += sumTerm;
1529 
1530  for (j = i+1; j < nAmp_; j++) {
1531 
1532  Bool_t gotKMRes2 = this->gotKMatrixMatch(j, propName);
1533  if (gotKMRes2 == kFALSE) {continue;}
1534 
1535  LauComplex AmpjConj = Amp_[j].conj();
1536  LauComplex AmpTerm = Amp_[i]*AmpjConj;
1537 
1538  Double_t crossTerm = 2.0*(AmpTerm*fifjSum_[i][j]).re()*fNorm_[i]*fNorm_[j];
1539  //Double_t crossEffTerm = 2.0*(AmpTerm*fifjEffSum_[i][j]).re()*fNorm_[i]*fNorm_[j];
1540 
1541  kMatrixTotFitFrac += crossTerm;
1542 
1543  }
1544 
1545  }
1546 
1547  kMatrixTotFitFrac /= fifjTot;
1548 
1549  TString parName("KMatrixTotFF_"); parName += propInt;
1550  extraParameters_[propInt].name( parName );
1551  extraParameters_[propInt] = kMatrixTotFitFrac;
1552  if (init) {
1553  extraParameters_[propInt].genValue(kMatrixTotFitFrac);
1554  extraParameters_[propInt].initValue(kMatrixTotFitFrac);
1555  }
1556 
1557  cout<<"Total K-matrix fit fraction for propagator "<<propName<<" is "<<kMatrixTotFitFrac<<endl;
1558 
1559  ++propInt;
1560 
1561  }
1562 
1563 
1564 }
1565 
1566 Bool_t LauIsobarDynamics::gotKMatrixMatch(UInt_t resAmpInt, const TString& propName) const
1567 {
1568 
1569  Bool_t gotMatch(kFALSE);
1570 
1571  if (resAmpInt >= nAmp_) {return kFALSE;}
1572 
1573  const LauAbsResonance* theResonance = sigResonances_[resAmpInt];
1574 
1575  if (theResonance == 0) {return kFALSE;}
1576 
1577  Int_t resModelInt = theResonance->getResonanceModel();
1578 
1579  if (resModelInt == LauAbsResonance::KMatrix) {
1580 
1581  TString resName = theResonance->getResonanceName();
1582 
1583  KMStringMap::const_iterator kMPropSetIter = kMatrixPropSet_.find(resName);
1584 
1585  if (kMPropSetIter != kMatrixPropSet_.end()) {
1586  TString kmPropString = kMPropSetIter->second;
1587  if (kmPropString == propName) {gotMatch = kTRUE;}
1588  }
1589 
1590  }
1591 
1592  return gotMatch;
1593 
1594 }
1595 
1597 {
1598  // Calculate the normalisation for the log-likelihood function.
1599  DPNorm_ = 0.0;
1600 
1601  for (UInt_t i = 0; i < nAmp_; i++) {
1602  // fifjEffSum is the contribution from the term involving the resonance
1603  // dynamics (f_i for resonance i) and the efficiency term.
1604  Double_t fifjEffSumReal = fifjEffSum_[i][i].re();
1605  // We need to normalise this contribution w.r.t. the complete dynamics in the DP.
1606  // Hence we scale by the fNorm_i factor (squared), which is calculated by the
1607  // initialise() function, when the normalisation integrals are calculated and cached.
1608  // We also include the complex amplitude squared to get the total normalisation
1609  // contribution from this resonance.
1610  DPNorm_ += Amp_[i].abs2()*fifjEffSumReal*fNorm_[i]*fNorm_[i];
1611  }
1612 
1613  // We now come to the cross-terms (between resonances i and j) in the normalisation.
1614  for (UInt_t i = 0; i < nAmp_; i++) {
1615  for (UInt_t j = i+1; j < nAmp_; j++) {
1616  LauComplex AmpjConj = Amp_[j].conj();
1617  LauComplex AmpTerm = Amp_[i]*AmpjConj;
1618  // Again, fifjEffSum is the contribution from the term involving the resonance
1619  // dynamics (f_i*f_j_conjugate) and the efficiency cross term.
1620  // Also include the relative normalisation between these two resonances w.r.t. the
1621  // total DP dynamical structure (fNorm_i and fNorm_j) and the complex
1622  // amplitude squared (mag,phase) part.
1623  DPNorm_ += 2.0*(AmpTerm*fifjEffSum_[i][j]).re()*fNorm_[i]*fNorm_[j];
1624  }
1625  }
1626 
1627  return DPNorm_;
1628 }
1629 
1631 {
1632  // Routine to generate a signal event according to the Dalitz plot
1633  // model we have defined.
1634 
1635  nSigGenLoop_ = 0;
1636  Bool_t generatedSig(kFALSE);
1637 
1638  while (generatedSig == kFALSE && nSigGenLoop_ < iterationsMax_) {
1639 
1640  // Generates uniform DP phase-space distribution
1641  Double_t m13Sq(0.0), m23Sq(0.0);
1642  kinematics_->genFlatPhaseSpace(m13Sq, m23Sq);
1643 
1644  // If we're in a symmetrical DP then we should only generate events in one half
1645  if ( symmetricalDP_ && m13Sq > m23Sq ) {
1646  Double_t tmpSq = m13Sq;
1647  m13Sq = m23Sq;
1648  m23Sq = tmpSq;
1649  }
1650 
1651  // calculates the amplitudes and total amplitude for the given DP point
1652  this->calcLikelihoodInfo(m13Sq, m23Sq);
1653 
1654  if (integralsDone_ == kTRUE) {
1655 
1656  // No need for efficiency correction here. It is already done in
1657  // the dynamics() function (unlike the Fortran version).
1658 
1659  // Very important line to avoid bias in MC generation for the accept/reject method.
1660  // Make sure that the total amplitude squared is below some number, given
1661  // by aSqMaxSet_. If it is, then the event is valid.
1662  // Otherwise, go through another toy MC loop until we can generate the event
1663  // OK, or until we reach the maximum iteration limit.
1664 
1665  Double_t randNo = LauRandom::randomFun()->Rndm();
1666 
1667  if (randNo > ASq_/aSqMaxSet_) {
1668  nSigGenLoop_++;
1669  } else {
1670  generatedSig = kTRUE;
1671  nSigGenLoop_ = 0;
1672  if (ASq_ > aSqMaxVar_) {aSqMaxVar_ = ASq_;}
1673  }
1674 
1675  } else {
1676  // For toy MC numerical integration only
1677  generatedSig = kTRUE;
1678  }
1679  } // while loop
1680 
1681  Bool_t sigGenOK(kTRUE);
1682  if (GenOK != this->checkToyMC(kFALSE,kFALSE)) {
1683  sigGenOK = kFALSE;
1684  }
1685 
1686  return sigGenOK;
1687 }
1688 
1689 LauAbsDPDynamics::ToyMCStatus LauIsobarDynamics::checkToyMC(Bool_t printErrorMessages, Bool_t printInfoMessages)
1690 {
1691  // Check whether we have generated the toy MC OK.
1692  ToyMCStatus ok(GenOK);
1693 
1694  if (nSigGenLoop_ >= iterationsMax_) {
1695  if (printErrorMessages) {
1696  cerr<<"Error in sigGen. More than "<<iterationsMax_<<" iterations required."<<endl;
1697  cerr<<"Try to decrease the maximum allowed value of the total amplitude squared using the "
1698  <<"LauIsobarDynamics::setASqMaxValue(Double_t) function and re-run."<<endl;
1699  cerr<<"This needs to be, perhaps significantly, less than "<<aSqMaxSet_<<endl;
1700  cerr<<"Maximum value of ASq so far = "<<aSqMaxVar_<<endl;
1701  }
1702  aSqMaxSet_ = 1.01 * aSqMaxVar_;
1703  cout<<"|A|^2 max reset to "<<aSqMaxSet_<<endl;
1704  ok = MaxIterError;
1705  } else if (aSqMaxVar_ > aSqMaxSet_) {
1706  if (printErrorMessages) {
1707  cerr<<"Warning in LauIsobarDynamics::checkToyMC : aSqMaxSet_ was set to "<<aSqMaxSet_<<" but actual aSqMax was "<<aSqMaxVar_<<endl;
1708  cerr<<"Run was invalid, as any generated MC will be biased, according to the accept/reject method!"<<endl;
1709  cerr<<"Please set aSqMaxSet >= "<<aSqMaxVar_<<" using the LauIsobarDynamics::setASqMaxValue(Double_t) function and re-run."<<endl;
1710  }
1711  aSqMaxSet_ = 1.01 * aSqMaxVar_;
1712  cout<<"|A|^2 max reset to "<<aSqMaxSet_<<endl;
1713  ok = ASqMaxError;
1714  } else if (printInfoMessages) {
1715  cout<<"aSqMaxSet = "<<aSqMaxSet_<<" and aSqMaxVar = "<<aSqMaxVar_<<endl;
1716  }
1717 
1718  return ok;
1719 }
1720 
1722 {
1729  scfFraction_ = currentEvent_->retrieveScfFraction(); // These two are necessary, even though the dynamics don't actually use scfFraction_ or jacobian_,
1730  jacobian_ = currentEvent_->retrieveJacobian(); // since this is at the heart of the caching mechanism.
1731 }
1732 
1734 {
1735  // Calculate the likelihood and associated info
1736  // for the given event using cached information
1737  evtLike_ = 0.0;
1738 
1739  // retrieve the cached dynamics from the tree:
1740  // realAmp, imagAmp for each resonance plus efficiency, scf fraction and jacobian
1741  this->setDataEventNo(iEvt);
1742 
1743  // use realAmp and imagAmp to create the resonance amplitudes
1744  for (UInt_t i = 0; i < nAmp_; i++) {
1746  }
1747 
1748  // Update the dynamics - calculates totAmp_ and then ASq_ = totAmp_.abs2() * eff_
1749  // All calculated using cached information on the individual amplitudes and efficiency.
1750  this->dynamics(kTRUE, 1.0);
1751 
1752  // Calculate the normalised matrix element squared value
1753  if (DPNorm_ > 1e-10) {
1754  evtLike_ = ASq_/DPNorm_;
1755  }
1756 }
1757 
1758 void LauIsobarDynamics::calcLikelihoodInfo(Double_t m13Sq, Double_t m23Sq)
1759 {
1760  this->calcLikelihoodInfo(m13Sq, m23Sq, -1);
1761 }
1762 
1763 void LauIsobarDynamics::calcLikelihoodInfo(Double_t m13Sq, Double_t m23Sq, Int_t tagCat)
1764 {
1765  // Calculate the likelihood and associated info
1766  // for the given point in the Dalitz plot
1767  // Also retrieves the SCF fraction in the bin where the event lies (done
1768  // here to cache it along with the the rest of the DP quantities, like eff)
1769  // The jacobian for the square DP is calculated here for the same reason.
1770  evtLike_ = 0.0;
1771 
1772  // update the kinematics for the specified DP point
1773  kinematics_->updateKinematics(m13Sq, m23Sq);
1774 
1775  // calculate the jacobian and the scfFraction to cache them later
1776  scfFraction_ = this->retrieveScfFraction(tagCat);
1777  if (kinematics_->squareDP() == kTRUE) {
1778  // If cacheResData == kFALSE, updateKinematics has been called before dynamics(). Then get Jacobian.
1780  }
1781 
1782  // calculates the ff_ terms and retrieves eff_ from the efficiency model
1783  // then calculates totAmp_ and finally ASq_ = totAmp_.abs2() * eff_
1784  this->dynamics(kFALSE, 1.0);
1785 
1786  // Calculate the normalised matrix element squared value
1787  if (DPNorm_ > 1e-10) {
1788  evtLike_ = ASq_/DPNorm_;
1789  }
1790 }
1791 
1793 {
1794  // In LauFitDataTree, the first two variables should always be m13^2 and m23^2.
1795  // Other variables follow thus: charge/flavour tag prob, etc.
1796 
1797  UInt_t nBranches = inputFitTree.nBranches();
1798 
1799  if (nBranches < 2) {
1800  cout<<"Error in LauIsobarDynamics::fillDataTree. Expecting at least 2 variables "
1801  <<"in input data tree, but there are "<<nBranches<<"! Make sure you have "
1802  <<"the right number of variables in your input data file!"<<endl;
1803  gSystem->Exit(-1);
1804  }
1805 
1806  // Data structure that will cache the variables required to
1807  // calculate the signal likelihood for this experiment
1808  for ( std::vector<LauCacheData*>::iterator iter = data_.begin(); iter != data_.end(); ++iter ) {
1809  delete (*iter);
1810  }
1811  data_.clear();
1812 
1813  Double_t m13Sq(0.0), m23Sq(0.0);
1814  Double_t mPrime(0.0), thPrime(0.0);
1815  Int_t tagCat(-1);
1816  std::vector<Double_t> realAmp(nAmp_), imagAmp(nAmp_);
1817  Double_t eff(0.0), scfFraction(0.0), jacobian(0.0);
1818 
1819  UInt_t nEvents = inputFitTree.nEvents() + inputFitTree.nFakeEvents();
1820 
1821  data_.reserve(nEvents);
1822 
1823  for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) {
1824 
1825  const LauFitData& dataValues = inputFitTree.getData(iEvt);
1826  LauFitData::const_iterator iter = dataValues.find("m13Sq");
1827  m13Sq = iter->second;
1828  iter = dataValues.find("m23Sq");
1829  m23Sq = iter->second;
1830 
1831  // is there more than one tagging category?
1832  // if so then we need to know the category from the data
1833  if (scfFractionModel_.size()>1) {
1834  iter = dataValues.find("tagCat");
1835  tagCat = static_cast<Int_t>(iter->second);
1836  }
1837 
1838  // calculates the amplitudes and total amplitude for the given DP point
1839  // tagging category not needed by dynamics, but to find out the scfFraction
1840  this->calcLikelihoodInfo(m13Sq, m23Sq, tagCat);
1841 
1842  // extract the real and imaginary parts of the ff_ terms for storage
1843  for (UInt_t i = 0; i < nAmp_; i++) {
1844  realAmp[i] = ff_[i].re();
1845  imagAmp[i] = ff_[i].im();
1846  }
1847 
1848  if ( kinematics_->squareDP() ) {
1849  mPrime = kinematics_->getmPrime();
1850  thPrime = kinematics_->getThetaPrime();
1851  }
1852 
1853  eff = this->getEvtEff();
1854  scfFraction = this->getEvtScfFraction();
1855  jacobian = this->getEvtJacobian();
1856 
1857  // store the data for each event in the list
1858  data_.push_back( new LauCacheData() );
1859  data_[iEvt]->storem13Sq(m13Sq);
1860  data_[iEvt]->storem23Sq(m23Sq);
1861  data_[iEvt]->storemPrime(mPrime);
1862  data_[iEvt]->storethPrime(thPrime);
1863  data_[iEvt]->storeEff(eff);
1864  data_[iEvt]->storeScfFraction(scfFraction);
1865  data_[iEvt]->storeJacobian(jacobian);
1866  data_[iEvt]->storeRealAmp(realAmp);
1867  data_[iEvt]->storeImagAmp(imagAmp);
1868  }
1869 }
1870 
1872 {
1873  // Select the event (kinematics_) using an accept/reject method based on the
1874  // ratio of the current value of ASq to the maximal value.
1875  Bool_t accepted(kFALSE);
1876 
1877  this->dynamics(kFALSE, 1.0, kFALSE);
1878 
1879  // Compare the ASq value with the maximal value (set by the user)
1880  if (LauRandom::randomFun()->Rndm() < ASq_/aSqMaxSet_) {
1881  accepted = kTRUE;
1882  }
1883 
1884  if (ASq_ > aSqMaxVar_) {aSqMaxVar_ = ASq_;}
1885 
1886  return accepted;
1887 }
1888 
1890 {
1891  // calculate the dynamics from the current kinematics
1892  this->dynamics(kFALSE, 1.0, kFALSE);
1893 
1894  if (ASq_ > aSqMaxVar_) {aSqMaxVar_ = ASq_;}
1895 
1896  // return the event weight = the value of the squared amplitude divided
1897  // by the user-defined ceiling
1898  return ASq_ / aSqMaxSet_;
1899 }
1900 
KMStringMap kMatrixPropSet_
The names of the M-matrix components in the model mapped to their propagators.
Bool_t changeLASSEffectiveRange_
Whether the default value of the LASS effective range has been changed.
virtual Double_t getEvtJacobian() const
Retrieve the Jacobian, for the transformation into square DP coordinates, for the current event...
Int_t nSigGenLoop_
The number of unsucessful attempts to generate an event so far.
TString getNameDaug1() const
Get name of the first daughter particle.
Double_t getm23Max() const
Get the m23 maximum defined as (mParent - m1)
Class for defining the abstract interface for signal Dalitz plot dynamics.
Double_t retrievemPrime() const
Retrieve the square Dalitz plot coordinate, m&#39;.
virtual void initialise()=0
Initialise the model.
Bool_t squareDP() const
Are the square Dalitz plot co-ordinates being calculated?
TRandom * randomFun()
Access the singleton random number generator with a particular seed.
Definition: LauRandom.cc:20
File containing declaration of LauKMatrixPropFactory class.
File containing declaration of LauFitDataTree class.
File containing declaration of LauNRAmplitude class.
Bool_t changeLASSResonancePhase_
Whether the default value of the LASS resonance phase has been changed.
Double_t retrievem23Sq() const
Retrieve the invariant mass squared of the second and third daugthers.
Double_t resBarrierRadius_
The radius of the resonance barrier factor for new amplitude components.
virtual LauComplex resAmp(Int_t index)
Calculate the dynamic part of the amplitude for a given component at the current point in the Dalitz ...
virtual ~LauIsobarDynamics()
Destructor.
Int_t getTypeDaug1() const
Get PDG code of the first daughter particle.
const TString & getResonanceName() const
Get the name of the resonance.
void setBarrierRadii(Double_t resRadius, Double_t parRadius, LauAbsResonance::BarrierType type)
Set the form factor model and parameters.
Double_t FlatteParameterg1_
The constant parameter g1 for new Flatte components.
Double_t m13BinWidth_
The bin width to use when integrating over m13.
Bool_t flipHelicity_
The helicity flip flag for new amplitude components.
Double_t LASSBackgroundMag_
The background magnitude for new LASS components.
virtual void calcDPPartialIntegral(Double_t minm13, Double_t maxm13, Double_t minm23, Double_t maxm23, Double_t m13BinWidth, Double_t m23BinWidth)
Calculate the Dalitz plot normalisation integrals over a given range.
std::vector< Int_t > resPairAmp_
The index of the daughter not produced by the resonance for each amplitude component.
LauParameter DPRate_
The overall Dalitz plot rate.
virtual void addResonance(const TString &resName, Int_t resPairAmpInt, const TString &resType="RelBW", Double_t newMass=-1.0, Double_t newWidth=-1.0, Int_t newSpin=-1)
Add a resonance to the Dalitz plot.
File containing declaration of LauKMatrixProdPole class.
const TString & name() const
The parameter name.
Double_t getmPrime() const
Get m&#39; value.
Int_t getNPoles() const
Get the number of poles.
Class that defines the particular 3-body decay under study.
Definition: LauDaughters.hh:33
const std::vector< Double_t > & retrieveImagAmp() const
Retrieve the imaginary parts of the amplitudes.
Int_t getTypeParent() const
Get PDG code of the parent particle.
Int_t getCharge(Int_t resPairAmpInt) const
Get charge of a particular two-daughter combination.
File containing declaration of LauBelleNR class.
virtual void defineKMatrixPropagator(const TString &propName, const TString &paramFileName, Int_t resPairAmpInt, Int_t nChannels, Int_t nPoles, Int_t rowIndex=1)
Define a new K-matrix Propagator.
Class for defininf the Gounaris-Sakurai resonance model.
Double_t LASSScatteringLength_
The scattering length for new LASS components.
File containing declaration of LauKMatrixProdSVP class.
virtual void setAlpha(Double_t alpha)
Set the parameter alpha, the effective range.
Definition: LauBelleNR.hh:66
Bool_t changeLASSResonanceMag_
Whether the default value of the LASS resonance magnitude has been changed.
virtual void dynamics(Bool_t cacheResData=kTRUE, Double_t weight=1.0, Bool_t useEff=kTRUE)
Calculate the total Dalitz plot amplitude at the current point in the Dalitz plot.
File containing declaration of LauDaughters class.
Double_t aSqMaxVar_
The maximum value of A squared that has been seen so far while generating.
Double_t scfFraction_
The fraction of events that are poorly reconstructed (the self cross feed fraction) at the current po...
LauAbsResonance * getResonance(const TString &resName, Int_t resPairAmpInt, const TString &resType)
Create a resonance.
Double_t m13Sq_
The invariant mass squared of the first and third daughters.
Bool_t changeLASSBackgroundMag_
Whether the default value of the LASS background magnitude has been changed.
Int_t iterationsMax_
The maximum allowed number of attempts when generating an event.
LauComplex totAmp_
The total amplitude for the current event.
File containing declaration of LauGounarisSakuraiRes class.
LauAbsResonance::BarrierType barrierType_
The type of the barrier factor for new amplitude components.
Double_t aSqMaxSet_
The maximum allowed value of A squared.
File containing declaration of LauBelleSymNR class.
virtual Double_t retrieveEfficiency()
Obtain the efficiency of the current event from the model.
Bool_t symmetricalDP_
Whether the Dalitz plot is symmetrical.
File containing declaration of LauPrint class.
Int_t getNChannels() const
Get the number of channels.
Double_t getm13Min() const
Get the m13 minimum defined as (m1 + m3)
Double_t parBarrierRadius_
The radius of the parent barrier factor for new amplitude components.
Double_t retrieveJacobian() const
Retrieve the Jacobian for the transformation into square-Dalitz-plot coordinates. ...
Class to contain cached data relating to an event.
Definition: LauCacheData.hh:30
Class for performing numerical integration routines.
Definition: LauIntegrals.hh:30
void setBarrierRadii(Double_t resRadius, Double_t parRadius, LauAbsResonance::BarrierType type)
Set the form factor model and parameters.
std::map< TString, Double_t > LauFitData
Type for holding event data.
virtual void calcDPNormalisation()
Calculate the Dalitz plot normalisation integrals across the whole Dalitz plot.
TString getNameDaug2() const
Get name of the second daughter particle.
UInt_t nBranches() const
Obtain the number of branches in the tree.
virtual void setFFTerm(UInt_t index, Double_t realPart, Double_t imagPart)
Set the dynamic part of the amplitude for a given amplitude component at the current point in the Dal...
std::vector< LauComplex > Amp_
The complex coefficients for the amplitude components.
LauKMatrixPropagator * getPropagator(const TString &name, const TString &paramFileName, Int_t resPairAmpInt, Int_t nChannels, Int_t nPoles, Int_t rowIndex)
Retrieve the propagator if it already exists, otherwise create one.
KMPropMap kMatrixPropagators_
The K-matrix propagators.
std::vector< Double_t > fSqSum_
The event-by-event running total of the dynamical amplitude squared for each amplitude component...
Int_t resTypeInt(const TString &name) const
Retrieve the integer index for the specified resonance.
Bool_t flipHelicity() const
Get the helicity flip flag.
File containing declaration of LauKinematics class.
File containing declaration of LauRelBreitWignerRes class.
virtual Double_t getEvtEff() const
Retrieve the efficiency for the current event.
File containing declaration of LauKMatrixPropagator class.
Double_t getm13Max() const
Get the m13 maximum defined as (mParent - m2)
std::vector< LauCacheData * > data_
The cached data for all events.
Class for defining a K-matrix production pole amplitude term.
Double_t getm12Min() const
Get the m12 minimum defined as (m1 + m2)
Double_t retrievem13Sq() const
Retrieve the invariant mass squared of the first and third daugthers.
Double_t m23Sq_
The invariant mass squared of the second and third daughters.
Double_t thPrime_
The square Dalitz plot coordinate theta&#39;.
Double_t evtLike_
The normalised likelihood for the current event.
File containing declaration of LauIsobarDynamics class.
std::vector< std::vector< LauComplex > > fifjEffSum_
The event-by-event running total of efficiency corrected amplitude cross terms for each pair of ampli...
void updateKinematics(Double_t m13Sq, Double_t m23Sq)
Update all kinematic quantities based on the DP co-ordinates m13Sq and m23Sq.
virtual void initSummary()
Print a summary of the model to be used.
virtual void fillDataTree(const LauFitDataTree &fitDataTree)
Fill the internal data structure that caches the resonance dynamics.
virtual LauAbsResonance * findResonance(const TString &name)
Retrieve the named resonance.
Double_t ASq_
The value of A squared for the current event.
Double_t abs2() const
Obtain the square of the absolute value of the complex number.
Definition: LauComplex.hh:229
Bool_t gotSymmetricalDP() const
Is Dalitz plot symmetric.
Definition: LauDaughters.hh:66
virtual void changeResonance(const TString &resName, Double_t newMass=-1.0, Double_t newWidth=-1.0, Int_t newSpin=-1)
Change the properties of a resonance particle within this model.
std::vector< LauParameter > extraParameters_
any extra parameters/quantities (e.g. K-matrix total fit fractions)
static LauKMatrixPropFactory * getInstance()
Get a static instance of this factory class. Only one is created per application. ...
File containing declaration of LauResonanceMaker class.
Double_t m23BinWidth_
The bin width to use when integrating over m23.
void calcGaussLegendreWeights(const Int_t numPoints, std::vector< Double_t > &abscissas, std::vector< Double_t > &weights)
Calculate the Gauss-Legendre weights.
Definition: LauIntegrals.cc:37
TString intFileName_
The name of the file to save integrals to.
Class for defining the Belle nonresonant model.
Definition: LauBelleNR.hh:33
Int_t getResPairAmpInt() const
Get the DP axis identifier.
UInt_t nAmp_
The number of amplitude components.
virtual LauResonanceModel getResonanceModel() const =0
Get the resonance model type.
Double_t getm12Max() const
Get the m12 maximum defined as (mParent - m3)
virtual void addKMatrixProdPole(const TString &poleName, const TString &propName, Int_t poleIndex)
Add a K-matrix production pole term to the model.
Bool_t changeFlatteParameterg2_
Whether the default value of the Flatte parameter g2 has been changed.
std::vector< TString > resTypAmp_
The resonance types of all of the amplitude components.
virtual void calcLikelihoodInfo(UInt_t iEvt)
Calculate the likelihood (and all associated information) for the given event number.
const std::vector< Double_t > & retrieveRealAmp() const
Retrieve the real parts of the amplitudes.
std::vector< Int_t > resIntAmp_
The index within the resonance maker for each amplitude component.
LauParArray fitFrac_
The fit fractions for the amplitude components.
virtual void initialiseVectors()
Initialise the internal storage for this model.
Double_t jacobian_
The Jacobian, for the transformation into square DP coordinates at the current point in the Dalitz pl...
void flipAndUpdateKinematics()
Flips the DP variables m13^2 &lt;-&gt; m23^2 and recalculates all kinematic quantities. ...
Double_t LASSEffectiveRange_
The effective range for new LASS components.
std::vector< std::vector< LauComplex > > fifjSum_
The event-by-event running total of the amplitude cross terms for each pair of amplitude components...
std::map< Int_t, LauEffModel * > LauTagCatScfFractionModelMap
The type used for containing multiple self cross feed fraction models for different categories (e...
Double_t getm23Min() const
Get the m23 minimum defined as (m2 + m3)
UInt_t getnDefinedResonances() const
Retrieve the number of defined resonances in the resonance maker.
Double_t retrieveScfFraction() const
Retrieve the fraction of poorly reconstructed events (the so-called self cross feed fraction) ...
Double_t retrievethPrime() const
Retrieve the square Dalitz plot coordinate, theta&#39;.
virtual Double_t retrieveScfFraction(Int_t tagCat)
Obtain the self cross feed fraction of the current event from the model.
Bool_t changeLASSCutOff_
Whether the default value of the LASS cutoff has been changed.
virtual Double_t calcSigDPNorm()
Calculate the normalisation factor for the log-likelihood function.
virtual Double_t getEvtScfFraction() const
Retrieve the fraction of events that are poorly reconstructed (the self cross feed fraction) for the ...
Class for defining the relativistic Breit-Wigner resonance model.
Double_t FlatteParameterg2_
The constant parameter g2 for new Flatte components.
Bool_t withinDPLimits(Double_t m13Sq, Double_t m23Sq) const
Check whether a given (m13Sq,m23Sq) point is within the kinematic limits of the Dalitz plot...
Bool_t changeFlatteParameterg1_
Whether the default value of the Flatte parameter g1 has been changed.
Int_t getTypeDaug3() const
Get PDG code of the third daughter particle.
LauDaughters * daughters_
The daughters of the decay.
const LauFitData & getData(UInt_t iEvt) const
Retrieve the data for a given event.
File containing LauRandom namespace.
virtual Bool_t generate()
Generate a toy MC signal event.
File containing declaration of LauEffModel class.
virtual void removeCharge(TString &string) const
Remove the charge from the given particle name.
std::vector< Int_t > typDaug_
The PDG codes of the daughters.
Double_t eff_
The efficiency at the current point in the Dalitz plot.
Class for defining the symmetric Belle Non Resonant model.
Int_t getTypeDaug2() const
Get PDG code of the second daughter particle.
UInt_t nFakeEvents() const
Retrieve the number of fake events.
Double_t LASSBackgroundPhase_
The background phase for new LASS components.
Class that implements the efficiency description across the signal Dalitz plot.
Definition: LauEffModel.hh:37
virtual Bool_t gotReweightedEvent()
Calculates whether an event with the current kinematics should be accepted in order to produce a dist...
Abstract class for defining type for resonance amplitude models (Breit-Wigner, Flatte etc...
TString getNameParent() const
Get name of the parent particle.
Double_t LASSCutOff_
The cutoff value for new LASS components.
virtual void initialise()
Dummy initialisation.
virtual void updateCoeffs(const std::vector< LauComplex > &coeffs)
Update the complex coefficients for the resonances.
virtual void addKMatrixProdSVP(const TString &SVPName, const TString &propName, Int_t channelIndex)
Add a K-matrix slowly-varying part (SVP) term to the model.
void rescale(Double_t scaleVal)
Scale this by a factor.
Definition: LauComplex.hh:282
Int_t getChargeParent() const
Get charge of the parent particle.
void zero()
Set both real and imaginary part to zero.
Definition: LauComplex.hh:318
Bool_t changeLASSScatteringLength_
Whether the default value of the LASS scattering length has been changed.
Double_t DPNorm_
The normalisation factor for the log-likelihood function.
Double_t initValue() const
The initial value of the parameter.
Class for defining signal dynamics using the isobar model.
virtual void setDataEventNo(UInt_t iEvt)
Load the data for a given event.
virtual Double_t getEventWeight()
Calculate the acceptance rate, for events with the current kinematics, when generating events accordi...
File containing LauConstants namespace.
ToyMCStatus
The possible statuses for toy MC generation.
virtual void writeIntegralsFile()
Write the results of the integrals (and related information) to a file.
virtual ToyMCStatus checkToyMC(Bool_t printErrorMessages=kTRUE, Bool_t printInfoMessages=kFALSE)
Check the status of the toy MC generation.
File containing declaration of LauAbsResonance class.
virtual void setIntegralBinWidths(Double_t m13BinWidth, Double_t m23BinWidth)
Set the widths of the bins to use when integrating across the Dalitz plot.
LauResonanceMaker * resonanceMaker_
Object to create resonances.
LauKinematics * kinematics_
The kinematics of the decay.
Double_t getThetaPrime() const
Get theta&#39; value.
Class for defining a complex number.
Definition: LauComplex.hh:47
LauIsobarDynamics(LauDaughters *daughters, LauEffModel *effModel, LauEffModel *scfFractionModel=0)
Constructor.
Bool_t changeLASSBackgroundPhase_
Whether the default value of the LASS background phase has been changed.
Bool_t gotKMatrixMatch(UInt_t resAmpInt, const TString &propName) const
Check whether a resonance is a K-matrix component of a given propagator.
void genFlatPhaseSpace(Double_t &m13Sq, Double_t &m23Sq) const
Routine to generate events flat in phase-space.
File containing declaration of LauIntegrals class.
TString getName() const
Get the propagator name.
virtual void setDataEventNo(UInt_t iEvt)
Load the data for a given event.
LauCacheData * currentEvent_
The cached data for the current event.
Double_t value() const
The value of the parameter.
UInt_t nEvents() const
Retrieve the number of events.
LauParameter meanDPEff_
The mean efficiency across the Dalitz plot.
void changeResonance(Double_t newMass, Double_t newWidth, Int_t newSpin)
Allow the mass, width and spin of the resonance to be changed.
Double_t calcSqDPJacobian()
Calculate the Jacobian for the transformation m23^2, m13^2 -&gt; m&#39;, theta&#39; (square DP) ...
Bool_t integralsDone_
Whether the integrals have been performed.
virtual LauComplex amplitude(const LauKinematics *kinematics)
Calculate the complex amplitude.
std::vector< LauComplex > ff_
The dynamic part of the amplitude for each amplitude component at the current point in the Dalitz plo...
Double_t LASSResonancePhase_
The resonance phase for new LASS components.
Double_t LASSResonanceMag_
The resonance magnitude for new LASS components.
virtual void initialise(const std::vector< LauComplex > &coeffs)
Initialise the Dalitz plot dynamics.
TString getNameDaug3() const
Get name of the third daughter particle.
Double_t genValue() const
The value generated for the parameter.
Double_t mPrime_
The square Dalitz plot coordinate, m&#39;.
std::vector< Double_t > fNorm_
The normalisation factors for the dynamic parts of the amplitude for each amplitude component...
Class for defining a K-matrix production &quot;slowly-varying part&quot; (SVP) amplitude.
virtual void setResonanceParameter(Double_t value, const TString &name)
Set the updated parameters from changeResonance.
Class to store the input fit variables.
Double_t retrieveEff() const
Retrieve the efficiency.
Double_t BelleNRAlpha_
The alpha parameter for new Belle non-resonant components.
std::vector< LauAbsResonance * > sigResonances_
The resonances in the model.
virtual void calcExtraInfo(Bool_t init=kFALSE)
Calculate the fit fractions, mean efficiency and total DP rate.
LauTagCatScfFractionModelMap scfFractionModel_
The self cross feed fraction models across the Dalitz plot.
Class for defining a K-matrix propagator.