MINT2
BsDsKpipi_TD_ampFit.cpp
Go to the documentation of this file.
1 // Toy studies with full TD amplitude PDF
2 // author: Philippe d'Argent
3 #include "Mint/FitParameter.h"
4 #include "Mint/NamedParameter.h"
5 #include "Mint/Minimiser.h"
6 #include "Mint/Neg2LL.h"
7 #include "Mint/Neg2LLSum.h"
8 #include "Mint/DalitzEventList.h"
9 #include "Mint/NamedDecayTreeList.h"
10 #include "Mint/DecayTree.h"
11 #include "Mint/DiskResidentEventList.h"
12 #include "Mint/CLHEPPhysicalConstants.h"
13 #include "Mint/CLHEPSystemOfUnits.h"
14 #include "Mint/PdfBase.h"
15 #include "Mint/DalitzPdfBase.h"
16 #include "Mint/DalitzPdfBaseFastInteg.h"
17 #include "Mint/DalitzPdfBaseFlexiFastInteg.h"
18 #include "Mint/FitAmplitude.h"
19 #include "Mint/FitAmpSum.h"
20 #include "Mint/FitAmpIncoherentSum.h"
21 #include "Mint/DalitzEvent.h"
22 #include "Mint/AmpRatios.h"
23 #include "Mint/IEventGenerator.h"
24 #include "Mint/DalitzBWBoxSet.h"
25 #include "Mint/DalitzBoxSet.h"
26 #include "Mint/SignalGenerator.h"
27 #include "Mint/FromFileGenerator.h"
28 #include "Mint/DalitzSumPdf.h"
29 #include "Mint/cexp.h"
30 #include "Mint/DalitzPdfNormChecker.h"
31 #include "Mint/IFastAmplitudeIntegrable.h"
32 #include "Mint/DalitzPdfSaveInteg.h"
33 #include "Mint/Chi2Binning.h"
34 #include "Mint/FitAmpIncoherentSum.h"
35 #include "Mint/FitAmpList.h"
36 #include "Mint/DalitzPdfBaseMCInteg.h"
37 #include "TGraph.h"
38 #include "TFile.h"
39 #include "TCanvas.h"
40 #include "TNtupleD.h"
41 #include "TTree.h"
42 #include "TFile.h"
43 #include "TLegend.h"
44 #include <TStyle.h>
45 #include "TRandom2.h"
46 #include "TRandom3.h"
47 #include <ctime>
48 #include <iostream>
49 #include <fstream>
50 #include <algorithm>
51 //#include <omp.h>
52 #include <TROOT.h>
53 #include "RooRealConstant.h"
54 #include "RooAbsReal.h"
55 #include "RooRealVar.h"
56 #include "RooDataSet.h"
57 #include "RooGlobalFunc.h"
58 #include "RooRealVar.h"
59 #include "RooDecay.h"
60 #include "RooBDecay.h"
61 #include "RooPlot.h"
62 #include "RooEffProd.h"
63 #include "RooGenericPdf.h"
64 #include "RooGaussModel.h"
65 #include "RooProdPdf.h"
66 #include "RooConstVar.h"
67 #include "RooCategory.h"
68 #ifndef __CINT__
69 #include "RooGlobalFunc.h"
70 #endif
71 #include "Mint/RooCubicSplineFun.h"
72 #include "Mint/RooCubicSplineKnot.h"
73 #include "Mint/RooGaussEfficiencyModel.h"
74 #include "Mint/DecRateCoeff_Bd.h"
75 
76 using namespace std;
77 using namespace RooFit ;
78 using namespace MINT;
79 
80 class AmpsPdfFlexiFast
81 : public DalitzPdfBaseFlexiFastInteg
82 {
83 protected:
84  TRandom* _localRnd;
85  SignalGenerator* _sgGen;
86  FromFileGenerator* _fileGen;
87  IEventGenerator<IDalitzEvent>* _chosenGen;
88  NamedParameter<std::string> _integratorSource;
89  std::string _integratorFileName;
90 public:
91  double un_normalised_noPs(IDalitzEvent& evt){
92  double ampSq = _amps->RealVal(evt);
93  return ampSq;// * getEvent()->phaseSpace();
94  }
95 
96  std::complex<double> ComplexVal_un_normalised_noPs(IDalitzEvent& evt){
97  return _amps->ComplexVal(evt);
98  }
99 
100  AmpsPdfFlexiFast(const DalitzEventPattern& pat
101  , IFastAmplitudeIntegrable* amps
102  , MinuitParameterSet* mps
103  , double precision=1.e-4
104  , std::string method="efficient"
105  , std::string fname = "SignalIntegrationEvents.root", bool genMoreEvents = false
106  )
107  : DalitzPdfBaseFlexiFastInteg(pat, 0, amps, precision, mps)
108  , _localRnd(0)
109  , _sgGen(0)
110  , _fileGen(0)
111  , _chosenGen(0)
112  , _integratorSource("IntegratorSource", (std::string) "new", (char*) 0)
113  , _integratorFileName(fname)
114  {
115  cout << " AmpsPdfFlexiFast with integ method " << method << endl;
116  bool nonFlat = "efficient" == method;
117  bool generateNew = ((string)_integratorSource == (string)"new");
118  if(nonFlat){
119  cout << "AmpsPdfFlexiFast uses nonFlat integration." << endl;
120  if(generateNew){
121  _sgGen = new SignalGenerator(pat);
122  _sgGen->setWeighted();
123  //_sgGen->setSaveEvents(_integratorFileName);
124  _chosenGen = _sgGen;
125  }else{
126  // here, SignalGenerator is used by FromFileGenerator, to fill
127  // up missing events in case more are needed than found in the
128  // file. Since we don't know with which random seed the
129  // events in the file were generated, we supply a random
130  // number generator with randomised seed.
131  _localRnd = new TRandom3(time(0));
132  _sgGen = new SignalGenerator(pat, _localRnd);
133  _sgGen->setWeighted();
134  _sgGen->dontSaveEvents();// saving events is done by FromFileGenerator
135  if(genMoreEvents) _fileGen = new FromFileGenerator(_integratorFileName, _sgGen);
136  else{
137  _fileGen = new FromFileGenerator(_integratorFileName, 0, "UPDATE");
138  cout << "not going to generate any more events" << endl;
139  }
140  _chosenGen = _fileGen;
141  }
142  this->setEventGenerator(_chosenGen);
143  }else{
144  cout << "AmpsPdfFlexiFast uses flat integration." << endl;
145  }
146  }
147 
148  IFastAmplitudeIntegrable* getAmpSum(){ return _amps;}
149 
150  ~AmpsPdfFlexiFast(){
151  if(0 != _fileGen) delete _fileGen;
152  if(0 != _sgGen) delete _sgGen;
153  if(0 != _localRnd) delete _localRnd;
154  }
155 };
156 
157 class AmpRatio : virtual public IReturnComplex{
158  FitParameter& _re;
159  FitParameter& _im;
160  int _f;
161 public:
162  AmpRatio(FitParameter& re, FitParameter& im,int f)
163  : _re(re), _im(im), _f(f) {}
164 
165  complex<double> ComplexVal(){
166  std::complex<double> result(_re,static_cast<double>(_f) * _im);
167  return result;
168  }
169 };
170 
171 class CPV_amp : virtual public IReturnComplex{
172  FitParameter& _re;
173  FitParameter& _im;
174  int _sign;
175 public:
176  CPV_amp(FitParameter& re, FitParameter& im, int sign)
177  : _re(re), _im(im), _sign(sign) {}
178 
179  complex<double> ComplexVal(){
180  std::complex<double> result((double) ( 1.+ _re * (double) _sign),(double) (_im * (double) _sign) );
181  return result;
182  }
183 };
184 
185 class CPV_amp_polar : virtual public IReturnComplex{
186  FitParameter& _r;
187  FitParameter& _delta;
188  int _sign;
189 public:
190  CPV_amp_polar(FitParameter& r, FitParameter& delta, int sign)
191  : _r(r), _delta(delta), _sign(sign) {}
192 
193  complex<double> ComplexVal(){
194  std::complex<double> result= polar((double) sqrt( 1.+ _r * (double) _sign),(double) (_delta/360.*2.*pi * (double) _sign) );
195  return result;
196  }
197 };
198 
199 void AddScaledAmpsToList(FitAmpSum& fas_tmp, FitAmpSum& fas, FitAmpSum& fasCC, std::string name, counted_ptr<IReturnComplex>& r_plus, counted_ptr<IReturnComplex>& r_minus){
200 
201  counted_ptr<FitAmpList> List = fas_tmp.GetCloneOfSubsetSameFitParameters(name);
202  FitAmpSum fas_2(*List);
203  FitAmpSum fasCC_2(*List);
204  fasCC_2.CPConjugateSameFitParameters();
205  fasCC_2.CConjugateFinalStateSameFitParameters();
206  fas_2.multiply(r_plus);
207  fasCC_2.multiply(r_minus);
208  fas.addAsList(fas_2,1.);
209  fasCC.addAsList(fasCC_2,1.);
210 }
211 
212 std::vector<double> coherenceFactor(FitAmpSum& fas, FitAmpSum& fas_bar, double r, double delta, DalitzEventList& eventList){
213 
214  cout << "Calculating coherence factor ..." << endl << endl;
215  //fas.print();
216  //fas_bar.print();
217 
218  std::complex<double> valK(0,0);
219  double val1 = 0;
220  double val2 = 0;
221 
222  const complex<double> phase_diff = polar(r, delta/360.*2*pi);
223 
224  for(unsigned int i=0; i<eventList.size(); i++){
225  const std::complex<double> amp = fas.getVal(eventList[i]) ;
226  const std::complex<double> amp_bar = fas_bar.getVal(eventList[i])*phase_diff ;
227  valK += amp_bar*conj(amp)*eventList[i].getWeight()/ eventList[i].getGeneratorPdfRelativeToPhaseSpace();
228  val1 += norm(amp)*eventList[i].getWeight()/ eventList[i].getGeneratorPdfRelativeToPhaseSpace();
229  val2 += norm(amp_bar)*eventList[i].getWeight()/ eventList[i].getGeneratorPdfRelativeToPhaseSpace();
230  }
231 
232  std::vector<double> result;
233  result.push_back(sqrt(val2/val1));
234  result.push_back(std::abs(valK)/sqrt(val1)/sqrt(val2));
235  result.push_back(std::arg(valK)/(2.*pi)*360.);
236 
237  cout << "r = " << result[0] << endl;
238  cout << "k = " << result[1] << endl;
239  cout << "d = " << result[2] << " [deg]" << endl << endl;
240 
241  return result;
242 }
243 
244 enum basisType { noBasis=0 , expBasis= 3
245  , sinBasis=13, cosBasis=23
246  , sinhBasis=63, coshBasis=53 };
247 
248 // Full time-dependent PDF
249 class FullAmpsPdfFlexiFastCPV : public MINT::PdfBase<IDalitzEvent>
250 , virtual public IDalitzPdf{
251 
252 protected:
253  AmpsPdfFlexiFast* _amps1;
254  AmpsPdfFlexiFast* _amps2;
255  AmpsPdfFlexiFast* _ampsSum;
256 
257  FitParameter& _r;
258  FitParameter& _delta;
259  FitParameter& _gamma;
260 
261  FitParameter& _tau;
262  FitParameter& _dGamma;
263  FitParameter& _dm;
264  FitParameter& _eff_tag;
265  FitParameter& _w;
266 
267  // cast of MINT parameters to B2DX parameters
268  RooRealVar* _r_t;
269  RooRealVar* _r_dt;
270  RooCategory* _r_q;
271  RooRealVar* _r_mistag;
272  RooCategory* _r_f;
273 
274  RooRealVar* _r_scale_mean_dt;
275  RooRealVar* _r_scale_sigma_dt;
276  RooAbsPdf* _pdf_sigma_t;
277  RooAbsPdf* _pdf_omega_os;
278  RooAbsPdf* _pdf_omega_ss;
279 
280  RooCubicSplineFun* _spline;
281  RooGaussEfficiencyModel* _efficiency;
282  NamedParameter<double> _min_TAU;
283  NamedParameter<double> _max_TAU;
284 
285  DecRateCoeff_Bd* _cos_coeff;
286  DecRateCoeff_Bd* _cosh_coeff;
287  DecRateCoeff_Bd* _sin_coeff;
288  DecRateCoeff_Bd* _sinh_coeff;
289 
290  double _intA;
291  double _intAbar;
292  complex<double> _intAAbar;
293 
294 public:
295  void parametersChanged(){
296  _ampsSum->parametersChanged();
297  _intA = _ampsSum->integralForMatchingPatterns(true,1);
298  _intAbar = _ampsSum->integralForMatchingPatterns(true,-1);
299  _intAAbar = _ampsSum->ComplexSumForMatchingPatterns(false);
300  }
301  void beginFit(){
302  _ampsSum->beginFit();
303  printIntegralVals();
304  }
305  void endFit(){
306  printIntegralVals();
307  _ampsSum->endFit();
308  }
309 
310  void printIntegralVals(){
311  cout << "intSum = " << _ampsSum->getIntegralValue() << endl;
312  cout << "intA = " << _intA << endl;
313  cout << "intAbar = " << _intAbar << endl;
314  cout << "intAAbar = " << _intAAbar.real() << endl;
315  }
316 
317  inline double getValForGeneration(IDalitzEvent& evt){
318  const double t = (double) evt.getValueFromVector(0);
319  const double dt = (double) evt.getValueFromVector(1);
320  const double q = static_cast<double>((int)evt.getValueFromVector(2));
321  const double w = (double) evt.getValueFromVector(3);
322  const double f = static_cast<double>((int)evt.getValueFromVector(4));
323 
324  if(t < _min_TAU || t > _max_TAU )return 0.;
325  _r_t->setVal(t);
326  _r_dt->setVal(dt);
327  _r_q->setIndex(q);
328  _r_mistag->setVal(w);
329  _r_f->setIndex(f);
330 
331  double r = (double)_r; // * sqrt(_intA/_intAbar);
332  const complex<double> phase_diff = polar((double)r,((double) _delta -(double)_gamma*f)/360.*2*pi);
333 
334  const std::complex<double> amp = _amps1->ComplexVal_un_normalised_noPs(evt) ;
335  const std::complex<double> amp_bar = _amps2->ComplexVal_un_normalised_noPs(evt) * phase_diff;
336 
337  const double val = exp(-fabs(t)/(double)_tau) *
338  (
339  (norm(amp) + norm(amp_bar)) *cosh((double)_dGamma/2.*t) * _cosh_coeff->evaluate()
340  +(norm(amp) - norm(amp_bar)) *cos((double)_dm*t) * _cos_coeff->evaluate()
341  +real(amp_bar*conj(amp)) *sinh((double)_dGamma/2.*t) * _sinh_coeff->evaluate()
342  +imag(amp_bar*conj(amp)) *sin((double)_dm*t) * _sin_coeff->evaluate()
343  ) *_pdf_sigma_t->getVal()*_spline->getVal();
344 
345  return val;
346  }
347 
348  inline void getSmearedTime(double& t, double& dt, TRandom3& r){
349  while(true){
350  t = r.Exp(_tau);
351  dt = r.Uniform(0.,0.25);
352  t = t + r.Gaus(0,_r_scale_sigma_dt->getVal()*dt);
353  if(_min_TAU< t && t < _max_TAU) return;
354  }
355  }
356 
357  inline double un_normalised_noPs(IDalitzEvent& evt){
358  const double t = (double) evt.getValueFromVector(0);
359  const double dt = (double) evt.getValueFromVector(1);
360  const double q = static_cast<double>((int)evt.getValueFromVector(2));
361  const double w = (double) evt.getValueFromVector(3);
362  const double f = static_cast<double>((int)evt.getValueFromVector(4));
363 
364  if(t < _min_TAU || t > _max_TAU )return 0.;
365  _r_t->setVal(t);
366  _r_dt->setVal(dt);
367  _r_q->setIndex(q);
368  _r_mistag->setVal(w);
369  _r_f->setIndex(f);
370 
371  double r = (double)_r; // * sqrt(_intA/_intAbar);
372  const complex<double> phase_diff = polar((double)r,((double) _delta -(double)_gamma*f)/360.*2*pi);
373 
374  const std::complex<double> amp = _amps1->ComplexVal_un_normalised_noPs(evt) ;
375  const std::complex<double> amp_bar = _amps2->ComplexVal_un_normalised_noPs(evt) * phase_diff;
376 
377  const double cosh_term = _efficiency->evaluate(coshBasis,_tau,_dm,_dGamma);
378  const double cos_term = _efficiency->evaluate(cosBasis,_tau,_dm,_dGamma);
379  const double sinh_term = _efficiency->evaluate(sinhBasis,_tau,_dm,_dGamma);
380  const double sin_term = _efficiency->evaluate(sinBasis,_tau,_dm,_dGamma);
381 
382  const double val =
383  (
384  (norm(amp) + norm(amp_bar)) *cosh_term * _cosh_coeff->evaluate()
385  + (norm(amp) - norm(amp_bar)) *cos_term * _cos_coeff->evaluate()
386  + real(amp_bar*conj(amp)) *sinh_term * _sinh_coeff->evaluate()
387  + imag(amp_bar*conj(amp)) *sin_term * _sin_coeff->evaluate()
388  )* _pdf_sigma_t->getVal();
389 
390  return val;
391  }
392 
393  virtual double getVal(IDalitzEvent& evt){
394 
395  const double val = un_normalised_noPs(evt);
396  double r = (double)_r; // * sqrt(_intA/_intAbar);
397  const double f = static_cast<double>((int)evt.getValueFromVector(4));
398 
399  const complex<double> phase_diff_m = polar((double)r,((double) _delta + (double)_gamma)/360.*2*pi);
400  const complex<double> int_interference_m = phase_diff_m * _intAAbar ;
401 
402  const complex<double> phase_diff_p = polar((double)r,((double) _delta -(double)_gamma)/360.*2*pi);
403  const complex<double> int_interference_p = phase_diff_p * _intAAbar ;
404 
405  const complex<double> phase_diff = polar((double)r,((double) _delta -(double)_gamma*f)/360.*2*pi);
406  const complex<double> int_interference = phase_diff * _intAAbar ;
407 
408  TString range = "";
409  //if(f == 1) range = "Range_p";
410  //else range = "Range_m";
411 
412 
413  if(_intA == -1 ){
414  cout << "AmpsPdfFlexiFastCPV:: _norm = -1, should not have happened." << endl;
415  throw "can't deal with that";
416  }
417 
418  double norm = (_intA + r* r * _intAbar) * _efficiency->analyticalIntegral(coshBasis,_tau,_dm,_dGamma) * _cosh_coeff->analyticalIntegral(1,range)
419  + (_intA - r* r * _intAbar) * _efficiency->analyticalIntegral(cosBasis,_tau,_dm,_dGamma) * _cos_coeff->analyticalIntegral(1,range)
420  //+ (int_interference_m.real() * _sinh_coeff->analyticalIntegral(1,"Range_m") + int_interference_p.real() * _sinh_coeff->analyticalIntegral(1,"Range_p")) * _efficiency->analyticalIntegral(sinhBasis,_tau,_dm,_dGamma)
421  //+ (int_interference_m.imag() * _sin_coeff->analyticalIntegral(1,"Range_m") + int_interference_p.imag() * _sin_coeff->analyticalIntegral(1,"Range_p")) * _efficiency->analyticalIntegral(sinBasis,_tau,_dm,_dGamma) ;
422  + int_interference.real()/2. * _sinh_coeff->analyticalIntegral(1,range) * _efficiency->analyticalIntegral(sinhBasis,_tau,_dm,_dGamma)
423  + int_interference.imag()/2. * _sin_coeff->analyticalIntegral(1,range) * _efficiency->analyticalIntegral(sinBasis,_tau,_dm,_dGamma) ;
424 
425  /*
426  const double t = (double) evt.getValueFromVector(0);
427  const double dt = (double) evt.getValueFromVector(1);
428  const double q = static_cast<double>((int)evt.getValueFromVector(2));
429  const double w = (double) evt.getValueFromVector(3);
430  cout << t << endl;
431  cout << dt << endl;
432  cout << q << endl;
433  cout << w << endl;
434  cout << f << endl << endl;
435 
436  cout << _cosh_coeff->evaluate() << endl;
437  cout << _cosh_coeff->analyticalIntegral(1) << endl;
438  cout << _cosh_coeff->analyticalIntegral(1,"Range_p") << endl;
439  cout << _cosh_coeff->analyticalIntegral(1,"Range_m") << endl << endl;
440 
441  cout << _cos_coeff->evaluate() << endl;
442  cout << _cos_coeff->analyticalIntegral(1) << endl;
443  cout << _cos_coeff->analyticalIntegral(1,"Range_p") << endl;
444  cout << _cos_coeff->analyticalIntegral(1,"Range_m") << endl << endl;
445 
446  cout << _sinh_coeff->evaluate() << endl;
447  cout << _sinh_coeff->analyticalIntegral(1) << endl;
448  cout << _sinh_coeff->analyticalIntegral(1,"Range_p") << endl;
449  cout << _sinh_coeff->analyticalIntegral(1,"Range_m") << endl << endl;
450 
451  cout << _sin_coeff->evaluate() << endl;
452  cout << _sin_coeff->analyticalIntegral(1) << endl;
453  cout << _sin_coeff->analyticalIntegral(1,"Range_p") << endl;
454  cout << _sin_coeff->analyticalIntegral(1,"Range_m") << endl << endl;
455 
456  throw "";
457  */
458  return val/norm*2.;
459  }
460 
461  virtual double getVal_withPs(IDalitzEvent& evt){return getVal(evt);}
462  virtual double getVal_noPs(IDalitzEvent& evt){return getVal(evt);}
463 
464  virtual double getVal(IDalitzEvent* evt){
465  if(0 == evt) return 0;
466  return getVal(*evt);
467  }
468  virtual double getVal_withPs(IDalitzEvent* evt){
469  if(0 == evt) return 0;
470  return getVal_withPs(*evt);
471  }
472  virtual double getVal_noPs(IDalitzEvent* evt){
473  if(0 == evt) return 0;
474  return getVal_noPs(*evt);
475  }
476 
477  virtual DalitzHistoSet histoSet(){return _ampsSum->histoSet();}
478 
479  void doFinalStatsAndSaveForAmp12(MINT::Minimiser* min=0,const std::string& fname = "FitAmpResults", const std::string& fnameROOT="fitFractions"){
480  _amps1->redoIntegrator();
481  _amps2->redoIntegrator();
482  _amps1->doFinalStatsAndSave(min,((string)fname+".txt").c_str(),((string)fnameROOT+".root").c_str());
483  _amps2->doFinalStatsAndSave(min,((string)fname+"_CC.txt").c_str(),((string)fnameROOT+"_CC.root").c_str());
484  }
485 
486  FullAmpsPdfFlexiFastCPV(AmpsPdfFlexiFast* amps1, AmpsPdfFlexiFast* amps2, AmpsPdfFlexiFast* ampsSum,
487  MINT::FitParameter& r, MINT::FitParameter& delta,MINT::FitParameter& gamma,
488  MINT::FitParameter& tau, MINT::FitParameter& dGamma, MINT::FitParameter& dm, MINT::FitParameter& eff_tag, MINT::FitParameter& w ):
489  _amps1(amps1),_amps2(amps2),_ampsSum(ampsSum),_r(r),_delta(delta),_gamma(gamma),_tau(tau),_dGamma(dGamma),_dm(dm),_eff_tag(eff_tag),_w(w),
490  _intA(-1),_intAbar(-1),_intAAbar(-1),_min_TAU("min_TAU", 0.4), _max_TAU("max_TAU", 10.)
491  {
492  // Init B2DX parameters
493  _r_t = new RooRealVar("t", "time", _min_TAU, _max_TAU);
494  _r_dt = new RooRealVar("dt", "per-candidate time resolution estimate",0., 0.25);
495  _r_mistag = new RooRealVar("mistag", "mistag",0.);
496  _r_f = new RooCategory("qf", "qf");
497  _r_f->defineType("h+", +1);
498  _r_f->defineType("h-", -1);
499  _r_f->setRange("Range_p","h+");
500  _r_f->setRange("Range_m","h-");
501  _r_q = new RooCategory("qt", "qt");
502  _r_q->defineType("B+", +1);
503  _r_q->defineType("B-", -1) ;
504  _r_q->defineType("untagged", 0);
505 
506  _cosh_coeff = new DecRateCoeff_Bd("cosh_coeff",
507  "cosh_coeff",
508  DecRateCoeff_Bd::kCosh ,
509  *_r_f,
510  RooRealConstant::value(1.),
511  RooRealConstant::value(1.),
512  *_r_q,
513  *_r_mistag,
514  RooRealConstant::value(0.),
515  RooRealConstant::value(1.),
516  RooRealConstant::value(0.),
517  RooRealConstant::value(0.),
518  RooRealConstant::value(0.),
519  RooRealConstant::value(_eff_tag),
520  RooRealConstant::value(0.),
521  RooRealConstant::value(0.),
522  RooRealConstant::value(0.));
523 
524  _cos_coeff = new DecRateCoeff_Bd("cos_coeff",
525  "cos_coeff",
526  DecRateCoeff_Bd::kCos ,
527  *_r_f,
528  RooRealConstant::value(1.),
529  RooRealConstant::value(-1.),
530  *_r_q,
531  *_r_mistag,
532  RooRealConstant::value(0.),
533  RooRealConstant::value(1.),
534  RooRealConstant::value(0.),
535  RooRealConstant::value(0.),
536  RooRealConstant::value(0.),
537  RooRealConstant::value(_eff_tag),
538  RooRealConstant::value(0.),
539  RooRealConstant::value(0.),
540  RooRealConstant::value(0.));
541 
542  _sinh_coeff = new DecRateCoeff_Bd("sinh_coeff",
543  "sinh_coeff",
544  DecRateCoeff_Bd::kSinh ,
545  *_r_f,
546  RooRealConstant::value(-2.),
547  RooRealConstant::value(-2.),
548  *_r_q,
549  *_r_mistag,
550  RooRealConstant::value(0.),
551  RooRealConstant::value(1.),
552  RooRealConstant::value(0.),
553  RooRealConstant::value(0.),
554  RooRealConstant::value(0.),
555  RooRealConstant::value(_eff_tag),
556  RooRealConstant::value(0.),
557  RooRealConstant::value(0.),
558  RooRealConstant::value(0.));
559 
560  _sin_coeff = new DecRateCoeff_Bd("sin_coeff",
561  "sin_coeff",
562  DecRateCoeff_Bd::kSin,
563  *_r_f,
564  RooRealConstant::value(2.),
565  RooRealConstant::value(-2.),
566  *_r_q,
567  *_r_mistag,
568  RooRealConstant::value(0.),
569  RooRealConstant::value(1.),
570  RooRealConstant::value(0.),
571  RooRealConstant::value(0.),
572  RooRealConstant::value(0.),
573  RooRealConstant::value(_eff_tag),
574  RooRealConstant::value(0.),
575  RooRealConstant::value(0.),
576  RooRealConstant::value(0.));
577 
578  _r_scale_mean_dt = new RooRealVar("scale_mean_dt", "scale_mean_dt", 0);
579  _r_scale_sigma_dt = new RooRealVar("scale_sigma_dt", "scale_sigma_dt", 1.2);
580 
581  //SPLINE KNOTS
582  NamedParameter<double> knot_positions("knot_positions", 0.5, 1., 1.5, 2., 3., 6., 9.5);
583  vector<double> myBinning = knot_positions.getVector();
584 
585  NamedParameter<double> knot_values("knot_values", 0.38,0.63,0.86,1.05,1.14,1.24,1.22);
586  vector<double> values = knot_values.getVector() ;
587 
588  //SPLINE COEFFICIENTS
589  RooArgList tacc_list;
590  for(int i= 0; i< values.size(); i++){
591  tacc_list.add(*(new RooRealVar(("coeff_"+anythingToString(i)).c_str(), ("coeff_"+anythingToString(i)).c_str(), values[i], 0.0, 10.0)));
592  }
593  tacc_list.add(*(new RooRealVar(("coeff_"+anythingToString((int)values.size())).c_str(), ("coeff_"+anythingToString((int)values.size())).c_str(), 1.0)));
594 
595  RooFormulaVar* coeff_last = new RooFormulaVar(("coeff_"+anythingToString((int)values.size()+1)).c_str(),("coeff_"+anythingToString((int)values.size()+1)).c_str(), "@0 + ((@0-@1)/(@2-@3)) * (@4 - @2)", RooArgList(RooRealConstant::value(1.0), *tacc_list.find(("coeff_"+anythingToString((int)values.size()-1)).c_str()) , RooRealConstant::value(myBinning[myBinning.size()-1]), RooRealConstant::value(myBinning[myBinning.size()-2]), RooRealConstant::value(_r_t->getMax()) ));
596 
597  tacc_list.add(*coeff_last);
598 
599  //CUBIC SPLINE FUNCTION
600  _spline = new RooCubicSplineFun("splinePdf", "splinePdf", *_r_t, myBinning, tacc_list);
601  _efficiency = new RooGaussEfficiencyModel("resmodel", "resmodel", *_r_t, *_spline, RooRealConstant::value(0.), *_r_dt, *_r_scale_mean_dt, *_r_scale_sigma_dt );
602 
603  // Plot acceptance
604  TH1F *h_spline = new TH1F("", "", 100, _min_TAU, _max_TAU);
605  for (int i = 1; i<=h_spline->GetNbinsX(); i++) {
606  _r_t->setVal(h_spline->GetXaxis()->GetBinCenter(i));
607  h_spline->SetBinContent(i,_spline->getVal());
608  }
609 
610  TCanvas* c = new TCanvas();
611  h_spline->SetLineColor(kRed);
612  h_spline->Draw("histc");
613  c->Print("spline.eps");
614  c->Print("spline.pdf");
615 
616  // Init pdf for sigma_t
617  _pdf_sigma_t = new RooGenericPdf("pdf_sigma_t","pow(7. / @1, 7) / 720. * pow(@0, 6) * exp(-7. * @0 / @1)",RooArgList(*_r_dt, RooRealConstant::value(0.04)));
618 
619  }
620 };
621 
622 class AmpsPdfFlexiFastCPV_mod : public MINT::PdfBase<IDalitzEvent>
623 , virtual public IDalitzPdf{
624 
625 protected:
626  AmpsPdfFlexiFast* _amps1;
627  AmpsPdfFlexiFast* _amps2;
628  AmpsPdfFlexiFast* _ampsSum;
629 
630  FitParameter& _r;
631  FitParameter& _delta;
632  FitParameter& _gamma;
633 
634  FitParameter& _tau;
635  FitParameter& _dGamma;
636  FitParameter& _dm;
637  FitParameter& _eff_tag;
638  FitParameter& _w;
639 
640  double _intA;
641  double _intAbar;
642  complex<double> _intAAbar;
643 
644  // cast of MINT parameters to B2DX parameters
645  RooRealVar* _r_t;
646  RooRealVar* _r_dt;
647  RooCategory* _r_q;
648  RooRealVar* _r_mistag;
649  RooCategory* _r_f;
650 
651  RooRealVar* _r_scale_mean_dt;
652  RooRealVar* _r_scale_sigma_dt;
653  RooAbsPdf* _pdf_sigma_t;
654  RooAbsPdf* _pdf_omega_os;
655  RooAbsPdf* _pdf_omega_ss;
656 
657  RooCubicSplineFun* _spline;
658  RooGaussEfficiencyModel* _efficiency;
659  NamedParameter<double> _min_TAU;
660  NamedParameter<double> _max_TAU;
661 
662  DecRateCoeff_Bd* _cos_coeff;
663  DecRateCoeff_Bd* _cosh_coeff;
664  DecRateCoeff_Bd* _sin_coeff;
665  DecRateCoeff_Bd* _sinh_coeff;
666 
667 public:
668  void parametersChanged(){
669  _ampsSum->parametersChanged();
670  _intA = _ampsSum->integralForMatchingPatterns(true,1);
671  _intAbar = _ampsSum->integralForMatchingPatterns(true,-1);
672  _intAAbar = _ampsSum->ComplexSumForMatchingPatterns(false);
673  }
674  void beginFit(){
675  _ampsSum->beginFit();
676  printIntegralVals();
677  }
678  void endFit(){
679  printIntegralVals();
680  _ampsSum->endFit();
681  }
682 
683  void printIntegralVals(){
684  cout << "intSum = " << _ampsSum->getIntegralValue() << endl;
685  cout << "intA = " << _intA << endl;
686  cout << "intAbar = " << _intAbar << endl;
687  cout << "intAAbar = " << _intAAbar.real() << endl;
688  }
689 
690  inline double un_normalised_noPs(IDalitzEvent& evt){
691  const double t = (double) evt.getValueFromVector(0);
692  const double dt = (double) evt.getValueFromVector(1);
693  const double q = static_cast<double>((int)evt.getValueFromVector(2));
694  const double w = (double) evt.getValueFromVector(3);
695  const double f = static_cast<double>((int)evt.getValueFromVector(4));
696  _r_t->setVal(t);
697  _r_dt->setVal(dt);
698  _r_q->setIndex(q);
699  _r_mistag->setVal(w);
700  _r_f->setIndex(f);
701 
702  //const double e_eff = fabs(q)*_eff_tag + (1.-fabs(q))*(1.-_eff_tag);
703  double r = (double)_r; // * sqrt(_intA/_intAbar);
704  const complex<double> phase_diff = polar((double)r,((double) _delta -(double)_gamma*f)/360.*2*pi);
705 
706  const std::complex<double> amp = _amps1->ComplexVal_un_normalised_noPs(evt) ;
707  const std::complex<double> amp_bar = _amps2->ComplexVal_un_normalised_noPs(evt) * phase_diff;
708 
709  const double val = exp(-fabs(t)/(double)_tau) *
710  (
711  (norm(amp) + norm(amp_bar))*cosh((double)_dGamma/2.*t)* _cosh_coeff->evaluate()
712  +(norm(amp) - norm(amp_bar)) *cos((double)_dm*t) * _cos_coeff->evaluate()
713  +real(amp_bar*conj(amp))*sinh((double)_dGamma/2.*t) * _sinh_coeff->evaluate()
714  +imag(amp_bar*conj(amp))*sin((double)_dm*t) * _sin_coeff->evaluate()
715  );
716 
717  return val;
718  }
719 
720  virtual double getVal(IDalitzEvent& evt){
721 
722  const double f = static_cast<double>((int)evt.getValueFromVector(4));
723  const double val = un_normalised_noPs(evt);
724 
725  double r = (double)_r; // * sqrt(_intA/_intAbar);
726  const double Gamma = 1./((double) _tau);
727  const complex<double> phase_diff = polar((double)r,((double) _delta -(double)_gamma*f)/360.*2*pi);
728  const double int_interference = (phase_diff*_intAAbar).real();
729 
730  const complex<double> phase_diff_m = polar((double)r,((double) _delta + (double)_gamma)/360.*2*pi);
731  const complex<double> int_interference_m = phase_diff_m * _intAAbar ;
732 
733  const complex<double> phase_diff_p = polar((double)r,((double) _delta -(double)_gamma)/360.*2*pi);
734  const complex<double> int_interference_p = phase_diff_p * _intAAbar ;
735 
736 
737  if(_intA == -1 ){
738  cout << "AmpsPdfFlexiFastCPV:: _norm = -1, should not have happened." << endl;
739  throw "can't deal with that";
740  }
741 
742  //return val/(( _cosh_coeff->analyticalIntegral(1) * (_intA + r* r * _intAbar) * 4.*Gamma + _sinh_coeff->analyticalIntegral(1) * (int_interference_m.real() + int_interference_p.real())/4. * 2. * _dGamma )/ (4.*Gamma*Gamma-_dGamma*_dGamma));
743 
744  return val/(( _cosh_coeff->analyticalIntegral(1) * (_intA + r* r * _intAbar) * 4.*Gamma + _sinh_coeff->analyticalIntegral(1) * int_interference/2. * 2. * _dGamma )/ (4.*Gamma*Gamma-_dGamma*_dGamma))*2.;
745  }
746 
747  virtual double getVal_withPs(IDalitzEvent& evt){return getVal(evt);}
748  virtual double getVal_noPs(IDalitzEvent& evt){return getVal(evt);}
749 
750  virtual double getVal(IDalitzEvent* evt){
751  if(0 == evt) return 0;
752  return getVal(*evt);
753  }
754  virtual double getVal_withPs(IDalitzEvent* evt){
755  if(0 == evt) return 0;
756  return getVal_withPs(*evt);
757  }
758  virtual double getVal_noPs(IDalitzEvent* evt){
759  if(0 == evt) return 0;
760  return getVal_noPs(*evt);
761  }
762 
763  virtual DalitzHistoSet histoSet(){return _ampsSum->histoSet();}
764 
765  void doFinalStatsAndSaveForAmp12(MINT::Minimiser* min=0,const std::string& fname = "FitAmpResults", const std::string& fnameROOT="fitFractions"){
766  _amps1->redoIntegrator();
767  _amps2->redoIntegrator();
768  _amps1->doFinalStatsAndSave(min,((string)fname+".txt").c_str(),((string)fnameROOT+".root").c_str());
769  _amps2->doFinalStatsAndSave(min,((string)fname+"_CC.txt").c_str(),((string)fnameROOT+"_CC.root").c_str());
770  }
771 
772  AmpsPdfFlexiFastCPV_mod(AmpsPdfFlexiFast* amps1, AmpsPdfFlexiFast* amps2, AmpsPdfFlexiFast* ampsSum,
773  MINT::FitParameter& r, MINT::FitParameter& delta,MINT::FitParameter& gamma,
774  MINT::FitParameter& tau, MINT::FitParameter& dGamma, MINT::FitParameter& dm, MINT::FitParameter& eff_tag, MINT::FitParameter& w ):
775  _amps1(amps1),_amps2(amps2),_ampsSum(ampsSum),_r(r),_delta(delta),_gamma(gamma),_tau(tau),_dGamma(dGamma),_dm(dm),_eff_tag(eff_tag),_w(w),
776  _intA(-1),_intAbar(-1),_intAAbar(-1),_min_TAU("min_TAU", 0.4), _max_TAU("max_TAU", 10.)
777  {
778  // Init B2DX parameters
779  _r_t = new RooRealVar("t", "time", _min_TAU, _max_TAU);
780  _r_dt = new RooRealVar("dt", "per-candidate time resolution estimate",0., 0.25);
781  _r_mistag = new RooRealVar("mistag", "mistag",0.);
782  _r_f = new RooCategory("qf", "qf");
783  _r_f->defineType("h+", +1);
784  _r_f->defineType("h-", -1);
785  _r_f->setRange("Range_p","h+");
786  _r_f->setRange("Range_m","h-");
787  _r_q = new RooCategory("qt", "qt");
788  _r_q->defineType("B+", +1);
789  _r_q->defineType("B-", -1) ;
790  _r_q->defineType("untagged", 0);
791 
792  _cosh_coeff = new DecRateCoeff_Bd("cosh_coeff",
793  "cosh_coeff",
794  DecRateCoeff_Bd::kCosh ,
795  *_r_f,
796  RooRealConstant::value(1.),
797  RooRealConstant::value(1.),
798  *_r_q,
799  *_r_mistag,
800  RooRealConstant::value(0.),
801  RooRealConstant::value(1.),
802  RooRealConstant::value(0.),
803  RooRealConstant::value(0.),
804  RooRealConstant::value(0.),
805  RooRealConstant::value(_eff_tag),
806  RooRealConstant::value(0.),
807  RooRealConstant::value(0.),
808  RooRealConstant::value(0.));
809 
810  _cos_coeff = new DecRateCoeff_Bd("cos_coeff",
811  "cos_coeff",
812  DecRateCoeff_Bd::kCos ,
813  *_r_f,
814  RooRealConstant::value(1.),
815  RooRealConstant::value(-1.),
816  *_r_q,
817  *_r_mistag,
818  RooRealConstant::value(0.),
819  RooRealConstant::value(1.),
820  RooRealConstant::value(0.),
821  RooRealConstant::value(0.),
822  RooRealConstant::value(0.),
823  RooRealConstant::value(_eff_tag),
824  RooRealConstant::value(0.),
825  RooRealConstant::value(0.),
826  RooRealConstant::value(0.));
827 
828  _sinh_coeff = new DecRateCoeff_Bd("sinh_coeff",
829  "sinh_coeff",
830  DecRateCoeff_Bd::kSinh ,
831  *_r_f,
832  RooRealConstant::value(-2.),
833  RooRealConstant::value(-2.),
834  *_r_q,
835  *_r_mistag,
836  RooRealConstant::value(0.),
837  RooRealConstant::value(1.),
838  RooRealConstant::value(0.),
839  RooRealConstant::value(0.),
840  RooRealConstant::value(0.),
841  RooRealConstant::value(_eff_tag),
842  RooRealConstant::value(0.),
843  RooRealConstant::value(0.),
844  RooRealConstant::value(0.));
845 
846  _sin_coeff = new DecRateCoeff_Bd("sin_coeff",
847  "sin_coeff",
848  DecRateCoeff_Bd::kSin,
849  *_r_f,
850  RooRealConstant::value(2.),
851  RooRealConstant::value(-2.),
852  *_r_q,
853  *_r_mistag,
854  RooRealConstant::value(0.),
855  RooRealConstant::value(1.),
856  RooRealConstant::value(0.),
857  RooRealConstant::value(0.),
858  RooRealConstant::value(0.),
859  RooRealConstant::value(_eff_tag),
860  RooRealConstant::value(0.),
861  RooRealConstant::value(0.),
862  RooRealConstant::value(0.));
863  // Init pdf for sigma_t
864  _pdf_sigma_t = new RooGenericPdf("pdf_sigma_t","pow(7. / @1, 7) / 720. * pow(@0, 6) * exp(-7. * @0 / @1)",RooArgList(*_r_dt, RooRealConstant::value(0.04)));
865 
866  }
867 
868 };
869 
870 // Time-dependent PDF
871 class AmpsPdfFlexiFastCPV : public MINT::PdfBase<IDalitzEvent>
872 , virtual public IDalitzPdf{
873 
874 protected:
875  AmpsPdfFlexiFast* _amps1;
876  AmpsPdfFlexiFast* _amps2;
877  AmpsPdfFlexiFast* _ampsSum;
878 
879  FitParameter& _r;
880  FitParameter& _delta;
881  FitParameter& _gamma;
882 
883  FitParameter& _tau;
884  FitParameter& _dGamma;
885  FitParameter& _dm;
886  FitParameter& _eff_tag;
887  FitParameter& _w;
888 
889  double _intA;
890  double _intAbar;
891  complex<double> _intAAbar;
892 
893 public:
894  void parametersChanged(){
895  _ampsSum->parametersChanged();
896  _intA = _ampsSum->integralForMatchingPatterns(true,1);
897  _intAbar = _ampsSum->integralForMatchingPatterns(true,-1);
898  _intAAbar = _ampsSum->ComplexSumForMatchingPatterns(false);
899  }
900  void beginFit(){
901  _ampsSum->beginFit();
902  printIntegralVals();
903  }
904  void endFit(){
905  printIntegralVals();
906  _ampsSum->endFit();
907  }
908 
909  void printIntegralVals(){
910  cout << "intSum = " << _ampsSum->getIntegralValue() << endl;
911  cout << "intA = " << _intA << endl;
912  cout << "intAbar = " << _intAbar << endl;
913  cout << "intAAbar = " << _intAAbar.real() << endl;
914  }
915 
916  inline double un_normalised_noPs(IDalitzEvent& evt){
917  const double t = (double) evt.getValueFromVector(0);
918  const double dt = (double) evt.getValueFromVector(1);
919  const double q = static_cast<double>((int)evt.getValueFromVector(2));
920  const double w = (double) evt.getValueFromVector(3);
921  const double f = static_cast<double>((int)evt.getValueFromVector(4));
922 
923  const double e_eff = fabs(q)*_eff_tag + (1.-fabs(q))*(1.-_eff_tag);
924 
925  double r = (double)_r; // * sqrt(_intA/_intAbar);
926  const complex<double> phase_diff = polar((double)r,((double) _delta -(double)_gamma*f)/360.*2*pi);
927 
928  const std::complex<double> amp = _amps1->ComplexVal_un_normalised_noPs(evt) ;
929  const std::complex<double> amp_bar = _amps2->ComplexVal_un_normalised_noPs(evt) * phase_diff;
930 
931  const double val = exp(-fabs(t)/(double)_tau) *
932  (
933  (2.-fabs(q))*(norm(amp) + norm(amp_bar))*cosh((double)_dGamma/2.*t)
934  +f*q*(1.-2.*w)*(norm(amp) - norm(amp_bar)) *cos((double)_dm*t)
935  -(2.-fabs(q))*2.0*real(amp_bar*conj(amp))*sinh((double)_dGamma/2.*t)
936  -f*2.0*q*(1.-2.*w)*imag(amp_bar*conj(amp))*sin((double)_dm*t)
937  )*e_eff;
938 
939  return val;
940  }
941 
942  virtual double getVal(IDalitzEvent& evt){
943 
944  const double f = static_cast<double>((int)evt.getValueFromVector(4));
945  const double val = un_normalised_noPs(evt);
946 
947  double r = (double)_r; // * sqrt(_intA/_intAbar);
948  const double Gamma = 1./((double) _tau);
949  const complex<double> phase_diff = polar((double)r,((double) _delta -(double)_gamma*f)/360.*2*pi);
950  const double int_interference = (phase_diff*_intAAbar).real();
951 
952  if(_intA == -1 ){
953  cout << "AmpsPdfFlexiFastCPV:: _norm = -1, should not have happened." << endl;
954  throw "can't deal with that";
955  }
956 
957  return val/(2.* ((_intA + r* r * _intAbar) * 4.*Gamma - int_interference * 2. * _dGamma )/ (4.*Gamma*Gamma-_dGamma*_dGamma));
958  }
959 
960  virtual double getVal_withPs(IDalitzEvent& evt){return getVal(evt);}
961  virtual double getVal_noPs(IDalitzEvent& evt){return getVal(evt);}
962 
963  virtual double getVal(IDalitzEvent* evt){
964  if(0 == evt) return 0;
965  return getVal(*evt);
966  }
967  virtual double getVal_withPs(IDalitzEvent* evt){
968  if(0 == evt) return 0;
969  return getVal_withPs(*evt);
970  }
971  virtual double getVal_noPs(IDalitzEvent* evt){
972  if(0 == evt) return 0;
973  return getVal_noPs(*evt);
974  }
975 
976  virtual DalitzHistoSet histoSet(){return _ampsSum->histoSet();}
977 
978  void doFinalStatsAndSaveForAmp12(MINT::Minimiser* min=0,const std::string& fname = "FitAmpResults", const std::string& fnameROOT="fitFractions"){
979  _amps1->redoIntegrator();
980  _amps2->redoIntegrator();
981  _amps1->doFinalStatsAndSave(min,((string)fname+".txt").c_str(),((string)fnameROOT+".root").c_str());
982  _amps2->doFinalStatsAndSave(min,((string)fname+"_CC.txt").c_str(),((string)fnameROOT+"_CC.root").c_str());
983  }
984 
985  AmpsPdfFlexiFastCPV(AmpsPdfFlexiFast* amps1, AmpsPdfFlexiFast* amps2, AmpsPdfFlexiFast* ampsSum,
986  MINT::FitParameter& r, MINT::FitParameter& delta,MINT::FitParameter& gamma,
987  MINT::FitParameter& tau, MINT::FitParameter& dGamma, MINT::FitParameter& dm, MINT::FitParameter& eff_tag, MINT::FitParameter& w ):
988  _amps1(amps1),_amps2(amps2),_ampsSum(ampsSum),_r(r),_delta(delta),_gamma(gamma),_tau(tau),_dGamma(dGamma),_dm(dm),_eff_tag(eff_tag),_w(w),
989  _intA(-1),_intAbar(-1),_intAAbar(-1)
990  {;}
991 };
992 
993 class AmpsPdfFlexiFastCPV_time_integrated : public MINT::PdfBase<IDalitzEvent>
994 , virtual public IDalitzPdf{
995 
996 protected:
997  AmpsPdfFlexiFast* _amps1;
998  AmpsPdfFlexiFast* _amps2;
999  AmpsPdfFlexiFast* _ampsSum;
1000 
1001  FitParameter& _r;
1002  FitParameter& _delta;
1003  FitParameter& _gamma;
1004 
1005  FitParameter& _tau;
1006  FitParameter& _dGamma;
1007  FitParameter& _dm;
1008  FitParameter& _eff_tag;
1009  FitParameter& _w;
1010 
1011  double _intA;
1012  double _intAbar;
1013  complex<double> _intAAbar;
1014 
1015 public:
1016  void parametersChanged(){
1017  _ampsSum->parametersChanged();
1018  _intA = _ampsSum->integralForMatchingPatterns(true,1);
1019  _intAbar = _ampsSum->integralForMatchingPatterns(true,-1);
1020  _intAAbar = _ampsSum->ComplexSumForMatchingPatterns(false);
1021  }
1022  void beginFit(){
1023  _ampsSum->beginFit();
1024  printIntegralVals();
1025  }
1026  void endFit(){
1027  printIntegralVals();
1028  _ampsSum->endFit();
1029  }
1030 
1031  void printIntegralVals(){
1032  cout << "intSum = " << _ampsSum->getIntegralValue() << endl;
1033  cout << "intA = " << _intA << endl;
1034  cout << "intAbar = " << _intAbar << endl;
1035  cout << "intAAbar = " << _intAAbar.real() << endl;
1036  }
1037 
1038  inline double un_normalised_noPs(IDalitzEvent& evt){
1039  const double t = (double) evt.getValueFromVector(0);
1040  const double dt = (double) evt.getValueFromVector(1);
1041  const double q = static_cast<double>((int)evt.getValueFromVector(2));
1042  const double w = (double) evt.getValueFromVector(3);
1043  const double f = static_cast<double>((int)evt.getValueFromVector(4));
1044 
1045  const double e_eff = fabs(q)*_eff_tag + (1.-fabs(q))*(1.-_eff_tag);
1046 
1047  double r = (double)_r; // * sqrt(_intA/_intAbar);
1048  const complex<double> phase_diff = polar((double)r,((double) _delta -(double)_gamma*f)/360.*2*pi);
1049 
1050  const std::complex<double> amp = _amps1->ComplexVal_un_normalised_noPs(evt) ;
1051  const std::complex<double> amp_bar = _amps2->ComplexVal_un_normalised_noPs(evt) * phase_diff;
1052 
1053  double Gamma = 1./((double) _tau);
1054 
1055  const double val =
1056  (
1057  (2.-fabs(q))*(norm(amp) + norm(amp_bar))*(4.*Gamma/(4.*Gamma*Gamma-_dGamma*_dGamma))
1058  +f*q*(1.-2.*w)*(norm(amp) - norm(amp_bar)) *(Gamma/(Gamma*Gamma+_dm*_dm))
1059  -(2.-fabs(q))*2.0*real(amp_bar*conj(amp))*(2.*_dGamma/(4.*Gamma*Gamma-_dGamma*_dGamma))
1060  -f*2.0*q*(1.-2.*w)*imag(amp_bar*conj(amp))*(_dm/(Gamma*Gamma+_dm*_dm))
1061  )*e_eff;
1062 
1063  return val;
1064  }
1065 
1066  virtual double getVal(IDalitzEvent& evt){
1067 
1068  const double f = static_cast<double>((int)evt.getValueFromVector(4));
1069  const double val = un_normalised_noPs(evt);
1070 
1071  double r = (double)_r; // * sqrt(_intA/_intAbar);
1072  const double Gamma = 1./((double) _tau);
1073  const complex<double> phase_diff = polar((double)r,((double) _delta -(double)_gamma*f)/360.*2*pi);
1074  const double int_interference = (phase_diff*_intAAbar).real();
1075 
1076  if(_intA == -1 ){
1077  cout << "AmpsPdfFlexiFastCPV:: _norm = -1, should not have happened." << endl;
1078  throw "can't deal with that";
1079  }
1080 
1081  return val/(2.* ((_intA + r* r * _intAbar) * 4.*Gamma - int_interference * 2. * _dGamma )/ (4.*Gamma*Gamma-_dGamma*_dGamma));
1082  }
1083 
1084  virtual double getVal_withPs(IDalitzEvent& evt){return getVal(evt);}
1085  virtual double getVal_noPs(IDalitzEvent& evt){return getVal(evt);}
1086 
1087  virtual double getVal(IDalitzEvent* evt){
1088  if(0 == evt) return 0;
1089  return getVal(*evt);
1090  }
1091  virtual double getVal_withPs(IDalitzEvent* evt){
1092  if(0 == evt) return 0;
1093  return getVal_withPs(*evt);
1094  }
1095  virtual double getVal_noPs(IDalitzEvent* evt){
1096  if(0 == evt) return 0;
1097  return getVal_noPs(*evt);
1098  }
1099 
1100  virtual DalitzHistoSet histoSet(){return _ampsSum->histoSet();}
1101 
1102 
1103  void doFinalStatsAndSaveForAmp12(MINT::Minimiser* min=0,const std::string& fname = "FitAmpResults", const std::string& fnameROOT="fitFractions"){
1104  _amps1->redoIntegrator();
1105  _amps2->redoIntegrator();
1106  _amps1->doFinalStatsAndSave(min,((string)fname+".txt").c_str(),((string)fnameROOT+"_TI.root").c_str());
1107  _amps2->doFinalStatsAndSave(min,((string)fname+"_CC.txt").c_str(),((string)fnameROOT+"_TI_CC.root").c_str());
1108  }
1109 
1110  AmpsPdfFlexiFastCPV_time_integrated(AmpsPdfFlexiFast* amps1, AmpsPdfFlexiFast* amps2, AmpsPdfFlexiFast* ampsSum,
1111  MINT::FitParameter& r, MINT::FitParameter& delta,MINT::FitParameter& gamma,
1112  MINT::FitParameter& tau, MINT::FitParameter& dGamma, MINT::FitParameter& dm, MINT::FitParameter& eff_tag, MINT::FitParameter& w ):
1113  _amps1(amps1),_amps2(amps2),_ampsSum(ampsSum),_r(r),_delta(delta),_gamma(gamma),_tau(tau),_dGamma(dGamma),_dm(dm),_eff_tag(eff_tag),_w(w),
1114  _intA(-1),_intAbar(-1),_intAAbar(-1)
1115  {;}
1116 };
1117 
1118 // DsK like time PDF with additional coherence factor
1119 class TimePdf : public MINT::PdfBase<IDalitzEvent>
1120 {
1121 protected:
1122  FitParameter& _C;
1123  FitParameter& _D;
1124  FitParameter& _D_bar;
1125  FitParameter& _S;
1126  FitParameter& _S_bar;
1127  FitParameter& _k;
1128 
1129  FitParameter& _tau;
1130  FitParameter& _dGamma;
1131  FitParameter& _dm;
1132  FitParameter& _eff_tag;
1133  FitParameter& _w;
1134 
1135 public:
1136  void parametersChanged(){
1137  }
1138  void beginFit(){
1139  }
1140  void endFit(){
1141  }
1142 
1143  inline double un_normalised(IDalitzEvent& evt){
1144  const double t = (double) evt.getValueFromVector(0);
1145  const double dt = (double) evt.getValueFromVector(1);
1146  const double q = static_cast<double>((int)evt.getValueFromVector(2));
1147  const double w = (double) evt.getValueFromVector(3);
1148  const double f = static_cast<double>((int)evt.getValueFromVector(4));
1149 
1150  const double e_eff = fabs(q)*_eff_tag + (1.-fabs(q))*(1.-_eff_tag);
1151 
1152  const double D = 1./2. * ((1.+f) * _D + (1.-f) * _D_bar);
1153  const double S = 1./2. * ((1.+f) * _S + (1.-f) * _S_bar);
1154 
1155  const double val = exp(-fabs(t)/(double)_tau) *
1156  (
1157  (2.-fabs(q))*cosh((double)_dGamma/2.*t)
1158  +f*q*(1.-2.*w)* _C *cos((double)_dm*t)
1159  -(2.-fabs(q))*2.0*_k* D *sinh((double)_dGamma/2.*t)
1160  -f*2.0*q*(1.-2.*w)*_k* S *sin((double)_dm*t)
1161  )*e_eff;
1162 
1163  return val;
1164  }
1165 
1166  virtual double getVal(IDalitzEvent& evt){
1167 
1168  const double f = static_cast<double>((int)evt.getValueFromVector(4));
1169  const double D = 1./2. * ((1.+f) * _D + (1.-f) * _D_bar);
1170 
1171  const double Gamma = 1./((double) _tau);
1172 
1173  double val = un_normalised(evt);
1174 
1175  double norm = 2.* ((4.*Gamma/(4.*Gamma*Gamma-_dGamma*_dGamma)) - 2.0*_k* D * (2.*_dGamma/(4.*Gamma*Gamma-_dGamma*_dGamma)));
1176 
1177  return val/norm;
1178  }
1179 
1180  virtual double getVal_withPs(IDalitzEvent& evt){return getVal(evt);}
1181  virtual double getVal_noPs(IDalitzEvent& evt){return getVal(evt);}
1182 
1183  virtual double getVal(IDalitzEvent* evt){
1184  if(0 == evt) return 0;
1185  return getVal(*evt);
1186  }
1187  virtual double getVal_withPs(IDalitzEvent* evt){
1188  if(0 == evt) return 0;
1189  return getVal_withPs(*evt);
1190  }
1191  virtual double getVal_noPs(IDalitzEvent* evt){
1192  if(0 == evt) return 0;
1193  return getVal_noPs(*evt);
1194  }
1195 
1196  TimePdf(MINT::FitParameter& C,MINT::FitParameter& D, MINT::FitParameter& D_bar,MINT::FitParameter& S, MINT::FitParameter& S_bar,MINT::FitParameter& k, MINT::FitParameter& tau, MINT::FitParameter& dGamma, MINT::FitParameter& dm, MINT::FitParameter& eff_tag,MINT::FitParameter& w ):
1197  _C(C),_D(D),_D_bar(D_bar),_S(S),_S_bar(S_bar),_k(k),_tau(tau),_dGamma(dGamma),_dm(dm),_eff_tag(eff_tag), _w(w) {;}
1198 
1199 };
1200 
1201 // Time PDF in term of r, gamma, delta instead of CP coefficients
1202 class TimePdf_mod : public MINT::PdfBase<IDalitzEvent>
1203 {
1204 protected:
1205  FitParameter& _r;
1206  FitParameter& _delta;
1207  FitParameter& _gamma;
1208  FitParameter& _k;
1209 
1210  FitParameter& _tau;
1211  FitParameter& _dGamma;
1212  FitParameter& _dm;
1213  FitParameter& _eff_tag;
1214  FitParameter& _w;
1215 
1216 public:
1217  void parametersChanged(){
1218  }
1219  void beginFit(){
1220  }
1221  void endFit(){
1222  }
1223 
1224  inline double un_normalised(IDalitzEvent& evt){
1225  const double t = (double) evt.getValueFromVector(0);
1226  const double dt = (double) evt.getValueFromVector(1);
1227  const double q = static_cast<double>((int)evt.getValueFromVector(2));
1228  const double w = (double) evt.getValueFromVector(3);
1229  const double f = static_cast<double>((int)evt.getValueFromVector(4));
1230 
1231  const double e_eff = fabs(q)*_eff_tag + (1.-fabs(q))*(1.-_eff_tag);
1232 
1233  const double C = (1.-(double)_r*(double)_r)/(1.+(double)_r*(double)_r);
1234  const double D = (double)_r * cos(((double)_delta-(double)_gamma*f)/360.*2*pi)/(1.+(double)_r*(double)_r);
1235  const double S = (double)_r * sin(((double)_delta-(double)_gamma*f)/360.*2*pi)/(1.+(double)_r*(double)_r);
1236 
1237  const double val = exp(-fabs(t)/(double)_tau) *
1238  (
1239  (2.-fabs(q))*cosh((double)_dGamma/2.*t)
1240  +f*q*(1.-2.*w)* C *cos((double)_dm*t)
1241  -(2.-fabs(q))*2.0*_k* D *sinh((double)_dGamma/2.*t)
1242  -f*2.0*q*(1.-2.*w)*_k* S *sin((double)_dm*t)
1243  )*e_eff;
1244 
1245  return val;
1246  }
1247 
1248  virtual double getVal(IDalitzEvent& evt){
1249 
1250  const double f = static_cast<double>((int)evt.getValueFromVector(4));
1251  const double D = (double)_r * cos(((double)_delta-(double)_gamma*f)/360.*2*pi)/(1.+(double)_r*(double)_r);
1252 
1253  const double Gamma = 1./((double) _tau);
1254 
1255  double val = un_normalised(evt);
1256 
1257  double norm = 2.* ((4.*Gamma/(4.*Gamma*Gamma-_dGamma*_dGamma)) - 2.0*_k* D * (2.*_dGamma/(4.*Gamma*Gamma-_dGamma*_dGamma)));
1258 
1259  return val/norm;
1260  }
1261 
1262  virtual double getVal_withPs(IDalitzEvent& evt){return getVal(evt);}
1263  virtual double getVal_noPs(IDalitzEvent& evt){return getVal(evt);}
1264 
1265  virtual double getVal(IDalitzEvent* evt){
1266  if(0 == evt) return 0;
1267  return getVal(*evt);
1268  }
1269  virtual double getVal_withPs(IDalitzEvent* evt){
1270  if(0 == evt) return 0;
1271  return getVal_withPs(*evt);
1272  }
1273  virtual double getVal_noPs(IDalitzEvent* evt){
1274  if(0 == evt) return 0;
1275  return getVal_noPs(*evt);
1276  }
1277 
1278  TimePdf_mod(MINT::FitParameter& r,MINT::FitParameter& delta, MINT::FitParameter& gamma,MINT::FitParameter& k, MINT::FitParameter& tau, MINT::FitParameter& dGamma, MINT::FitParameter& dm, MINT::FitParameter& eff_tag,MINT::FitParameter& w ):
1279  _r(r),_delta(delta),_gamma(gamma),_k(k),_tau(tau),_dGamma(dGamma),_dm(dm),_eff_tag(eff_tag), _w(w) {;}
1280 
1281 };
1282 
1283 
1284 int ampFit(int step=0){
1285 
1286  // Set random seed, important for toy studies
1287  TRandom3 ranLux;
1288  NamedParameter<int> RandomSeed("RandomSeed", 0);
1289  int seed = RandomSeed + step;
1290  ranLux.SetSeed((int)seed);
1291  gRandom = &ranLux;
1292 
1293  // Generate list of amplitudes
1294  FitAmplitude::AutogenerateFitFile();
1295 
1296  // Read options from steering file
1297  NamedParameter<string> InputFileName("InputFileName", (std::string) "");
1298  std::string inputFile = InputFileName;
1299  bool generateNew = (std::string) InputFileName == "";
1300  std::cout << "InputFileName: " << InputFileName << std::endl;
1301 
1302  NamedParameter<string> OutputDir("OutputDir", (std::string) "", (char*) 0);
1303  NamedParameter<string> OutputRootFile("OutputRootFile", (std::string) "OutputRootFile.root", (char*) 0);
1304 
1305  NamedParameter<int> EvtGenTest("EvtGenTest", 0);
1306  NamedParameter<string> IntegratorEventFile("IntegratorEventFile", (std::string) "SignalIntegrationEvents.root", (char*) 0);
1307  TString integratorEventFile = (string) IntegratorEventFile;
1308  if(EvtGenTest) integratorEventFile.ReplaceAll(".root",("_" + anythingToString(seed) + ".root").c_str() );
1309 
1310 
1311  NamedParameter<double> integPrecision("IntegPrecision", 1.e-2);
1312  NamedParameter<std::string> integMethod("IntegMethod", (std::string)"efficient");
1313 
1314  NamedParameter<int> EventPattern("Event Pattern", 521, 321, 211, -211, 443);
1315  DalitzEventPattern pat(EventPattern.getVector());
1316  cout << " got event pattern: " << pat << endl;
1317 
1318  NamedParameter<int> Nevents("Nevents", 100);
1319  NamedParameter<int> saveEvents("saveEvents", 1);
1320  NamedParameter<double> pdf_max("pdf_max", 100);
1321 
1322  NamedParameter<int> doPlots("doPlots", 1);
1323  NamedParameter<int> do2DScan("do2DScan", 0);
1324  NamedParameter<int> doTimeFit("doTimeFit", 1);
1325  NamedParameter<int> doDalitzFit("doDalitzFit", 1);
1326 
1327  vector<int> s123;
1328  s123.push_back(1);
1329  s123.push_back(2);
1330  s123.push_back(3);
1331 
1332  vector<int> s234;
1333  s234.push_back(2);
1334  s234.push_back(3);
1335  s234.push_back(4);
1336 
1337  vector<int> s134;
1338  s134.push_back(1);
1339  s134.push_back(3);
1340  s134.push_back(4);
1341 
1342  vector<int> s124;
1343  s124.push_back(1);
1344  s124.push_back(2);
1345  s124.push_back(4);
1346 
1347  // Fit parameters for time-dependent part
1348  FitParameter tau("tau");
1349  FitParameter dGamma("dGamma");
1350  FitParameter dm("dm");
1351 
1352  FitParameter r("ratio");
1353  FitParameter delta("delta");
1354  FitParameter gamma("gamma");
1355 
1356  FitParameter eff_tag("eff_tag");
1357  FitParameter mistag("mistag");
1358 
1359  // Define amplitude model
1360  DalitzEventList eventListPhsp,eventList;
1361  DalitzEventList eventList_f, eventList_fbar;
1362 
1363  eventListPhsp.generatePhaseSpaceEvents(100000,pat);
1364 
1365  FitAmpSum fas_tmp((DalitzEventPattern)pat);
1366  fas_tmp.getVal(eventListPhsp[0]);
1367  fas_tmp.normalizeAmps(eventListPhsp);
1368 
1369  counted_ptr<FitAmpList> List_1 = fas_tmp.GetCloneOfSubsetSameFitParameters("K(1)(1270)+");
1370  FitAmpSum fas(*List_1);
1371  FitAmpSum fasCC(*List_1);
1372  fasCC.CPConjugateSameFitParameters();
1373  fasCC.CConjugateFinalStateSameFitParameters();
1374  FitParameter r_K1_re("r_K1_Re",2,0,0.01);
1375  FitParameter r_K1_im("r_K1_Im",2,0,0.01);
1376  counted_ptr<IReturnComplex> r_K1_plus = new CPV_amp_polar(r_K1_re,r_K1_im,1);
1377  counted_ptr<IReturnComplex> r_K1_minus = new CPV_amp_polar(r_K1_re,r_K1_im,-1);
1378  fas.multiply(r_K1_plus);
1379  fasCC.multiply(r_K1_minus);
1380 
1381  AddScaledAmpsToList(fas_tmp, fas, fasCC, "K(1)(1400)+", r_K1_plus, r_K1_minus );
1382 
1383  FitParameter r_2_re("r_2_Re",2,0,0.01);
1384  FitParameter r_2_im("r_2_Im",2,0,0.01);
1385  counted_ptr<IReturnComplex> r_2_plus = new CPV_amp_polar(r_2_re,r_2_im,1);
1386  counted_ptr<IReturnComplex> r_2_minus = new CPV_amp_polar(r_2_re,r_2_im,-1);
1387  AddScaledAmpsToList(fas_tmp, fas, fasCC, "K*(1410)+", r_2_plus, r_2_minus );
1388 
1389  FitParameter r_3_re("r_3_Re",2,0,0.01);
1390  FitParameter r_3_im("r_3_Im",2,0,0.01);
1391  counted_ptr<IReturnComplex> r_3_plus = new CPV_amp_polar(r_3_re,r_3_im,1);
1392  counted_ptr<IReturnComplex> r_3_minus = new CPV_amp_polar(r_3_re,r_3_im,-1);
1393  AddScaledAmpsToList(fas_tmp, fas, fasCC, "NonResV0(->Ds-,pi+),K*(892)0(->K+,pi-)", r_3_plus, r_3_minus );
1394  AddScaledAmpsToList(fas_tmp, fas, fasCC, "NonResS0(->Ds-,pi+),K*(892)0(->K+,pi-)", r_3_plus, r_3_minus );
1395  AddScaledAmpsToList(fas_tmp, fas, fasCC, "NonResV0(->Ds-,K+),sigma10(->pi+,pi-)", r_3_plus, r_3_minus );
1396 
1397  /*
1398  FitParameter r_4_re("r_4_Re",2,0,0.01);
1399  FitParameter r_4_im("r_4_Im",2,0,0.01);
1400  counted_ptr<IReturnComplex> r_4_plus = new CPV_amp_polar(r_4_re,r_4_im,1);
1401  counted_ptr<IReturnComplex> r_4_minus = new CPV_amp_polar(r_4_re,r_4_im,-1);
1402  AddScaledAmpsToList(fas_tmp, fas, fasCC, "NonResA0(->sigma10(->pi+,pi-),Ds-)", r_4_plus, r_4_minus );
1403  AddScaledAmpsToList(fas_tmp, fas, fasCC, "NonResA0(->rho(770)0(->pi+,pi-),Ds-),K+", r_4_plus, r_4_minus );
1404  */
1405  vector<double> k_gen = coherenceFactor(fas,fasCC,(double)r, (double)delta,eventListPhsp);
1406 
1407  counted_ptr<FitAmpList> sumList = fas.GetCloneSameFitParameters();
1408  FitAmpSum fas_sum(*sumList);
1409  fas_sum.addAsList(fasCC,1.);
1410  fas_sum.getVal(eventListPhsp[0]);
1411 
1412  AmpsPdfFlexiFast ampsSig(pat, &fas, 0, integPrecision,integMethod, (std::string) integratorEventFile);
1413  AmpsPdfFlexiFast ampsSigCC(pat, &fasCC, 0, integPrecision,integMethod, (std::string) integratorEventFile);
1414  AmpsPdfFlexiFast ampsSum(pat, &fas_sum, 0, integPrecision,integMethod, (std::string) integratorEventFile);
1415 
1416  // Make full time-dependent PDF
1417  AmpsPdfFlexiFastCPV pdf(&ampsSig,&ampsSigCC,&ampsSum, r, delta, gamma, tau, dGamma, dm, eff_tag, mistag );
1418 
1419  // Generate toys
1420  if(generateNew){
1421  time_t startTime = time(0);
1422 
1423  TFile* fileTD= TFile::Open("dummy.root","RECREATE");
1424  fileTD->cd();
1425 
1426  TTree *tree = new TTree("TD","TD");
1427  double t,dt,w;
1428  int q,f;
1429  tree->Branch("t",&t,"t/D");
1430  tree->Branch("dt",&dt,"dt/D");
1431  tree->Branch("q",&q,"q/I");
1432  tree->Branch("w",&w,"w/D");
1433  tree->Branch("f",&f,"f/I");
1434 
1435  // This returns the incoherent sum of the amplitudes, ie, the sum of their magnitudes, rather than
1436  // the magnitude of the sum of their complex values, like FitAmpSum.
1437  // So, gives the maximum possible value of the PDF at any decay time?
1438  // But it's only for one flavour?
1439  FitAmpIncoherentSum fasGen((DalitzEventPattern)pat);
1440  fasGen.getVal(eventListPhsp[0]);
1441  fasGen.normalizeAmps(eventListPhsp);
1442  SignalGenerator sg(pat,&fasGen);
1443  fasGen.print();
1444 
1445  //simple hit and miss
1446  for(int i = 0; i < Nevents; i++){
1447  while(true){
1448  // Simple exponential. Where does the mixing come in? At the accept/reject stage?
1449  t = ranLux.Exp(tau);
1450  dt = ranLux.Uniform(0.,0.25);
1451  //pdf.getSmearedTime(t,dt,ranLux);
1452  // Initial flavour.
1453  const double q_rand = ranLux.Uniform();
1454  q = 0;
1455  if (q_rand < 1./3.) q = -1;
1456  if (q_rand > 2./3.) q = 1;
1457  // Mistag
1458  w = mistag;
1459  // Final state (Ds-K+pipi or Ds+K-pipi)
1460  const double f_rand = ranLux.Uniform();
1461  if(f_rand < 0.5)f = 1;
1462  else f = -1;
1463  //if(i < Nevents/2)f = 1;
1464  //else f = -1;
1465 
1466  counted_ptr<IDalitzEvent> evtPtr(sg.newEvent());
1467  DalitzEvent evt(evtPtr.get());
1468  //if(!(sqrt(evt.sij(s234)/(GeV*GeV)) < 1.95 && sqrt(evt.s(2,4)/(GeV*GeV)) < 1.2 && sqrt(evt.s(3,4)/(GeV*GeV)) < 1.2))continue;
1469  // pdf_max is user configurable, seems weird ... I guess it saves computationally finding the maximum?
1470  // What is evt.getGeneratorPdfRelativeToPhaseSpace() returning?
1471  double maxVal = evt.getGeneratorPdfRelativeToPhaseSpace()*exp(-fabs(t)/(tau))/(tau)*pdf_max;
1472 
1473  // These are accessed by AmpsPdfFlexiFastCPV::un_normalised_noPs to calculate the PDF value.
1474  evt.setValueInVector(0, t);
1475  evt.setValueInVector(1, dt);
1476  evt.setValueInVector(2, q);
1477  evt.setValueInVector(3, w);
1478  evt.setValueInVector(4, f);
1479 
1480  // This just returns the decay rate using Eq 2.1 of the note.
1481  const double pdfVal = pdf.un_normalised_noPs(evt);
1482  //const double pdfVal = pdf.getValForGeneration(evt);
1483 
1484  // Random number to determine if it's accepted.
1485  const double height = ranLux.Uniform(0,maxVal);
1486  //Safety check on the maxmimum generated height
1487  if( pdfVal > maxVal ){
1488  std::cout << "ERROR: PDF above determined maximum." << std::endl;
1489  std::cout << pdfVal << " > " << maxVal << std::endl;
1490  //exit(1);
1491  pdf_max = pdf_max * 2.;
1492  }
1493 
1494  //Hit-and-miss
1495  if( height < pdfVal ){
1496  eventList.Add(evt);
1497  if(f==1)eventList_f.Add(evt);
1498  else eventList_fbar.Add(evt);
1499  tree->Fill();
1500  break;
1501  }
1502  }
1503  }
1504 
1505  cout << " Generated " << Nevents << " Events. Took " << (time(0) - startTime)/60 << " mins "<< endl;
1506 
1507  if(saveEvents){
1508  // Very messy workaround, needs to be fixded
1509  // event vector should be saved with evt
1510  eventList.saveAsNtuple(OutputRootFile);
1511  TFile* file= TFile::Open(((std::string) OutputRootFile).c_str(),"UPDATE");
1512  file->cd();
1513  tree->CloneTree()->Write();
1514  file->Close();
1515  fileTD->Close();
1516  }
1517  }
1518  else{
1519 
1520  cout << "reading events from file " << inputFile << endl;
1521 
1522  TFile *_InputFile = TFile::Open(inputFile.c_str());
1523  TTree* in_tree, *in_treeTD;
1524  in_tree=dynamic_cast<TTree*>(_InputFile->Get("DalitzEventList"));
1525  eventList.fromNtuple(in_tree,1);
1526  in_treeTD=dynamic_cast<TTree*>(_InputFile->Get("TD"));
1527  double t,dt,w;
1528  int q,f;
1529  in_treeTD->SetBranchAddress("t",&t);
1530  in_treeTD->SetBranchAddress("dt",&dt);
1531  in_treeTD->SetBranchAddress("q",&q);
1532  in_treeTD->SetBranchAddress("w",&w);
1533  in_treeTD->SetBranchAddress("f",&f);
1534 
1535  // Very messy workaround, needs to be fixded
1536  // event vector should be saved with evt
1537  if(in_treeTD->GetEntries() != in_tree->GetEntries()){cout << "ERROR: Different number of DalitzEvents and time information " << endl; throw "crash"; }
1538 
1539  for (unsigned int i = 0; i < in_treeTD->GetEntries(); i++) {
1540  in_treeTD->GetEntry(i);
1541 
1542  eventList[i].setValueInVector(0, t);
1543  eventList[i].setValueInVector(1, dt);
1544  eventList[i].setValueInVector(2, q);
1545  eventList[i].setValueInVector(3, w);
1546  eventList[i].setValueInVector(4, f);
1547 
1548  if(f==1)eventList_f.Add(eventList[i]);
1549  else eventList_fbar.Add(eventList[i]);
1550  }
1551 
1552  cout << " I've got " << eventList.size() << " events." << endl;
1553  _InputFile->Close();
1554  }
1555 
1556  // Fit
1557  Neg2LL fcn(pdf, eventList_f);
1558  Neg2LL fcn_bar(pdf, eventList_fbar);
1559 
1560  Neg2LL neg2LL(pdf, eventList);
1561  //Neg2LLSum neg2LLSum(&neg2LL);
1562  Neg2LLSum neg2LLSum(&fcn,&fcn_bar);
1563  //neg2LLSum.addConstraints();
1564  Minimiser mini(&neg2LLSum);
1565  mini.doFit();
1566  mini.printResultVsInput();
1567 
1568  // Calculate pulls
1569  gDirectory->cd();
1570  TFile* paraFile = new TFile(((string)OutputDir+"pull_"+anythingToString((int)seed)+".root").c_str(), "RECREATE");
1571  paraFile->cd();
1572  TNtupleD* ntp=0;
1573  TTree* tree = new TTree("Coherence","Coherence");
1574  double r_val,delta_val,gamma_val,k_val,n2ll;
1575  TBranch* br_r = tree->Branch( "r", &r_val, "r_val/D" );
1576  TBranch* br_delta = tree->Branch( "delta", &delta_val, "delta_val/D" );
1577  TBranch* br_gamma = tree->Branch( "gamma", &gamma_val, "gamma_val/D" );
1578  TBranch* br_k = tree->Branch( "k", &k_val, "k_val/D" );
1579  TBranch* br_n2ll = tree->Branch( "n2ll", &n2ll, "n2ll/D" );
1580 
1581  double r_val_t,delta_val_t,gamma_val_t,k_val_t,n2ll_t;
1582  TBranch* br_r_t = tree->Branch( "r_t", &r_val_t, "r_val_t/D" );
1583  TBranch* br_delta_t = tree->Branch( "delta_t", &delta_val_t, "delta_val_t/D" );
1584  TBranch* br_gamma_t = tree->Branch( "gamma_t", &gamma_val_t, "gamma_val_t/D" );
1585  TBranch* br_k_t = tree->Branch( "k_t", &k_val_t, "k_val_t/D" );
1586  TBranch* br_n2ll_t = tree->Branch( "n2ll_t", &n2ll_t, "n2ll_t/D" );
1587 
1588  double r_val_t_fixed_k,delta_val_t_fixed_k,gamma_val_t_fixed_k,k_val_t_fixed_k,n2ll_t_fixed_k;
1589  TBranch* br_r_t_fixed_k = tree->Branch( "r_t_fixed_k", &r_val_t_fixed_k, "r_val_t_fixed_k/D" );
1590  TBranch* br_delta_t_fixed_k = tree->Branch( "delta_t_fixed_k", &delta_val_t_fixed_k, "delta_val_t_fixed_k/D" );
1591  TBranch* br_gamma_t_fixed_k = tree->Branch( "gamma_t_fixed_k", &gamma_val_t_fixed_k, "gamma_val_t_fixed_k/D" );
1592  TBranch* br_k_t_fixed_k = tree->Branch( "k_t_fixed_k", &k_val_t_fixed_k, "k_val_t_fixed_k/D" );
1593  TBranch* br_n2ll_t_fixed_k = tree->Branch( "n2ll_t_fixed_k", &n2ll_t_fixed_k, "n2ll_t_fixed_k/D" );
1594 
1595  double r_val_time_integrated,delta_val_time_integrated,gamma_val_time_integrated,k_val_time_integrated,n2ll_time_integrated;
1596  TBranch* br_r_time_integrated = tree->Branch( "r_time_integrated", &r_val_time_integrated, "r_val_time_integrated/D" );
1597  TBranch* br_delta_time_integrated = tree->Branch( "delta_time_integrated", &delta_val_time_integrated, "delta_val_time_integrated/D" );
1598  TBranch* br_gamma_time_integrated = tree->Branch( "gamma_time_integrated", &gamma_val_time_integrated, "gamma_val_time_integrated/D" );
1599  TBranch* br_k_time_integrated = tree->Branch( "k_time_integrated", &k_val_time_integrated, "k_val_time_integrated/D" );
1600  TBranch* br_n2ll_time_integrated = tree->Branch( "n2ll_time_integrated", &n2ll_time_integrated, "n2ll_time_integrated/D" );
1601 
1602  TBranch* br_seed = tree->Branch( "seed", &seed, "seed/I" );
1603 
1604  MinuitParameterSet::getDefaultSet()->fillNtp(paraFile, ntp);
1605  ntp->AutoSave();
1606 
1607  n2ll = neg2LLSum.getVal();
1608 
1609  for(unsigned int i=0; i < MinuitParameterSet::getDefaultSet()->size(); i++){
1610  if(0 == MinuitParameterSet::getDefaultSet()->getParPtr(i)) continue;
1611  if(A_is_in_B("gamma",MinuitParameterSet::getDefaultSet()->getParPtr(i)->name()))gamma_val = MinuitParameterSet::getDefaultSet()->getParPtr(i)->mean();
1612  }
1613 
1614  vector<double> k_fit = coherenceFactor(fas,fasCC,(double)r, (double)delta,eventListPhsp);
1615  r_val = k_fit[0];
1616  k_val = k_fit[1];
1617  delta_val = k_fit[2];
1618 
1619  // Calculate fit fractions
1620  pdf.doFinalStatsAndSaveForAmp12(&mini,((string)OutputDir+"FitAmpResults_rand_"+anythingToString((int)seed)).c_str(),((string)OutputDir+"fitFractions_"+anythingToString((int)seed)).c_str());
1621 
1622  if(doDalitzFit==1){
1623 
1624  for(unsigned int i=0; i < MinuitParameterSet::getDefaultSet()->size(); i++){
1625  if(0 == MinuitParameterSet::getDefaultSet()->getParPtr(i)) continue;
1626  MinuitParameterSet::getDefaultSet()->getParPtr(i)->setCurrentFitVal(MinuitParameterSet::getDefaultSet()->getParPtr(i)->meanInit());
1627  }
1628 
1629  AmpsPdfFlexiFastCPV_time_integrated pdf_time_integrated(&ampsSig,&ampsSigCC,&ampsSum, r, delta, gamma, tau, dGamma, dm, eff_tag, mistag );
1630  Neg2LL fcn_time_integrated(pdf_time_integrated, eventList_f);
1631  Neg2LL fcn_bar_time_integrated(pdf_time_integrated, eventList_fbar);
1632  Neg2LLSum neg2LLSum_time_integrated(&fcn_time_integrated,&fcn_bar_time_integrated);
1633  Minimiser mini_time_integrated(&neg2LLSum_time_integrated);
1634  mini_time_integrated.doFit();
1635  mini_time_integrated.printResultVsInput();
1636 
1637  n2ll_time_integrated = neg2LLSum_time_integrated.getVal();
1638 
1639  for(unsigned int i=0; i < MinuitParameterSet::getDefaultSet()->size(); i++){
1640  if(0 == MinuitParameterSet::getDefaultSet()->getParPtr(i)) continue;
1641  if(A_is_in_B("gamma",MinuitParameterSet::getDefaultSet()->getParPtr(i)->name()))gamma_val_time_integrated = MinuitParameterSet::getDefaultSet()->getParPtr(i)->mean();
1642  }
1643 
1644  vector<double> k_fit_time_integrated = coherenceFactor(fas,fasCC,(double)r, (double)delta,eventListPhsp);
1645  r_val_time_integrated = k_fit_time_integrated[0];
1646  k_val_time_integrated = k_fit_time_integrated[1];
1647  delta_val_time_integrated = k_fit_time_integrated[2];
1648  }
1649 
1650  if(doTimeFit==1){
1651  for(unsigned int i=0; i < MinuitParameterSet::getDefaultSet()->size(); i++){
1652  if(0 == MinuitParameterSet::getDefaultSet()->getParPtr(i)) continue;
1653  MinuitParameterSet::getDefaultSet()->getParPtr(i)->setCurrentFitVal(MinuitParameterSet::getDefaultSet()->getParPtr(i)->meanInit());
1654  }
1655 
1656  for(unsigned int i=0; i < MinuitParameterSet::getDefaultSet()->size(); i++){
1657  if(0 == MinuitParameterSet::getDefaultSet()->getParPtr(i)) continue;
1658  if(A_is_in_B("_Re",MinuitParameterSet::getDefaultSet()->getParPtr(i)->name())){
1659  MinuitParameterSet::getDefaultSet()->unregister(MinuitParameterSet::getDefaultSet()->getParPtr(i));
1660  i= 0;
1661  }
1662  else if(A_is_in_B("_Im",MinuitParameterSet::getDefaultSet()->getParPtr(i)->name())){
1663  MinuitParameterSet::getDefaultSet()->unregister(MinuitParameterSet::getDefaultSet()->getParPtr(i));
1664  i = 0;
1665  }
1666  else if(A_is_in_B("ratio",MinuitParameterSet::getDefaultSet()->getParPtr(i)->name())){
1667  MinuitParameterSet::getDefaultSet()->unregister(MinuitParameterSet::getDefaultSet()->getParPtr(i));
1668  i = 0;
1669  }
1670  }
1671 
1672  FitParameter r2("ratio2",0,k_gen[0],0.01);
1673  FitParameter k("kappa",0,k_gen[1],0.05,0.,1.);
1674  TimePdf_mod t_pdf(r2,delta,gamma,k,tau,dGamma,dm,eff_tag,mistag );
1675  Neg2LL fcn_t(t_pdf, eventList_f);
1676  Neg2LL fcn_t_bar(t_pdf, eventList_fbar);
1677  Neg2LLSum neg2LLSum_t(&fcn_t,&fcn_t_bar);
1678 
1679  Minimiser mini_t(&neg2LLSum_t);
1680  mini_t.doFit();
1681  mini_t.printResultVsInput();
1682  n2ll_t = neg2LLSum_t.getVal();
1683 
1684  for(unsigned int i=0; i < MinuitParameterSet::getDefaultSet()->size(); i++){
1685  if(0 == MinuitParameterSet::getDefaultSet()->getParPtr(i)) continue;
1686  if(A_is_in_B("gamma",MinuitParameterSet::getDefaultSet()->getParPtr(i)->name()))gamma_val_t= MinuitParameterSet::getDefaultSet()->getParPtr(i)->mean();
1687  if(A_is_in_B("ratio2",MinuitParameterSet::getDefaultSet()->getParPtr(i)->name()))r_val_t= MinuitParameterSet::getDefaultSet()->getParPtr(i)->mean();
1688  if(A_is_in_B("delta",MinuitParameterSet::getDefaultSet()->getParPtr(i)->name()))delta_val_t= MinuitParameterSet::getDefaultSet()->getParPtr(i)->mean();
1689  if(A_is_in_B("kappa",MinuitParameterSet::getDefaultSet()->getParPtr(i)->name()))k_val_t= MinuitParameterSet::getDefaultSet()->getParPtr(i)->mean();
1690  }
1691 
1692  for(unsigned int i=0; i < MinuitParameterSet::getDefaultSet()->size(); i++){
1693  if(0 == MinuitParameterSet::getDefaultSet()->getParPtr(i)) continue;
1694  MinuitParameterSet::getDefaultSet()->getParPtr(i)->setCurrentFitVal(MinuitParameterSet::getDefaultSet()->getParPtr(i)->meanInit());
1695  }
1696  k.setCurrentFitVal(k_gen[1]);
1697  k.fix();
1698  mini_t.doFit();
1699  mini_t.printResultVsInput();
1700  n2ll_t_fixed_k = neg2LLSum_t.getVal();
1701 
1702  for(unsigned int i=0; i < MinuitParameterSet::getDefaultSet()->size(); i++){
1703  if(0 == MinuitParameterSet::getDefaultSet()->getParPtr(i)) continue;
1704  if(A_is_in_B("gamma",MinuitParameterSet::getDefaultSet()->getParPtr(i)->name()))gamma_val_t_fixed_k= MinuitParameterSet::getDefaultSet()->getParPtr(i)->mean();
1705  if(A_is_in_B("ratio2",MinuitParameterSet::getDefaultSet()->getParPtr(i)->name()))r_val_t_fixed_k= MinuitParameterSet::getDefaultSet()->getParPtr(i)->mean();
1706  if(A_is_in_B("delta",MinuitParameterSet::getDefaultSet()->getParPtr(i)->name()))delta_val_t_fixed_k= MinuitParameterSet::getDefaultSet()->getParPtr(i)->mean();
1707  if(A_is_in_B("kappa",MinuitParameterSet::getDefaultSet()->getParPtr(i)->name()))k_val_t_fixed_k= MinuitParameterSet::getDefaultSet()->getParPtr(i)->mean();
1708  }
1709  }
1710 
1711  tree->Fill();
1712  paraFile->cd();
1713  tree->SetDirectory(paraFile);
1714  tree->Write();
1715  paraFile->Close();
1716  delete paraFile;
1717 
1718  if(doPlots==1){
1719  cout << "Now plotting:" << endl;
1720  //DalitzHistoSet datH = eventList.histoSet();
1721  //DalitzHistoSet fitH = ampsSig.histoSet();
1722  //datH.drawWithFitNorm(fitH, ((string)OutputDir+(string)"datFit_"+anythingToString((int)RandomSeed)+"_").c_str(),"eps");
1723 
1724  int nBinst = 60;
1725  int nBins = 50;
1726 
1727  TH1D* h_t = new TH1D("",";t (ps);Events (norm.) ",nBinst,0,6);
1728  TH1D* s_Kpipi = new TH1D("",";#left[m^{2}(K^{+} #pi^{+} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0.8,12);
1729  TH1D* s_Kpi = new TH1D("",";#left[m^{2}(K^{+} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0.,4);
1730  TH1D* s_pipi = new TH1D("",";#left[m^{2}(#pi^{+} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0,2);
1731  TH1D* s_Dspipi = new TH1D("",";#left[m^{2}(D_{s}^{-} #pi^{+} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0,30);
1732  TH1D* s_DsK = new TH1D("",";#left[m^{2}(D_{s}^{-} K^{+})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0,30) ;
1733  TH1D* s_DsKpi = new TH1D("",";#left[m^{2}(D_{s}^{-} K^{+} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,5,30);
1734  TH1D* s_Dspi = new TH1D("",";#left[m^{2}(D_{s}^{-} #pi^{+})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0,25);
1735  TH1D* s_Dspim = new TH1D("",";#left[m^{2}(D_{s}^{-} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0,25);
1736 
1737  for (int i=0; i<eventList.size(); i++) {
1738  h_t->Fill(eventList[i].getValueFromVector(0));
1739  s_Kpipi->Fill(eventList[i].sij(s234)/(GeV*GeV));
1740  s_Kpi->Fill(eventList[i].s(2,4)/(GeV*GeV));
1741  s_pipi->Fill(eventList[i].s(3,4)/(GeV*GeV));
1742  s_Dspipi->Fill(eventList[i].sij(s134)/(GeV*GeV));
1743  s_DsK->Fill(eventList[i].s(1,2)/(GeV*GeV));
1744  s_DsKpi->Fill(eventList[i].sij(s124)/(GeV*GeV));
1745  s_Dspi->Fill(eventList[i].s(1,3)/(GeV*GeV));
1746  s_Dspim->Fill(eventList[i].s(1,4)/(GeV*GeV));
1747  }
1748 
1749  TH1D* h_t_fit = new TH1D("",";t",nBinst,0,6);
1750  TH1D* h_t_fit_mp = new TH1D("",";t",nBinst,0,6);
1751  TH1D* h_t_fit_0p = new TH1D("",";t",nBinst,0,6);
1752  TH1D* h_t_fit_pp = new TH1D("",";t",nBinst,0,6);
1753  TH1D* h_t_fit_mm = new TH1D("",";t",nBinst,0,6);
1754  TH1D* h_t_fit_0m = new TH1D("",";t",nBinst,0,6);
1755  TH1D* h_t_fit_pm = new TH1D("",";t",nBinst,0,6);
1756 
1757  TH1D* s_Kpipi_fit = new TH1D("",";#left[m^{2}(K^{+} #pi^{+} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0.8,12);
1758  TH1D* s_Kpi_fit = new TH1D("",";#left[m^{2}(K^{+} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0.,4);
1759  TH1D* s_pipi_fit = new TH1D("",";#left[m^{2}(K^{+} #pi^{+} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0,2);
1760  TH1D* s_Dspipi_fit = new TH1D("",";#left[m^{2}(D_{s}^{-} #pi^{+} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0,30);
1761  TH1D* s_DsK_fit = new TH1D("",";#left[m^{2}(D_{s}^{-} K^{+})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0,30);
1762  TH1D* s_DsKpi_fit = new TH1D("",";#left[m^{2}(D_{s}^{-} K^{+} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,5,30);
1763  TH1D* s_Dspi_fit = new TH1D("",";#left[m^{2}(D_{s}^{-} #pi^{+})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0,25);
1764  TH1D* s_Dspim_fit = new TH1D("",";#left[m^{2}(D_{s}^{-} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0,25);
1765 
1766  TH1D* s_Kpipi_fitBs = new TH1D("",";#left[m^{2}(K^{+} #pi^{+} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0.8,3.5);
1767  TH1D* s_Kpi_fitBs = new TH1D("",";#left[m^{2}(K^{+} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0.3,2.);
1768  TH1D* s_pipi_fitBs = new TH1D("",";#left[m^{2}(K^{+} #pi^{+} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0,2);
1769  TH1D* s_Dspipi_fitBs = new TH1D("",";#left[m^{2}(D_{s}^{-} #pi^{+} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",20,4,25);
1770  TH1D* s_DsK_fitBs = new TH1D("",";#left[m^{2}(D_{s}^{-} K^{+})#right] (GeV^{2}/c^{4});Events (norm.) ",20,4,26);
1771 
1772  TH1D* s_Kpipi_fitBsbar = new TH1D("",";#left[m^{2}(K^{+} #pi^{+} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0.8,3.5);
1773  TH1D* s_Kpi_fitBsbar = new TH1D("",";#left[m^{2}(K^{+} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0.3,2.);
1774  TH1D* s_pipi_fitBsbar = new TH1D("",";#left[m^{2}(K^{+} #pi^{+} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0,2);
1775  TH1D* s_Dspipi_fitBsbar = new TH1D("",";#left[m^{2}(D_{s}^{-} #pi^{+} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",20,4,25);
1776  TH1D* s_DsK_fitBsbar = new TH1D("",";#left[m^{2}(D_{s}^{-} K^{+})#right] (GeV^{2}/c^{4});Events (norm.) ",20,4,26);
1777 
1778  TH1D* s_Kpipi_fit_notag = new TH1D("",";#left[m^{2}(K^{+} #pi^{+} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0.8,3.5);
1779  TH1D* s_Kpi_fit_notag = new TH1D("",";#left[m^{2}(K^{+} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0.3,2.);
1780  TH1D* s_pipi_fit_notag = new TH1D("",";#left[m^{2}(K^{+} #pi^{+} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",nBins,0,2);
1781  TH1D* s_Dspipi_fit_notag = new TH1D("",";#left[m^{2}(D_{s}^{-} #pi^{+} #pi^{-})#right] (GeV^{2}/c^{4});Events (norm.) ",20,4,25);
1782  TH1D* s_DsK_fit_notag = new TH1D("",";#left[m^{2}(D_{s}^{-} K^{+})#right] (GeV^{2}/c^{4});Events (norm.) ",20,4,26);
1783 
1784  SignalGenerator sg(pat,&fas);
1785  sg.setWeighted();
1786 
1787  for(int i = 0; i < 2000000; i++){
1788 
1789  const double t = ranLux.Exp(tau);
1790  const double dt = ranLux.Uniform(0.00001,0.25);
1791  const double w = mistag;
1792 
1793  counted_ptr<IDalitzEvent> evtPtr(sg.newEvent());
1794  DalitzEvent evt(evtPtr.get());
1795 
1796 // if(!(sqrt(evt.sij(s234)/(GeV*GeV)) < 1.95 && sqrt(evt.s(2,4)/(GeV*GeV)) < 1.2 && sqrt(evt.s(3,4)/(GeV*GeV)) < 1.2))continue;
1797 
1798  evt.setValueInVector(0, t);
1799  evt.setValueInVector(1, dt);
1800  evt.setValueInVector(3, w);
1801 
1802  evt.setValueInVector(2, -1);
1803  evt.setValueInVector(4, 1);
1804  double weight_mp = pdf.un_normalised_noPs(evt)*evt.getWeight()/evt.getGeneratorPdfRelativeToPhaseSpace()/exp(-fabs(t)/(tau));
1805 
1806  evt.setValueInVector(2, 0);
1807  evt.setValueInVector(4, 1);
1808  double weight_0p = pdf.un_normalised_noPs(evt)*evt.getWeight()/evt.getGeneratorPdfRelativeToPhaseSpace()/exp(-fabs(t)/(tau));
1809 
1810  evt.setValueInVector(2, 1);
1811  evt.setValueInVector(4, 1);
1812  double weight_pp = pdf.un_normalised_noPs(evt)*evt.getWeight()/evt.getGeneratorPdfRelativeToPhaseSpace()/exp(-fabs(t)/(tau));
1813 
1814  evt.setValueInVector(2, -1);
1815  evt.setValueInVector(4, -1);
1816  double weight_mm = pdf.un_normalised_noPs(evt)*evt.getWeight()/evt.getGeneratorPdfRelativeToPhaseSpace()/exp(-fabs(t)/(tau));
1817 
1818  evt.setValueInVector(2, 0);
1819  evt.setValueInVector(4, -1);
1820  double weight_0m = pdf.un_normalised_noPs(evt)*evt.getWeight()/evt.getGeneratorPdfRelativeToPhaseSpace()/exp(-fabs(t)/(tau));
1821 
1822  evt.setValueInVector(2, 1);
1823  evt.setValueInVector(4, -1);
1824  double weight_pm = pdf.un_normalised_noPs(evt)*evt.getWeight()/evt.getGeneratorPdfRelativeToPhaseSpace()/exp(-fabs(t)/(tau));
1825 
1826  double weight = weight_mp + weight_0p + weight_pp + weight_mm + weight_0m + weight_pm ;
1827 
1828  h_t_fit->Fill(t,weight);
1829  h_t_fit_mp->Fill(t,weight_mp);
1830  h_t_fit_0p->Fill(t,weight_0p);
1831  h_t_fit_pp->Fill(t,weight_pp);
1832  h_t_fit_mm->Fill(t,weight_mm);
1833  h_t_fit_0m->Fill(t,weight_0m);
1834  h_t_fit_pm->Fill(t,weight_pm);
1835 
1836  double weight_B = fas.RealVal(evt)*evt.getWeight()/evt.getGeneratorPdfRelativeToPhaseSpace();
1837  double weight_Bbar = fasCC.RealVal(evt)*evt.getWeight()/evt.getGeneratorPdfRelativeToPhaseSpace();
1838  double weight_notag = weight_Bbar+ (double)r * (double)r * weight_Bbar;
1839 
1840  const complex<double> phase_diff = polar((double)r, (double)delta/360.*2*pi);
1841  double weight_coherence = std::abs(conj(fas.getVal(evt)) * fasCC.getVal(evt)*phase_diff) *evt.getWeight()/evt.getGeneratorPdfRelativeToPhaseSpace();
1842 
1843  s_Kpipi_fit->Fill(evt.sij(s234)/(GeV*GeV),weight);
1844  s_Kpi_fit->Fill(evt.s(2,4)/(GeV*GeV),weight);
1845  s_pipi_fit->Fill(evt.s(3,4)/(GeV*GeV),weight);
1846  s_Dspipi_fit->Fill(evt.sij(s134)/(GeV*GeV),weight);
1847  s_DsK_fit->Fill(evt.s(1,2)/(GeV*GeV),weight);
1848  s_DsKpi_fit->Fill(evt.sij(s124)/(GeV*GeV),weight);
1849  s_Dspi_fit->Fill(evt.s(1,3)/(GeV*GeV),weight);
1850  s_Dspim_fit->Fill(evt.s(1,4)/(GeV*GeV),weight);
1851 
1852  s_Kpipi_fitBs->Fill(evt.sij(s234)/(GeV*GeV),weight_B);
1853  s_Kpi_fitBs->Fill(evt.s(2,4)/(GeV*GeV),weight_B);
1854  s_pipi_fitBs->Fill(evt.s(3,4)/(GeV*GeV),weight_B);
1855  s_Dspipi_fitBs->Fill(evt.sij(s134)/(GeV*GeV),weight_B);
1856  s_DsK_fitBs->Fill(evt.s(1,2)/(GeV*GeV),weight_B);
1857 
1858  s_Kpipi_fitBsbar->Fill(evt.sij(s234)/(GeV*GeV),weight_Bbar);
1859  s_Kpi_fitBsbar->Fill(evt.s(2,4)/(GeV*GeV),weight_Bbar);
1860  s_pipi_fitBsbar->Fill(evt.s(3,4)/(GeV*GeV),weight_Bbar);
1861  s_Dspipi_fitBsbar->Fill(evt.sij(s134)/(GeV*GeV),weight_Bbar);
1862  s_DsK_fitBsbar->Fill(evt.s(1,2)/(GeV*GeV),weight_Bbar);
1863 
1864  s_Kpipi_fit_notag->Fill(evt.sij(s234)/(GeV*GeV),weight_notag);
1865  s_Kpi_fit_notag->Fill(evt.s(2,4)/(GeV*GeV),weight_notag);
1866  s_pipi_fit_notag->Fill(evt.s(3,4)/(GeV*GeV),weight_notag);
1867  s_Dspipi_fit_notag->Fill(evt.sij(s134)/(GeV*GeV),weight_notag);
1868  s_DsK_fit_notag->Fill(evt.s(1,2)/(GeV*GeV),weight_notag);
1869 
1870  }
1871 
1872  TCanvas* c = new TCanvas();
1873 
1874  h_t->SetLineColor(kBlack);
1875  h_t->DrawNormalized("e1",1);
1876  h_t_fit->SetLineColor(kBlue);
1877  h_t_fit->SetLineWidth(3);
1878  h_t_fit->DrawNormalized("histcsame",1);
1879 
1880  c->Print(((string)OutputDir+"h_t.eps").c_str());
1881  gPad->SetLogy(1);
1882  c->Print(((string)OutputDir+"h_t_log.eps").c_str());
1883  gPad->SetLogy(0);
1884 
1885  h_t->SetLineColor(kBlack);
1886  h_t->DrawNormalized("e1",1);
1887  h_t_fit->SetLineColor(kBlack);
1888  h_t_fit->SetLineWidth(3);
1889  h_t_fit->DrawNormalized("histcsame",1);
1890 
1891  h_t_fit_pm->SetLineColor(kBlue);
1892  h_t_fit_pm->SetLineWidth(3);
1893  h_t_fit_pm->DrawNormalized("histcsame",0.5*0.5*eff_tag);
1894  h_t_fit_mm->SetLineColor(kBlue);
1895  h_t_fit_mm->SetLineWidth(3);
1896  h_t_fit_mm->SetLineStyle(kDashed);
1897  h_t_fit_mm->DrawNormalized("histcsame",0.5*0.5*eff_tag);
1898 
1899  h_t_fit_pp->SetLineColor(kRed);
1900  h_t_fit_pp->SetLineWidth(3);
1901  h_t_fit_pp->DrawNormalized("histcsame",0.5*0.5*eff_tag);
1902  h_t_fit_mp->SetLineColor(kRed);
1903  h_t_fit_mp->SetLineWidth(3);
1904  h_t_fit_mp->SetLineStyle(kDashed);
1905  h_t_fit_mp->DrawNormalized("histcsame",0.5*0.5*eff_tag);
1906 
1907  h_t_fit_0p->SetLineColor(kGreen);
1908  h_t_fit_0p->SetLineWidth(3);
1909  h_t_fit_0p->DrawNormalized("histcsame",0.5*0.5*(1.-eff_tag));
1910  h_t_fit_0m->SetLineColor(kGreen);
1911  h_t_fit_0m->SetLineWidth(3);
1912  h_t_fit_0m->SetLineStyle(kDashed);
1913  h_t_fit_0m->DrawNormalized("histcsame",0.5*0.5*(1.-eff_tag));
1914 
1915  c->Print(((string)OutputDir+"h_t_2.eps").c_str());
1916  gPad->SetLogy(1);
1917  c->Print(((string)OutputDir+"h_t_2_log.eps").c_str());
1918  gPad->SetLogy(0);
1919 
1920  TLegend leg_t(0,0,1,1,"");
1921  leg_t.SetLineStyle(0);
1922  leg_t.SetLineColor(0);
1923  leg_t.SetFillColor(0);
1924  //leg_t.SetTextFont(22);
1925  leg_t.SetTextColor(1);
1926  //leg_t.SetTextSize(0.05);
1927  leg_t.SetTextAlign(12);
1928 
1929  leg_t.AddEntry(h_t_fit_pm,"B_{s} #rightarrow D_{s}^{-} K^{+} #pi^{+} #pi^{-}","l");
1930  leg_t.AddEntry(h_t_fit_mm,"#bar{B}_{s} #rightarrow D_{s}^{-} K^{+} #pi^{+} #pi^{-}","l");
1931  leg_t.AddEntry(h_t_fit_0m,"Untagged #rightarrow D_{s}^{-} K^{+} #pi^{+} #pi^{-}","l");
1932  leg_t.AddEntry(h_t_fit_mp,"#bar{B}_{s} #rightarrow D_{s}^{+} K^{-} #pi^{-} #pi^{+}","l");
1933  leg_t.AddEntry(h_t_fit_pp,"B_{s} #rightarrow D_{s}^{+} K^{-} #pi^{-} #pi^{+}","l");
1934  leg_t.AddEntry(h_t_fit_0p,"Untagged #rightarrow D_{s}^{+} K^{-} #pi^{-} #pi^{+}","l");
1935  leg_t.Draw();
1936  c->Print(((string)OutputDir+"leg_t.eps").c_str());
1937 
1938  s_Kpipi->SetLineColor(kBlack);
1939  s_Kpipi->DrawNormalized("e1",1);
1940  s_Kpipi_fit->SetLineColor(kBlue);
1941  s_Kpipi_fit->SetLineWidth(3);
1942  s_Kpipi_fit->DrawNormalized("histcsame",1);
1943  c->Print(((string)OutputDir+"s_Kpipi.eps").c_str());
1944 
1945  s_Kpipi->SetLineColor(kBlack);
1946  s_Kpipi->DrawNormalized("e1",1);
1947  s_Kpipi_fit->SetLineColor(kBlack);
1948  s_Kpipi_fit->SetLineWidth(3);
1949  s_Kpipi_fit->DrawNormalized("histcsame",1);
1950  s_Kpipi_fitBs->SetLineColor(kRed);
1951  s_Kpipi_fitBs->SetLineWidth(3);
1952  s_Kpipi_fitBs->SetLineStyle(kDashed);
1953  s_Kpipi_fitBs->DrawNormalized("histcsame",1./(1.+k_fit[0]*k_fit[0])*eff_tag);
1954  s_Kpipi_fitBsbar->SetLineWidth(3);
1955  s_Kpipi_fitBsbar->SetLineColor(kBlue);
1956  s_Kpipi_fitBsbar->SetLineStyle(kDashed);
1957  s_Kpipi_fitBsbar->DrawNormalized("histcsame",k_fit[0]*k_fit[0]/(1.+k_fit[0]*k_fit[0])*eff_tag);
1958  s_Kpipi_fit_notag->SetLineWidth(3);
1959  s_Kpipi_fit_notag->SetLineColor(kMagenta);
1960  s_Kpipi_fit_notag->SetLineStyle(kDashed);
1961  s_Kpipi_fit_notag->DrawNormalized("histcsame",(1.-eff_tag));
1962  c->Print(((string)OutputDir+"s_Kpipi_2.eps").c_str());
1963 
1964 
1965  s_Kpi->SetLineColor(kBlack);
1966  s_Kpi->DrawNormalized("e1",1);
1967  s_Kpi_fit->SetLineColor(kBlue);
1968  s_Kpi_fit->SetLineWidth(3);
1969  s_Kpi_fit->DrawNormalized("histcsame",1);
1970  c->Print(((string)OutputDir+"s_Kpi.eps").c_str());
1971 
1972  s_Kpi->SetLineColor(kBlack);
1973  s_Kpi->DrawNormalized("e1",1);
1974  s_Kpi_fit->SetLineColor(kBlack);
1975  s_Kpi_fit->SetLineWidth(3);
1976  s_Kpi_fit->DrawNormalized("histcsame",1);
1977  s_Kpi_fitBs->SetLineColor(kRed);
1978  s_Kpi_fitBs->SetLineWidth(3);
1979  s_Kpi_fitBs->SetLineStyle(kDashed);
1980  s_Kpi_fitBs->DrawNormalized("histcsame",1./(1.+k_fit[0]*k_fit[0])*eff_tag);
1981  s_Kpi_fitBsbar->SetLineWidth(3);
1982  s_Kpi_fitBsbar->SetLineColor(kBlue);
1983  s_Kpi_fitBsbar->SetLineStyle(kDashed);
1984  s_Kpi_fitBsbar->DrawNormalized("histcsame",k_fit[0]*k_fit[0]/(1.+k_fit[0]*k_fit[0])*eff_tag);
1985  s_Kpi_fit_notag->SetLineWidth(3);
1986  s_Kpi_fit_notag->SetLineColor(kMagenta);
1987  s_Kpi_fit_notag->SetLineStyle(kDashed);
1988  s_Kpi_fit_notag->DrawNormalized("histcsame",(1.-eff_tag));
1989  c->Print(((string)OutputDir+"s_Kpi_2.eps").c_str());
1990 
1991  s_pipi->SetLineColor(kBlack);
1992  s_pipi->DrawNormalized("e1",1);
1993  s_pipi_fit->SetLineColor(kBlue);
1994  s_pipi_fit->SetLineWidth(3);
1995  s_pipi_fit->DrawNormalized("histcsame",1);
1996  c->Print(((string)OutputDir+"s_pipi.eps").c_str());
1997 
1998  s_pipi->SetLineColor(kBlack);
1999  s_pipi->DrawNormalized("e1",1);
2000  s_pipi_fit->SetLineColor(kBlack);
2001  s_pipi_fit->SetLineWidth(3);
2002  s_pipi_fit->DrawNormalized("histcsame",1);
2003  s_pipi_fitBs->SetLineColor(kRed);
2004  s_pipi_fitBs->SetLineWidth(3);
2005  s_pipi_fitBs->SetLineStyle(kDashed);
2006  s_pipi_fitBs->DrawNormalized("histcsame",1./(1.+k_fit[0]*k_fit[0])*eff_tag);
2007  s_pipi_fitBsbar->SetLineWidth(3);
2008  s_pipi_fitBsbar->SetLineColor(kBlue);
2009  s_pipi_fitBsbar->SetLineStyle(kDashed);
2010  s_pipi_fitBsbar->DrawNormalized("histcsame",k_fit[0]*k_fit[0]/(1.+k_fit[0]*k_fit[0])*eff_tag);
2011  s_pipi_fit_notag->SetLineWidth(3);
2012  s_pipi_fit_notag->SetLineColor(kMagenta);
2013  s_pipi_fit_notag->SetLineStyle(kDashed);
2014  s_pipi_fit_notag->DrawNormalized("histcsame",(1.-eff_tag));
2015  c->Print(((string)OutputDir+"s_pipi_2.eps").c_str());
2016 
2017  s_Dspipi->SetLineColor(kBlack);
2018  s_Dspipi->DrawNormalized("e1",1);
2019  s_Dspipi_fit->SetLineColor(kBlue);
2020  s_Dspipi_fit->SetLineWidth(3);
2021  s_Dspipi_fit->DrawNormalized("histcsame",1);
2022  c->Print(((string)OutputDir+"s_Dspipi.eps").c_str());
2023 
2024  s_Dspipi->SetLineColor(kBlack);
2025  s_Dspipi->DrawNormalized("e1",1);
2026  s_Dspipi_fit->SetLineColor(kBlack);
2027  s_Dspipi_fit->SetLineWidth(3);
2028  s_Dspipi_fit->DrawNormalized("histcsame",1);
2029  s_Dspipi_fitBs->SetLineColor(kRed);
2030  s_Dspipi_fitBs->SetLineWidth(3);
2031  s_Dspipi_fitBs->SetLineStyle(kDashed);
2032  s_Dspipi_fitBs->DrawNormalized("histcsame",1./(1.+k_fit[0]*k_fit[0])*eff_tag);
2033  s_Dspipi_fitBsbar->SetLineWidth(3);
2034  s_Dspipi_fitBsbar->SetLineColor(kBlue);
2035  s_Dspipi_fitBsbar->SetLineStyle(kDashed);
2036  s_Dspipi_fitBsbar->DrawNormalized("histcsame",k_fit[0]*k_fit[0]/(1.+k_fit[0]*k_fit[0])*eff_tag);
2037  s_Dspipi_fit_notag->SetLineWidth(3);
2038  s_Dspipi_fit_notag->SetLineColor(kMagenta);
2039  s_Dspipi_fit_notag->SetLineStyle(kDashed);
2040  s_Dspipi_fit_notag->DrawNormalized("histcsame",(1.-eff_tag));
2041  c->Print(((string)OutputDir+"s_Dspipi_2.eps").c_str());
2042 
2043 
2044  s_DsK->SetLineColor(kBlack);
2045  s_DsK->DrawNormalized("e1",1);
2046  s_DsK_fit->SetLineColor(kBlue);
2047  s_DsK_fit->SetLineWidth(3);
2048  s_DsK_fit->DrawNormalized("histcsame",1);
2049  c->Print(((string)OutputDir+"s_DsK.eps").c_str());
2050 
2051  s_DsK->SetLineColor(kBlack);
2052  s_DsK->DrawNormalized("e1",1);
2053  s_DsK_fit->SetLineColor(kBlack);
2054  s_DsK_fit->SetLineWidth(3);
2055  s_DsK_fit->DrawNormalized("histcsame",1);
2056  s_DsK_fitBs->SetLineColor(kRed);
2057  s_DsK_fitBs->SetLineWidth(3);
2058  s_DsK_fitBs->SetLineStyle(kDashed);
2059  s_DsK_fitBs->DrawNormalized("histcsame",1./(1.+k_fit[0]*k_fit[0])*eff_tag);
2060  s_DsK_fitBsbar->SetLineWidth(3);
2061  s_DsK_fitBsbar->SetLineColor(kBlue);
2062  s_DsK_fitBsbar->SetLineStyle(kDashed);
2063  s_DsK_fitBsbar->DrawNormalized("histcsame",k_fit[0]*k_fit[0]/(1.+k_fit[0]*k_fit[0])*eff_tag);
2064  s_DsK_fit_notag->SetLineWidth(3);
2065  s_DsK_fit_notag->SetLineColor(kMagenta);
2066  s_DsK_fit_notag->SetLineStyle(kDashed);
2067  s_DsK_fit_notag->DrawNormalized("histcsame",(1.-eff_tag));
2068  c->Print(((string)OutputDir+"s_DsK_2.eps").c_str());
2069 
2070  s_DsKpi->SetLineColor(kBlack);
2071  s_DsKpi->DrawNormalized("e1",1);
2072  s_DsKpi_fit->SetLineColor(kBlue);
2073  s_DsKpi_fit->SetLineWidth(3);
2074  s_DsKpi_fit->DrawNormalized("histcsame",1);
2075  c->Print(((string)OutputDir+"s_DsKpi.eps").c_str());
2076 
2077  s_Dspi->SetLineColor(kBlack);
2078  s_Dspi->DrawNormalized("e1",1);
2079  s_Dspi_fit->SetLineColor(kBlue);
2080  s_Dspi_fit->SetLineWidth(3);
2081  s_Dspi_fit->DrawNormalized("histcsame",1);
2082  c->Print(((string)OutputDir+"s_Dspi.eps").c_str());
2083 
2084  s_Dspim->SetLineColor(kBlack);
2085  s_Dspim->DrawNormalized("e1",1);
2086  s_Dspim_fit->SetLineColor(kBlue);
2087  s_Dspim_fit->SetLineWidth(3);
2088  s_Dspim_fit->DrawNormalized("histcsame",1);
2089  c->Print(((string)OutputDir+"s_Dspim.eps").c_str());
2090 
2091  }
2092 
2093  if(do2DScan == 1){
2094  cout << "Now doing 2D scan:" << endl;
2095  int scanBins=20;
2096  double scanMin=0, scanMax=360;
2097  double nSigmaZoom = 2;
2098  double scanMinGammaZoom=min(gamma.meanInit(), gamma.mean()) - nSigmaZoom*gamma.err();
2099  double scanMaxGammaZoom=max(gamma.meanInit(), gamma.mean()) + nSigmaZoom*gamma.err();
2100  double scanMinDeltaZoom=min(delta.meanInit(), delta.mean()) - nSigmaZoom*delta.err();
2101  double scanMaxDeltaZoom=max(delta.meanInit(), delta.mean()) + nSigmaZoom*delta.err();
2102  double gammaZoomRange = scanMaxGammaZoom - scanMinGammaZoom;
2103  double deltaZoomRange = scanMaxDeltaZoom - scanMinDeltaZoom;
2104 
2105  TFile* scanFile = new TFile("scan.root", "RECREATE");
2106  TH2D* scanHisto = new TH2D("scan", "; #gamma [deg]; #delta [deg]", scanBins, scanMin, scanMax, scanBins, scanMin, scanMax);
2107  TH2D* scanHistoP = new TH2D("scanP", "; #gamma [deg]; #delta [deg]", scanBins, scanMin, scanMax, scanBins, scanMin, scanMax);
2108  TH2D* scanHistoM = new TH2D("scanM" , "; #gamma [deg]; #delta [deg]", scanBins, scanMin, scanMax, scanBins, scanMin, scanMax);
2109  TH2D* scanZoomHisto = new TH2D("scanZoom", "; #gamma [deg]; #delta [deg]", scanBins, scanMinGammaZoom, scanMaxGammaZoom, scanBins, scanMinDeltaZoom, scanMaxDeltaZoom);
2110 
2111  double scanMinLL=-9999;
2112  double scanMinLLP=-9999;
2113  double scanMinLLM=-9999;
2114  double scanMinLLZ=-9999;
2115 
2116  for(int i=0; i < scanBins; i++){
2117  double gamma_value = ((double)i+0.5)*360/((double)scanBins);
2118  gamma.setCurrentFitVal(gamma_value);
2119  for(int j=0; j < scanBins; j++){
2120  double delta_value = ((double)j+0.5)*360/((double)scanBins);
2121  delta.setCurrentFitVal(delta_value);
2122 
2123  double v = neg2LLSum.getNewVal();
2124  if( (i==0 && j==0) || v < scanMinLL) scanMinLL=v;
2125  scanHisto->Fill(gamma_value, delta_value, v);
2126 
2127  double vP = fcn.getNewVal();
2128  if( (i==0 && j==0) || vP < scanMinLLP) scanMinLLP=vP;
2129  scanHistoP->Fill(gamma_value, delta_value, vP);
2130 
2131  double vM = fcn_bar.getNewVal();
2132  if( (i==0 && j==0) || vM < scanMinLLM) scanMinLLM=vM;
2133  scanHistoM->Fill(gamma_value, delta_value, vM);
2134  }
2135  }
2136  for(int i=0; i < scanBins; i++){
2137  double gamma_value = scanMinGammaZoom + ((double)i+0.5) * gammaZoomRange/((double)scanBins);
2138  gamma.setCurrentFitVal(gamma_value);
2139  for(int j=0; j < scanBins; j++){
2140  double delta_value = scanMinDeltaZoom + ((double)j+0.5) * deltaZoomRange/((double)scanBins);
2141  delta.setCurrentFitVal(delta_value);
2142  double v = neg2LLSum.getNewVal();
2143 
2144  if( (i==0 && j==0) || v < scanMinLLZ) scanMinLLZ=v;
2145 
2146  scanZoomHisto->Fill(gamma_value, delta_value, v);
2147  }
2148  }
2149 
2150  for(int i=0; i < scanBins; i++){
2151  double gamma_value = ((double)i+0.5)*360/((double)scanBins);
2152  for(int j=0; j < scanBins; j++){
2153  double delta_value = ((double)j+0.5)*360/((double)scanBins);
2154  scanHisto->Fill(gamma_value, delta_value, -scanMinLL);
2155  scanHistoP->Fill(gamma_value, delta_value, -scanMinLLP);
2156  scanHistoM->Fill(gamma_value, delta_value, -scanMinLLM);
2157  }
2158  }
2159  for(int i=0; i < scanBins; i++){
2160  double gamma_value = scanMinGammaZoom + ((double)i+0.5) * gammaZoomRange/((double)scanBins);
2161  for(int j=0; j < scanBins; j++){
2162  double delta_value = scanMinDeltaZoom + ((double)j+0.5) * deltaZoomRange/((double)scanBins);
2163  scanZoomHisto->Fill(gamma_value, delta_value, -scanMinLLZ);
2164  }
2165  }
2166  scanFile->cd();
2167  scanHisto->Write();
2168  scanHistoP->Write();
2169  scanHistoM->Write();
2170  scanZoomHisto->Write();
2171  scanFile->Close();
2172 
2173  cout<< "done 2-D scan" << endl;
2174  }
2175 
2176  return 0;
2177 }
2178 
2179 // Coherence factor studies
2180 void calculateCoherence(int step=0){
2181 
2182  TRandom3 ranLux;
2183  NamedParameter<int> RandomSeed("RandomSeed", 0);
2184  int seed = RandomSeed + step;
2185  ranLux.SetSeed((int)seed);
2186  gRandom = &ranLux;
2187 
2188  FitAmplitude::AutogenerateFitFile();
2189 
2190  NamedParameter<string> IntegratorEventFile("IntegratorEventFile", (std::string) "SignalIntegrationEvents.root", (char*) 0);
2191  NamedParameter<int> EventPattern("Event Pattern", 521, 321, 211, -211, 443);
2192  DalitzEventPattern pat(EventPattern.getVector());
2193  cout << " got event pattern: " << pat << endl;
2194 
2195  NamedParameter<string> OutputDir("OutputDir", (std::string) "", (char*) 0);
2196  TFile* paraFile = new TFile(((string)OutputDir+"coherence_"+anythingToString((int)seed)+".root").c_str(), "RECREATE");
2197  paraFile->cd();
2198  TTree* tree = new TTree("coherence","coherence");
2199  double r_val,k,delta_val;
2200  double x,y;
2201  TBranch* br_r = tree->Branch( "r", &r_val, "r_val/D" );
2202  TBranch* br_k = tree->Branch( "k", &k, "k/D" );
2203  TBranch* br_delta = tree->Branch( "delta", &delta_val, "delta_val/D" );
2204 
2205  double r_Ks_val,k_Ks,delta_Ks_val;
2206  double r_rho_val,k_rho,delta_rho_val;
2207  double r_Ks_rho_val,k_Ks_rho,delta_Ks_rho_val;
2208  double r_K1_val,k_K1,delta_K1_val;
2209 
2210  TBranch* br_r_Ks = tree->Branch( "r_Ks", &r_Ks_val, "r_Ks_val/D" );
2211  TBranch* br_k_Ks = tree->Branch( "k_Ks", &k_Ks, "k_Ks/D" );
2212  TBranch* br_delta_Ks = tree->Branch( "delta_Ks", &delta_Ks_val, "delta_Ks_val/D" );
2213  TBranch* br_r_Ks_rho = tree->Branch( "r_Ks_rho", &r_Ks_rho_val, "r_Ks_rho_val/D" );
2214  TBranch* br_k_Ks_rho = tree->Branch( "k_Ks_rho", &k_Ks_rho, "k_Ks_rho/D" );
2215  TBranch* br_delta_Ks_rho = tree->Branch( "delta_Ks_rho", &delta_Ks_rho_val, "delta_Ks_rho_val/D" );
2216  TBranch* br_r_rho = tree->Branch( "r_rho", &r_rho_val, "r_rho_val/D" );
2217  TBranch* br_k_rho = tree->Branch( "k_rho", &k_rho, "k_rho/D" );
2218  TBranch* br_delta_rho = tree->Branch( "delta_rho", &delta_rho_val, "delta_rho_val/D" );
2219  TBranch* br_r_K1 = tree->Branch( "r_K1", &r_K1_val, "r_K1_val/D" );
2220  TBranch* br_k_K1 = tree->Branch( "k_K1", &k_K1, "k_K1/D" );
2221  TBranch* br_delta_K1 = tree->Branch( "delta_K1", &delta_K1_val, "delta_K1_val/D" );
2222 
2223  double r_Ks_val_2,k_Ks_2,delta_Ks_val_2;
2224  double r_rho_val_2,k_rho_2,delta_rho_val_2;
2225  double r_Ks_rho_val_2,k_Ks_rho_2,delta_Ks_rho_val_2;
2226  double r_K1_val_2,k_K1_2,delta_K1_val_2;
2227 
2228  TBranch* br_r_Ks_2 = tree->Branch( "r_Ks_2", &r_Ks_val_2, "r_Ks_val_2/D" );
2229  TBranch* br_k_Ks_2 = tree->Branch( "k_Ks_2", &k_Ks_2, "k_Ks_2/D" );
2230  TBranch* br_delta_Ks_2 = tree->Branch( "delta_Ks_2", &delta_Ks_val_2, "delta_Ks_val_2/D" );
2231  TBranch* br_r_Ks_rho_2 = tree->Branch( "r_Ks_rho_2", &r_Ks_rho_val_2, "r_Ks_rho_val_2/D" );
2232  TBranch* br_k_Ks_rho_2 = tree->Branch( "k_Ks_rho_2", &k_Ks_rho_2, "k_Ks_rho_2/D" );
2233  TBranch* br_delta_Ks_rho_2 = tree->Branch( "delta_Ks_rho_2", &delta_Ks_rho_val_2, "delta_Ks_rho_val_2/D" );
2234  TBranch* br_r_rho_2 = tree->Branch( "r_rho_2", &r_rho_val_2, "r_rho_val_2/D" );
2235  TBranch* br_k_rho_2 = tree->Branch( "k_rho_2", &k_rho_2, "k_rho_2/D" );
2236  TBranch* br_delta_rho_2 = tree->Branch( "delta_rho_2", &delta_rho_val_2, "delta_rho_val_2/D" );
2237  TBranch* br_r_K1_2 = tree->Branch( "r_K1_2", &r_K1_val_2, "r_K1_val_2/D" );
2238  TBranch* br_k_K1_2 = tree->Branch( "k_K1_2", &k_K1_2, "k_K1_2/D" );
2239  TBranch* br_delta_K1_2 = tree->Branch( "delta_K1_2", &delta_K1_val_2, "delta_K1_val_2/D" );
2240 
2241  TBranch* br_x = tree->Branch( "x", &x, "x/D" );
2242  TBranch* br_y = tree->Branch( "y", &y, "y/D" );
2243 
2244  FitParameter r("r");
2245  FitParameter delta("delta");
2246 
2247  DalitzEventList eventListPhsp,eventList;
2248  eventListPhsp.generatePhaseSpaceEvents(100000,pat);
2249 
2250  /*
2251  FitAmpIncoherentSum fasInco((DalitzEventPattern)pat);
2252  fasInco.normalizeAmps(eventListPhsp);
2253  SignalGenerator sg(pat,&fasInco);
2254  sg.FillEventList(eventList, 200000);
2255  */
2256 
2257  TFile *FileMC = TFile::Open(((string) IntegratorEventFile).c_str());
2258  TTree* treeMC = dynamic_cast<TTree*>(FileMC->Get("DalitzEventList"));
2259  eventList.fromNtuple(treeMC,1);
2260  FileMC->Close();
2261 
2262  vector<int> s234;
2263  s234.push_back(2);
2264  s234.push_back(3);
2265  s234.push_back(4);
2266 
2267  DalitzEventList eventList_Ks,eventList_rho,eventList_Ks_rho,eventList_K1;
2268  DalitzEventList eventList_Ks_2,eventList_rho_2,eventList_Ks_rho_2,eventList_K1_2;
2269 
2270  for(unsigned int i=0; i<eventList.size(); i++){
2271  if(abs(sqrt(eventList[i].s(2,4)/(GeV*GeV))-0.89166)<0.0474) eventList_Ks.Add(eventList[i]);
2272  if(abs(sqrt(eventList[i].sij(s234)/(GeV*GeV))-1.272)<0.09) eventList_K1.Add(eventList[i]);
2273  if(abs(sqrt(eventList[i].s(3,4)/(GeV*GeV))-0.7755)<0.1494) eventList_rho.Add(eventList[i]);
2274  if(abs(sqrt(eventList[i].s(2,4)/(GeV*GeV))-0.89166)<0.0474 || abs(sqrt(eventList[i].s(3,4)/(GeV*GeV))-0.7755)<0.1494) eventList_Ks_rho.Add(eventList[i]);
2275 
2276  if(abs(sqrt(eventList[i].s(2,4)/(GeV*GeV))-0.89166)<0.0474/2.) eventList_Ks_2.Add(eventList[i]);
2277  if(abs(sqrt(eventList[i].sij(s234)/(GeV*GeV))-1.272)<0.09/2.) eventList_K1_2.Add(eventList[i]);
2278  if(abs(sqrt(eventList[i].s(3,4)/(GeV*GeV))-0.7755)<0.1494/2.) eventList_rho_2.Add(eventList[i]);
2279  if(abs(sqrt(eventList[i].s(2,4)/(GeV*GeV))-0.89166)<0.0474/2. || abs(sqrt(eventList[i].s(3,4)/(GeV*GeV))-0.7755)<0.1494/2.) eventList_Ks_rho_2.Add(eventList[i]);
2280  }
2281 
2282  FitAmpSum fas_tmp((DalitzEventPattern)pat);
2283  fas_tmp.getVal(eventListPhsp[0]);
2284  fas_tmp.normalizeAmps(eventListPhsp);
2285 
2286  counted_ptr<FitAmpList> List_1 = fas_tmp.GetCloneOfSubsetSameFitParameters("K(1)(1270)+");
2287  FitAmpSum fas(*List_1);
2288  FitAmpSum fasCC(*List_1);
2289  fasCC.CPConjugateSameFitParameters();
2290  fasCC.CConjugateFinalStateSameFitParameters();
2291  x = ranLux.Uniform(-1,1);
2292  y = ranLux.Uniform(-180,180);
2293  FitParameter r_K1_re("r_K1_Re",2,x,0.01);
2294  FitParameter r_K1_im("r_K1_Im",2,y,0.01);
2295  counted_ptr<IReturnComplex> r_K1_plus = new CPV_amp_polar(r_K1_re,r_K1_im,1);
2296  counted_ptr<IReturnComplex> r_K1_minus = new CPV_amp_polar(r_K1_re,r_K1_im,-1);
2297  fas.multiply(r_K1_plus);
2298  fasCC.multiply(r_K1_minus);
2299 
2300  AddScaledAmpsToList(fas_tmp, fas, fasCC, "K(1)(1400)+", r_K1_plus, r_K1_minus );
2301 
2302  x = ranLux.Uniform(-1,1);
2303  y = ranLux.Uniform(-180,180);
2304  FitParameter r_2_re("r_2_Re",2,x,0.01);
2305  FitParameter r_2_im("r_2_Im",2,y,0.01);
2306  counted_ptr<IReturnComplex> r_2_plus = new CPV_amp_polar(r_2_re,r_2_im,1);
2307  counted_ptr<IReturnComplex> r_2_minus = new CPV_amp_polar(r_2_re,r_2_im,-1);
2308  AddScaledAmpsToList(fas_tmp, fas, fasCC, "K*(1410)+", r_2_plus, r_2_minus );
2309 
2310  x = ranLux.Uniform(-1,1);
2311  y = ranLux.Uniform(-180,180);
2312  FitParameter r_3_re("r_3_Re",2,x,0.01);
2313  FitParameter r_3_im("r_3_Im",2,y,0.01);
2314  counted_ptr<IReturnComplex> r_3_plus = new CPV_amp_polar(r_3_re,r_3_im,1);
2315  counted_ptr<IReturnComplex> r_3_minus = new CPV_amp_polar(r_3_re,r_3_im,-1);
2316  AddScaledAmpsToList(fas_tmp, fas, fasCC, "NonResV0(->Ds-,pi+),K*(892)0(->K+,pi-)", r_3_plus, r_3_minus );
2317  AddScaledAmpsToList(fas_tmp, fas, fasCC, "NonResS0(->Ds-,pi+),K*(892)0(->K+,pi-)", r_3_plus, r_3_minus );
2318  AddScaledAmpsToList(fas_tmp, fas, fasCC, "NonResV0(->Ds-,K+),sigma10(->pi+,pi-)", r_3_plus, r_3_minus );
2319 
2320  x = ranLux.Uniform(-1,1);
2321  y = ranLux.Uniform(-180,180);
2322  FitParameter r_4_re("r_4_Re",2,x,0.01);
2323  FitParameter r_4_im("r_4_Im",2,y,0.01);
2324  counted_ptr<IReturnComplex> r_4_plus = new CPV_amp_polar(r_4_re,r_4_im,1);
2325  counted_ptr<IReturnComplex> r_4_minus = new CPV_amp_polar(r_4_re,r_4_im,-1);
2326  AddScaledAmpsToList(fas_tmp, fas, fasCC, "NonResA0(->sigma10(->pi+,pi-),Ds-)", r_4_plus, r_4_minus );
2327  AddScaledAmpsToList(fas_tmp, fas, fasCC, "NonResA0(->rho(770)0(->pi+,pi-),Ds-),K+", r_4_plus, r_4_minus );
2328 
2329 
2330  fas.getVal(eventListPhsp[0]);
2331  fasCC.getVal(eventListPhsp[0]);
2332 
2333  cout << " Have " << eventList.size() << " integration events " << endl;
2334  cout << " Have " << eventList_Ks.size() << " Ks integration events " << endl;
2335  cout << " Have " << eventList_rho.size() << " rho integration events " << endl;
2336  cout << " Have " << eventList_Ks_rho.size() << " Ks/rho integration events " << endl;
2337  cout << " Have " << eventList_K1.size() << " K1 integration events " << endl;
2338 
2339  cout << " Have " << eventList_Ks_2.size() << " Ks_2 integration events " << endl;
2340  cout << " Have " << eventList_rho_2.size() << " rho_2 integration events " << endl;
2341  cout << " Have " << eventList_Ks_rho_2.size() << " Ks/rho_2 integration events " << endl;
2342  cout << " Have " << eventList_K1_2.size() << " K1_2 integration events " << endl;
2343 
2344  //vector<double> kappa = coherenceFactor(fas,fasCC,(double)r, (double)delta,eventListPhsp);
2345  vector<double> kappa = coherenceFactor(fas,fasCC,(double)r, (double)delta,eventList);
2346 
2347  vector<double> kappa_Ks = coherenceFactor(fas,fasCC,(double)r, (double)delta,eventList_Ks);
2348  vector<double> kappa_rho = coherenceFactor(fas,fasCC,(double)r, (double)delta,eventList_rho);
2349  vector<double> kappa_Ks_rho = coherenceFactor(fas,fasCC,(double)r, (double)delta,eventList_Ks_rho);
2350  vector<double> kappa_K1 = coherenceFactor(fas,fasCC,(double)r, (double)delta,eventList_K1);
2351 
2352  vector<double> kappa_Ks_2 = coherenceFactor(fas,fasCC,(double)r, (double)delta,eventList_Ks_2);
2353  vector<double> kappa_rho_2 = coherenceFactor(fas,fasCC,(double)r, (double)delta,eventList_rho_2);
2354  vector<double> kappa_Ks_rho_2 = coherenceFactor(fas,fasCC,(double)r, (double)delta,eventList_Ks_rho_2);
2355  vector<double> kappa_K1_2 = coherenceFactor(fas,fasCC,(double)r, (double)delta,eventList_K1_2);
2356 
2357  r_val = kappa[0];
2358  k = kappa[1];
2359  delta_val = kappa[2];
2360 
2361  r_Ks_val = kappa_Ks[0];
2362  k_Ks = kappa_Ks[1];
2363  delta_Ks_val = kappa_Ks[2];
2364 
2365  r_Ks_rho_val = kappa_Ks_rho[0];
2366  k_Ks_rho = kappa_Ks_rho[1];
2367  delta_Ks_rho_val = kappa_Ks_rho[2];
2368 
2369  r_rho_val = kappa_rho[0];
2370  k_rho = kappa_rho[1];
2371  delta_rho_val = kappa_rho[2];
2372 
2373  r_K1_val = kappa_K1[0];
2374  k_K1 = kappa_K1[1];
2375  delta_K1_val = kappa_K1[2];
2376 
2377  r_Ks_val_2 = kappa_Ks_2[0];
2378  k_Ks_2 = kappa_Ks_2[1];
2379  delta_Ks_val_2 = kappa_Ks_2[2];
2380 
2381  r_Ks_rho_val_2 = kappa_Ks_rho_2[0];
2382  k_Ks_rho_2 = kappa_Ks_rho_2[1];
2383  delta_Ks_rho_val_2 = kappa_Ks_rho_2[2];
2384 
2385  r_rho_val_2 = kappa_rho_2[0];
2386  k_rho_2 = kappa_rho_2[1];
2387  delta_rho_val_2 = kappa_rho_2[2];
2388 
2389  r_K1_val_2 = kappa_K1_2[0];
2390  k_K1_2 = kappa_K1_2[1];
2391  delta_K1_val_2 = kappa_K1_2[2];
2392 
2393  tree->Fill();
2394  paraFile->cd();
2395  tree->SetDirectory(paraFile);
2396  tree->Write();
2397  paraFile->Close();
2398  delete paraFile;
2399 
2400  return;
2401 }
2402 
2403 void coherence_plot(){
2404 
2405  NamedParameter<string> OutputDir("OutputDir", (std::string) "", (char*) 0);
2406 
2407  TFile *File = TFile::Open(((string)OutputDir+"coherence.root").c_str());
2408  TTree* tree;
2409  tree=dynamic_cast<TTree*>(File->Get("coherence"));
2410  Double_t k,k_Ks,k_rho,k_K1;
2411  tree->SetBranchAddress("k",&k) ;
2412  tree->SetBranchAddress("k_Ks",&k_Ks) ;
2413  tree->SetBranchAddress("k_rho",&k_rho) ;
2414  tree->SetBranchAddress("k_K1",&k_K1) ;
2415 
2416  TH1D* h_k = new TH1D("",";#kappa; Number of experiments ",50,0,1);
2417 
2418  TH1D* h_k_Ks = new TH1D("","K^{*0}(892) region ;#kappa; Number of experiments ",50,0,1);
2419  TH1D* h_k_rho = new TH1D("","#rho^{0}(770) region ;#kappa; Number of experiments ",50,0,1);
2420  TH1D* h_k_K1 = new TH1D("","K_{1}(1270) region ;#kappa; Number of experiments ",50,0,1);
2421 
2422  for(int i = 0; i<tree->GetEntries(); i++){
2423  tree->GetEntry(i);
2424  h_k->Fill(k);
2425  h_k_Ks->Fill(k_Ks);
2426  h_k_rho->Fill(k_rho);
2427  h_k_K1->Fill(k_K1);
2428  }
2429 
2430  TCanvas *c = new TCanvas();
2431  h_k->SetLineColor(kBlue);
2432  h_k->Draw("hist");
2433  c->Print(((string)OutputDir+"k.eps").c_str());
2434 
2435  h_k_Ks->SetLineColor(kBlue);
2436  h_k_Ks->Draw("hist");
2437  c->Print(((string)OutputDir+"k_Ks.eps").c_str());
2438 
2439  h_k_rho->SetLineColor(kBlue);
2440  h_k_rho->Draw("hist");
2441  c->Print(((string)OutputDir+"k_rho.eps").c_str());
2442 
2443  h_k_K1->SetLineColor(kBlue);
2444  h_k_K1->Draw("hist");
2445  c->Print(((string)OutputDir+"k_K1.eps").c_str());
2446 
2447  DalitzEventList eventList;
2448  TFile *_InputFile = TFile::Open("/auto/data/dargent/Bs2DsKpipi/MINT2/MINT_data_3sigma.root");
2449  TTree* in_tree;
2450  in_tree=dynamic_cast<TTree*>(_InputFile->Get("DalitzEventList"));
2451  eventList.fromNtuple(in_tree,1);
2452  _InputFile->Close();
2453 
2454  cout << " I've got " << eventList.size() << " events." << endl;
2455 
2456  DalitzEventList eventList_Ks,eventList_rho,eventList_Ks_rho,eventList_K1;
2457  DalitzEventList eventList_Ks_2,eventList_rho_2,eventList_Ks_rho_2,eventList_K1_2;
2458  vector<int> s234;
2459  s234.push_back(2);
2460  s234.push_back(3);
2461  s234.push_back(4);
2462 
2463  for(unsigned int i=0; i<eventList.size(); i++){
2464  if(abs(sqrt(eventList[i].s(2,4)/(GeV*GeV))-0.89166)<0.0474) eventList_Ks.Add(eventList[i]);
2465  if(abs(sqrt(eventList[i].sij(s234)/(GeV*GeV))-1.272)<0.09) eventList_K1.Add(eventList[i]);
2466  if(abs(sqrt(eventList[i].s(3,4)/(GeV*GeV))-0.7755)<0.1494) eventList_rho.Add(eventList[i]);
2467  if(abs(sqrt(eventList[i].s(2,4)/(GeV*GeV))-0.89166)<0.0474 || abs(sqrt(eventList[i].s(3,4)/(GeV*GeV))-0.7755)<0.1494) eventList_Ks_rho.Add(eventList[i]);
2468 
2469  if(abs(sqrt(eventList[i].s(2,4)/(GeV*GeV))-0.89166)<0.0474/2.) eventList_Ks_2.Add(eventList[i]);
2470  if(abs(sqrt(eventList[i].sij(s234)/(GeV*GeV))-1.272)<0.09/2.) eventList_K1_2.Add(eventList[i]);
2471  if(abs(sqrt(eventList[i].s(3,4)/(GeV*GeV))-0.7755)<0.1494/2.) eventList_rho_2.Add(eventList[i]);
2472  if(abs(sqrt(eventList[i].s(2,4)/(GeV*GeV))-0.89166)<0.0474/2. || abs(sqrt(eventList[i].s(3,4)/(GeV*GeV))-0.7755)<0.1494/2.) eventList_Ks_rho_2.Add(eventList[i]);
2473  }
2474 
2475  cout << " Have " << eventList_Ks.size()/(double)eventList.size() << " Ks integration events " << endl;
2476  cout << " Have " << eventList_rho.size()/(double)eventList.size() << " rho integration events " << endl;
2477  cout << " Have " << eventList_Ks_rho.size()/(double)eventList.size() << " Ks/rho integration events " << endl;
2478  cout << " Have " << eventList_K1.size()/(double)eventList.size() << " K1 integration events " << endl;
2479 
2480  cout << " Have " << eventList_Ks_2.size()/(double)eventList.size() << " Ks_2 integration events " << endl;
2481  cout << " Have " << eventList_rho_2.size()/(double)eventList.size() << " rho_2 integration events " << endl;
2482  cout << " Have " << eventList_Ks_rho_2.size()/(double)eventList.size() << " Ks/rho_2 integration events " << endl;
2483  cout << " Have " << eventList_K1_2.size()/(double)eventList.size() << " K1_2 integration events " << endl;
2484 
2485  return;
2486 }
2487 
2488 void makeIntegratorFile(){
2489 
2490  FitAmplitude::AutogenerateFitFile();
2491 
2492  NamedParameter<string> IntegratorEventFile("IntegratorEventFile", (std::string) "SignalIntegrationEvents.root", (char*) 0);
2493  NamedParameter<int> EventPattern("Event Pattern", 521, 321, 211, -211, 443);
2494  DalitzEventPattern pat(EventPattern.getVector());
2495  cout << " got event pattern: " << pat << endl;
2496 
2497  NamedParameter<int> IntegratorEvents("IntegratorEvents", 300000);
2498 
2499  DalitzEventList eventListPhsp,eventList,eventList_cut;
2500 
2501  eventListPhsp.generatePhaseSpaceEvents(100000,pat);
2502 
2503  FitAmpIncoherentSum fas((DalitzEventPattern)pat);
2504  fas.print();
2505  fas.getVal(eventListPhsp[0]);
2506  fas.normalizeAmps(eventListPhsp);
2507 
2508  SignalGenerator sg(pat,&fas);
2509 
2510  sg.FillEventList(eventList, IntegratorEvents);
2511  vector<int> s234;
2512  s234.push_back(2);
2513  s234.push_back(3);
2514  s234.push_back(4);
2515 
2516  for(int i = 0; i < eventList.size(); i++){
2517  if(sqrt(eventList[i].sij(s234)/(GeV*GeV)) < 1.95 && sqrt(eventList[i].s(2,4)/(GeV*GeV)) < 1.2 && sqrt(eventList[i].s(3,4)/(GeV*GeV)) < 1.2)eventList_cut.Add(eventList[i]);
2518  }
2519 
2520  cout << "Generated " << eventList_cut.size() << " events inside selected phasespace region" << endl;
2521 
2522  eventList_cut.saveAsNtuple(IntegratorEventFile);
2523  return;
2524 }
2525 
2526 
2527 int main(int argc, char** argv){
2528 
2529  time_t startTime = time(0);
2530 
2531  TH1::SetDefaultSumw2();
2532  TH2::SetDefaultSumw2();
2533  gROOT->ProcessLine(".x ../lhcbStyle.C");
2534 
2535  NamedParameter<int> EvtGenTest("EvtGenTest", 0);
2536  NamedParameter<string> IntegratorEventFile("IntegratorEventFile", (std::string) "SignalIntegrationEvents.root", (char*) 0);
2537 
2538  if(!EvtGenTest) if(! std::ifstream(((string)IntegratorEventFile).c_str()).good()) makeIntegratorFile();
2539 
2540  ampFit(atoi(argv[1]));
2541 
2542  cout << "==============================================" << endl;
2543  cout << " Done. " << " Total time since start " << (time(0) - startTime)/60.0 << " min." << endl;
2544  cout << "==============================================" << endl;
2545 
2546  return 0;
2547 }
2548 //
RooGaussEfficiencyModel::evaluate
virtual Double_t evaluate(Int_t basisCodeInt, Double_t tau, Double_t omega, Double_t dGamma) const
Definition: RooGaussEfficiencyModel.cpp:245
AmpsPdfFlexiFastCPV_time_integrated::getVal_withPs
virtual double getVal_withPs(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:1084
TimePdf::un_normalised
double un_normalised(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:1143
MINT::EventList::Add
virtual bool Add(const EVENT_TYPE &evt)
Definition: EventList.h:63
CPV_amp_polar::_delta
FitParameter & _delta
Definition: BsDsKpipi_TD_ampFit.cpp:187
AmpsPdfFlexiFastCPV_mod::_gamma
FitParameter & _gamma
Definition: BsDsKpipi_TD_ampFit.cpp:632
AmpsPdfFlexiFastCPV_mod::_r_q
RooCategory * _r_q
Definition: BsDsKpipi_TD_ampFit.cpp:647
FitAmpSum::GetCloneOfSubsetSameFitParameters
virtual MINT::counted_ptr< FitAmpListBase > GetCloneOfSubsetSameFitParameters(std::string name) const
Definition: FitAmpSum.cpp:147
AddScaledAmpsToList
void AddScaledAmpsToList(FitAmpSum &fas_tmp, FitAmpSum &fas, FitAmpSum &fasCC, std::string name, counted_ptr< IReturnComplex > &r_plus, counted_ptr< IReturnComplex > &r_minus)
Definition: BsDsKpipi_TD_ampFit.cpp:199
AmpsPdfFlexiFastCPV_time_integrated::getVal_noPs
virtual double getVal_noPs(IDalitzEvent *evt)
Definition: BsDsKpipi_TD_ampFit.cpp:1095
AmpsPdfFlexiFastCPV_mod::_eff_tag
FitParameter & _eff_tag
Definition: BsDsKpipi_TD_ampFit.cpp:637
AmpsPdfFlexiFastCPV_mod::_pdf_sigma_t
RooAbsPdf * _pdf_sigma_t
Definition: BsDsKpipi_TD_ampFit.cpp:653
FullAmpsPdfFlexiFastCPV::getVal_withPs
virtual double getVal_withPs(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:461
AmpsPdfFlexiFast::getAmpSum
IFastAmplitudeIntegrable * getAmpSum()
Definition: BsDsKpipi_TD_ampFit.cpp:148
TimePdf_mod::getVal_noPs
virtual double getVal_noPs(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:1263
FitAmplitude::AutogenerateFitFile
static void AutogenerateFitFile(const std::string &fname="protoFitAmplitudeFile.txt", const DalitzEventPattern &pat=DalitzEventPattern::NoPattern)
Definition: FitAmplitude.cpp:50
AmpsPdfFlexiFastCPV_time_integrated::_ampsSum
AmpsPdfFlexiFast * _ampsSum
Definition: BsDsKpipi_TD_ampFit.cpp:999
AmpsPdfFlexiFastCPV::_intAAbar
complex< double > _intAAbar
Definition: BsDsKpipi_TD_ampFit.cpp:891
MINT::Neg2LL::getNewVal
virtual double getNewVal()
Definition: Neg2LL.h:327
AmpsPdfFlexiFastCPV_time_integrated::_w
FitParameter & _w
Definition: BsDsKpipi_TD_ampFit.cpp:1009
TimePdf_mod::getVal_withPs
virtual double getVal_withPs(IDalitzEvent *evt)
Definition: BsDsKpipi_TD_ampFit.cpp:1269
AmpsPdfFlexiFastCPV_time_integrated::AmpsPdfFlexiFastCPV_time_integrated
AmpsPdfFlexiFastCPV_time_integrated(AmpsPdfFlexiFast *amps1, AmpsPdfFlexiFast *amps2, AmpsPdfFlexiFast *ampsSum, MINT::FitParameter &r, MINT::FitParameter &delta, MINT::FitParameter &gamma, MINT::FitParameter &tau, MINT::FitParameter &dGamma, MINT::FitParameter &dm, MINT::FitParameter &eff_tag, MINT::FitParameter &w)
Definition: BsDsKpipi_TD_ampFit.cpp:1110
AmpsPdfFlexiFastCPV_mod::_dGamma
FitParameter & _dGamma
Definition: BsDsKpipi_TD_ampFit.cpp:635
AmpsPdfFlexiFastCPV_mod::_r
FitParameter & _r
Definition: BsDsKpipi_TD_ampFit.cpp:630
FullAmpsPdfFlexiFastCPV::getSmearedTime
void getSmearedTime(double &t, double &dt, TRandom3 &r)
Definition: BsDsKpipi_TD_ampFit.cpp:348
CPV_amp::_sign
int _sign
Definition: BsDsKpipi_TD_ampFit.cpp:174
MINT::PolymorphVector::push_back
void push_back(const T &c)
Definition: PolymorphVector.h:50
AmpsPdfFlexiFastCPV::un_normalised_noPs
double un_normalised_noPs(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:916
AmpsPdfFlexiFastCPV::getVal_noPs
virtual double getVal_noPs(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:961
CPV_amp_polar::ComplexVal
complex< double > ComplexVal()
Definition: BsDsKpipi_TD_ampFit.cpp:193
AmpsPdfFlexiFastCPV_time_integrated::_eff_tag
FitParameter & _eff_tag
Definition: BsDsKpipi_TD_ampFit.cpp:1008
TimePdf_mod::_dGamma
FitParameter & _dGamma
Definition: BsDsKpipi_TD_ampFit.cpp:1211
MINT::Minimisable::getNewVal
double getNewVal()
Definition: Minimisable.h:29
expBasis
Definition: BsDsKpipi_TD_ampFit.cpp:244
FitAmpListBase::CPConjugateSameFitParameters
virtual bool CPConjugateSameFitParameters()
Definition: FitAmpListBase.cpp:146
SignalGenerator::newEvent
virtual MINT::counted_ptr< IDalitzEvent > newEvent()
Definition: SignalGenerator.cpp:130
AmpsPdfFlexiFastCPV_time_integrated::doFinalStatsAndSaveForAmp12
void doFinalStatsAndSaveForAmp12(MINT::Minimiser *min=0, const std::string &fname="FitAmpResults", const std::string &fnameROOT="fitFractions")
Definition: BsDsKpipi_TD_ampFit.cpp:1103
AmpsPdfFlexiFastCPV_mod::un_normalised_noPs
double un_normalised_noPs(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:690
TimePdf_mod::_gamma
FitParameter & _gamma
Definition: BsDsKpipi_TD_ampFit.cpp:1207
DecRateCoeff_Bd::analyticalIntegral
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Definition: DecRateCoeff_Bd.cpp:177
MINT::FitParameter::meanInit
double meanInit() const
Definition: FitParameter.cpp:244
pi
static const double pi
Definition: CLHEPPhysicalConstants.h:44
FitAmpSum::GetCloneSameFitParameters
virtual MINT::counted_ptr< FitAmpListBase > GetCloneSameFitParameters() const
Definition: FitAmpSum.cpp:116
AmpsPdfFlexiFastCPV_mod::_intAAbar
complex< double > _intAAbar
Definition: BsDsKpipi_TD_ampFit.cpp:642
FullAmpsPdfFlexiFastCPV::getVal
virtual double getVal(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:393
AmpsPdfFlexiFastCPV_time_integrated::getVal
virtual double getVal(IDalitzEvent *evt)
Definition: BsDsKpipi_TD_ampFit.cpp:1087
FullAmpsPdfFlexiFastCPV::_min_TAU
NamedParameter< double > _min_TAU
Definition: BsDsKpipi_TD_ampFit.cpp:282
TimePdf_mod::TimePdf_mod
TimePdf_mod(MINT::FitParameter &r, MINT::FitParameter &delta, MINT::FitParameter &gamma, MINT::FitParameter &k, MINT::FitParameter &tau, MINT::FitParameter &dGamma, MINT::FitParameter &dm, MINT::FitParameter &eff_tag, MINT::FitParameter &w)
Definition: BsDsKpipi_TD_ampFit.cpp:1278
TimePdf_mod::getVal
virtual double getVal(IDalitzEvent *evt)
Definition: BsDsKpipi_TD_ampFit.cpp:1265
AmpsPdfFlexiFastCPV_mod::getVal_withPs
virtual double getVal_withPs(IDalitzEvent *evt)
Definition: BsDsKpipi_TD_ampFit.cpp:754
s
static const double s
Definition: CLHEPSystemOfUnits.h:126
TimePdf_mod::_dm
FitParameter & _dm
Definition: BsDsKpipi_TD_ampFit.cpp:1212
DecRateCoeff_Bd::evaluate
Double_t evaluate() const
Definition: DecRateCoeff_Bd.cpp:146
AmpsPdfFlexiFastCPV_mod::_amps1
AmpsPdfFlexiFast * _amps1
Definition: BsDsKpipi_TD_ampFit.cpp:626
AmpsPdfFlexiFastCPV_mod::_max_TAU
NamedParameter< double > _max_TAU
Definition: BsDsKpipi_TD_ampFit.cpp:660
TimePdf_mod::_tau
FitParameter & _tau
Definition: BsDsKpipi_TD_ampFit.cpp:1210
DalitzEvent::setValueInVector
virtual void setValueInVector(unsigned int i, double value)
Definition: DalitzEvent.cpp:415
AmpsPdfFlexiFastCPV_mod::_dm
FitParameter & _dm
Definition: BsDsKpipi_TD_ampFit.cpp:636
FullAmpsPdfFlexiFastCPV::getVal_noPs
virtual double getVal_noPs(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:462
TimePdf::_D
FitParameter & _D
Definition: BsDsKpipi_TD_ampFit.cpp:1123
coherence_plot
void coherence_plot()
Definition: BsDsKpipi_TD_ampFit.cpp:2403
CPV_amp_polar::CPV_amp_polar
CPV_amp_polar(FitParameter &r, FitParameter &delta, int sign)
Definition: BsDsKpipi_TD_ampFit.cpp:190
TimePdf::_dm
FitParameter & _dm
Definition: BsDsKpipi_TD_ampFit.cpp:1131
coherenceFactor
std::vector< double > coherenceFactor(FitAmpSum &fas, FitAmpSum &fas_bar, double r, double delta, DalitzEventList &eventList)
Definition: BsDsKpipi_TD_ampFit.cpp:212
AmpsPdfFlexiFastCPV_mod::_pdf_omega_os
RooAbsPdf * _pdf_omega_os
Definition: BsDsKpipi_TD_ampFit.cpp:654
FullAmpsPdfFlexiFastCPV::un_normalised_noPs
double un_normalised_noPs(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:357
CPV_amp
Definition: BsDsKpipi_TD_ampFit.cpp:171
FitAmpSum::getVal
virtual std::complex< double > getVal(IDalitzEvent &evt)
Definition: FitAmpSum.cpp:182
CPV_amp::CPV_amp
CPV_amp(FitParameter &re, FitParameter &im, int sign)
Definition: BsDsKpipi_TD_ampFit.cpp:176
AmpsPdfFlexiFastCPV_time_integrated::_amps1
AmpsPdfFlexiFast * _amps1
Definition: BsDsKpipi_TD_ampFit.cpp:997
AmpsPdfFlexiFastCPV_mod::_r_dt
RooRealVar * _r_dt
Definition: BsDsKpipi_TD_ampFit.cpp:646
AmpsPdfFlexiFastCPV_mod::_cos_coeff
DecRateCoeff_Bd * _cos_coeff
Definition: BsDsKpipi_TD_ampFit.cpp:662
AmpsPdfFlexiFastCPV_mod::_w
FitParameter & _w
Definition: BsDsKpipi_TD_ampFit.cpp:638
AmpsPdfFlexiFastCPV_time_integrated
Definition: BsDsKpipi_TD_ampFit.cpp:993
main
int main(int argc, char **argv)
Definition: BsDsKpipi_TD_ampFit.cpp:2527
AmpsPdfFlexiFastCPV::getVal_noPs
virtual double getVal_noPs(IDalitzEvent *evt)
Definition: BsDsKpipi_TD_ampFit.cpp:971
basisType
basisType
Definition: TimePdfMaster.h:43
AmpRatio::AmpRatio
AmpRatio(FitParameter &re, FitParameter &im, int f)
Definition: BsDsKpipi_TD_ampFit.cpp:162
MINT::Minimiser::printResultVsInput
void printResultVsInput(std::ostream &os=std::cout) const
Definition: Minimiser.cpp:501
FitAmpListBase::addAsList
virtual int addAsList(const FitAmpListBase &other, double factor=1)
Definition: FitAmpListBase.cpp:69
noBasis
Definition: BsDsKpipi_TD_ampFit.cpp:244
AmpsPdfFlexiFastCPV_mod::_ampsSum
AmpsPdfFlexiFast * _ampsSum
Definition: BsDsKpipi_TD_ampFit.cpp:628
AmpsPdfFlexiFastCPV_time_integrated::_dm
FitParameter & _dm
Definition: BsDsKpipi_TD_ampFit.cpp:1007
TimePdf_mod::un_normalised
double un_normalised(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:1224
MINT::Neg2LLSum::getVal
double getVal()
Definition: Neg2LLSum.cpp:145
AmpsPdfFlexiFastCPV_mod::_r_scale_sigma_dt
RooRealVar * _r_scale_sigma_dt
Definition: BsDsKpipi_TD_ampFit.cpp:652
AmpsPdfFlexiFastCPV_time_integrated::_intA
double _intA
Definition: BsDsKpipi_TD_ampFit.cpp:1011
AmpsPdfFlexiFastCPV_time_integrated::printIntegralVals
void printIntegralVals()
Definition: BsDsKpipi_TD_ampFit.cpp:1031
AmpsPdfFlexiFastCPV_mod::_intAbar
double _intAbar
Definition: BsDsKpipi_TD_ampFit.cpp:641
AmpsPdfFlexiFastCPV::AmpsPdfFlexiFastCPV
AmpsPdfFlexiFastCPV(AmpsPdfFlexiFast *amps1, AmpsPdfFlexiFast *amps2, AmpsPdfFlexiFast *ampsSum, MINT::FitParameter &r, MINT::FitParameter &delta, MINT::FitParameter &gamma, MINT::FitParameter &tau, MINT::FitParameter &dGamma, MINT::FitParameter &dm, MINT::FitParameter &eff_tag, MINT::FitParameter &w)
Definition: BsDsKpipi_TD_ampFit.cpp:985
sinhBasis
Definition: BsDsKpipi_TD_ampFit.cpp:246
TimePdf::_D_bar
FitParameter & _D_bar
Definition: BsDsKpipi_TD_ampFit.cpp:1124
BaseGenerator::dontSaveEvents
void dontSaveEvents()
Definition: BaseGenerator.cpp:80
AmpsPdfFlexiFastCPV_mod::getVal_withPs
virtual double getVal_withPs(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:747
AmpsPdfFlexiFastCPV_mod::printIntegralVals
void printIntegralVals()
Definition: BsDsKpipi_TD_ampFit.cpp:683
AmpsPdfFlexiFastCPV_time_integrated::_dGamma
FitParameter & _dGamma
Definition: BsDsKpipi_TD_ampFit.cpp:1006
AmpsPdfFlexiFastCPV::getVal
virtual double getVal(IDalitzEvent *evt)
Definition: BsDsKpipi_TD_ampFit.cpp:963
TimePdf_mod::getVal
virtual double getVal(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:1248
AmpsPdfFlexiFastCPV_time_integrated::_tau
FitParameter & _tau
Definition: BsDsKpipi_TD_ampFit.cpp:1005
AmpsPdfFlexiFastCPV_mod::endFit
void endFit()
Definition: BsDsKpipi_TD_ampFit.cpp:678
AmpsPdfFlexiFastCPV_time_integrated::getVal_withPs
virtual double getVal_withPs(IDalitzEvent *evt)
Definition: BsDsKpipi_TD_ampFit.cpp:1091
DecRateCoeff_Bd::kCosh
Definition: DecRateCoeff_Bd.h:27
CPV_amp::_im
FitParameter & _im
Definition: BsDsKpipi_TD_ampFit.cpp:173
AmpsPdfFlexiFast::_integratorSource
NamedParameter< std::string > _integratorSource
Definition: BsDsKpipi_TD_ampFit.cpp:88
AmpRatio::ComplexVal
complex< double > ComplexVal()
Definition: BsDsKpipi_TD_ampFit.cpp:165
AmpsPdfFlexiFastCPV_mod
Definition: BsDsKpipi_TD_ampFit.cpp:622
FullAmpsPdfFlexiFastCPV::doFinalStatsAndSaveForAmp12
void doFinalStatsAndSaveForAmp12(MINT::Minimiser *min=0, const std::string &fname="FitAmpResults", const std::string &fnameROOT="fitFractions")
Definition: BsDsKpipi_TD_ampFit.cpp:479
AmpsPdfFlexiFastCPV_mod::getVal_noPs
virtual double getVal_noPs(IDalitzEvent *evt)
Definition: BsDsKpipi_TD_ampFit.cpp:758
sinBasis
Definition: BsDsKpipi_TD_ampFit.cpp:245
MINT::A_is_in_B
bool A_is_in_B(const std::string &a, const std::string &b)
Definition: Utils.cpp:34
RooGaussEfficiencyModel::analyticalIntegral
virtual Double_t analyticalIntegral(Int_t code, const char *rangeName) const
Definition: RooGaussEfficiencyModel.cpp:331
DalitzEventList::fromNtuple
bool fromNtuple(TTree *ntp)
Definition: MyDalitzEventList.h:352
DalitzPdfBaseFlexiFastInteg::ComplexSumForMatchingPatterns
std::complex< double > ComplexSumForMatchingPatterns(bool match)
Definition: DalitzPdfBaseFlexiFastInteg.h:193
FullAmpsPdfFlexiFastCPV::getVal_noPs
virtual double getVal_noPs(IDalitzEvent *evt)
Definition: BsDsKpipi_TD_ampFit.cpp:472
TimePdf::_S
FitParameter & _S
Definition: BsDsKpipi_TD_ampFit.cpp:1125
TimePdf::_C
FitParameter & _C
Definition: BsDsKpipi_TD_ampFit.cpp:1122
FitAmpIncoherentSum::getVal
double getVal(IDalitzEvent &evt)
Definition: FitAmpIncoherentSum.cpp:139
makeIntegratorFile
void makeIntegratorFile()
Definition: BsDsKpipi_TD_ampFit.cpp:2488
MINT::FitParameter::mean
double mean() const
Definition: FitParameter.cpp:312
MINT::counted_ptr
Definition: counted_ptr.h:39
BaseGenerator::setWeighted
void setWeighted(bool w=true)
Definition: BaseGenerator.h:62
MINT::FitParameter::setCurrentFitVal
virtual void setCurrentFitVal(double fv)
Definition: FitParameter.cpp:297
AmpsPdfFlexiFast::AmpsPdfFlexiFast
AmpsPdfFlexiFast(const DalitzEventPattern &pat, IFastAmplitudeIntegrable *amps, MinuitParameterSet *mps, double precision=1.e-4, std::string method="efficient", std::string fname="SignalIntegrationEvents.root", bool genMoreEvents=false)
Definition: BsDsKpipi_TD_ampFit.cpp:100
TimePdf_mod::beginFit
void beginFit()
Definition: BsDsKpipi_TD_ampFit.cpp:1219
AmpsPdfFlexiFastCPV_mod::_sin_coeff
DecRateCoeff_Bd * _sin_coeff
Definition: BsDsKpipi_TD_ampFit.cpp:664
ampFit
int ampFit(int step=0)
Definition: BsDsKpipi_TD_ampFit.cpp:1284
CPV_amp::_re
FitParameter & _re
Definition: BsDsKpipi_TD_ampFit.cpp:172
DecRateCoeff_Bd::kSinh
Definition: DecRateCoeff_Bd.h:26
AmpsPdfFlexiFastCPV_mod::getVal
virtual double getVal(IDalitzEvent *evt)
Definition: BsDsKpipi_TD_ampFit.cpp:750
AmpsPdfFlexiFastCPV_mod::beginFit
void beginFit()
Definition: BsDsKpipi_TD_ampFit.cpp:674
TimePdf_mod::getVal_withPs
virtual double getVal_withPs(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:1262
TimePdf::getVal
virtual double getVal(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:1166
AmpsPdfFlexiFastCPV_time_integrated::_gamma
FitParameter & _gamma
Definition: BsDsKpipi_TD_ampFit.cpp:1003
TimePdf::getVal_withPs
virtual double getVal_withPs(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:1180
AmpsPdfFlexiFastCPV::getVal_withPs
virtual double getVal_withPs(IDalitzEvent *evt)
Definition: BsDsKpipi_TD_ampFit.cpp:967
DalitzEvent::getGeneratorPdfRelativeToPhaseSpace
virtual double getGeneratorPdfRelativeToPhaseSpace() const
Definition: DalitzEvent.cpp:396
FitAmpListBase::normalizeAmps
void normalizeAmps(DalitzEventList &evtList)
Definition: FitAmpListBase.cpp:463
FullAmpsPdfFlexiFastCPV::getValForGeneration
double getValForGeneration(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:317
AmpsPdfFlexiFastCPV_mod::_delta
FitParameter & _delta
Definition: BsDsKpipi_TD_ampFit.cpp:631
MINT::NamedParameter::getVector
const std::vector< T > & getVector() const
Definition: NamedParameter.h:216
FullAmpsPdfFlexiFastCPV::_efficiency
RooGaussEfficiencyModel * _efficiency
Definition: BsDsKpipi_TD_ampFit.cpp:281
AmpsPdfFlexiFastCPV_time_integrated::_intAAbar
complex< double > _intAAbar
Definition: BsDsKpipi_TD_ampFit.cpp:1013
AmpsPdfFlexiFastCPV_mod::_amps2
AmpsPdfFlexiFast * _amps2
Definition: BsDsKpipi_TD_ampFit.cpp:627
AmpsPdfFlexiFastCPV_mod::_tau
FitParameter & _tau
Definition: BsDsKpipi_TD_ampFit.cpp:634
IDalitzEvent::getValueFromVector
virtual double getValueFromVector(unsigned int i) const =0
AmpsPdfFlexiFastCPV_mod::_r_t
RooRealVar * _r_t
Definition: BsDsKpipi_TD_ampFit.cpp:645
AmpsPdfFlexiFastCPV_mod::getVal_noPs
virtual double getVal_noPs(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:748
cosBasis
Definition: BsDsKpipi_TD_ampFit.cpp:245
AmpsPdfFlexiFast::_chosenGen
IEventGenerator< IDalitzEvent > * _chosenGen
Definition: BsDsKpipi_TD_ampFit.cpp:87
GeV
static const double GeV
Definition: CLHEPSystemOfUnits.h:152
DalitzPdfBaseFlexiFastInteg::integralForMatchingPatterns
double integralForMatchingPatterns(bool match, int pattern_sign)
Definition: DalitzPdfBaseFlexiFastInteg.h:188
AmpsPdfFlexiFastCPV_time_integrated::_r
FitParameter & _r
Definition: BsDsKpipi_TD_ampFit.cpp:1001
AmpsPdfFlexiFastCPV_time_integrated::getVal_noPs
virtual double getVal_noPs(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:1085
FullAmpsPdfFlexiFastCPV::_max_TAU
NamedParameter< double > _max_TAU
Definition: BsDsKpipi_TD_ampFit.cpp:283
MINT::EventList::size
virtual unsigned int size() const
Definition: EventList.h:59
MINT::anythingToString
std::string anythingToString(const T &anything)
Definition: Utils.h:62
DalitzPdfBaseFlexiFastInteg::histoSet
virtual DalitzHistoSet histoSet() const
Definition: DalitzPdfBaseFlexiFastInteg.cpp:272
AmpsPdfFlexiFastCPV_mod::_r_scale_mean_dt
RooRealVar * _r_scale_mean_dt
Definition: BsDsKpipi_TD_ampFit.cpp:651
AmpsPdfFlexiFastCPV_time_integrated::endFit
void endFit()
Definition: BsDsKpipi_TD_ampFit.cpp:1026
TimePdf_mod::_eff_tag
FitParameter & _eff_tag
Definition: BsDsKpipi_TD_ampFit.cpp:1213
AmpsPdfFlexiFastCPV::getVal
virtual double getVal(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:942
AmpsPdfFlexiFastCPV_mod::AmpsPdfFlexiFastCPV_mod
AmpsPdfFlexiFastCPV_mod(AmpsPdfFlexiFast *amps1, AmpsPdfFlexiFast *amps2, AmpsPdfFlexiFast *ampsSum, MINT::FitParameter &r, MINT::FitParameter &delta, MINT::FitParameter &gamma, MINT::FitParameter &tau, MINT::FitParameter &dGamma, MINT::FitParameter &dm, MINT::FitParameter &eff_tag, MINT::FitParameter &w)
Definition: BsDsKpipi_TD_ampFit.cpp:772
TimePdf::TimePdf
TimePdf(MINT::FitParameter &C, MINT::FitParameter &D, MINT::FitParameter &D_bar, MINT::FitParameter &S, MINT::FitParameter &S_bar, MINT::FitParameter &k, MINT::FitParameter &tau, MINT::FitParameter &dGamma, MINT::FitParameter &dm, MINT::FitParameter &eff_tag, MINT::FitParameter &w)
Definition: BsDsKpipi_TD_ampFit.cpp:1196
AmpsPdfFlexiFastCPV_mod::_efficiency
RooGaussEfficiencyModel * _efficiency
Definition: BsDsKpipi_TD_ampFit.cpp:658
FitAmpIncoherentSum::print
virtual void print(std::ostream &os=std::cout) const
Definition: FitAmpIncoherentSum.cpp:265
TimePdf::getVal
virtual double getVal(IDalitzEvent *evt)
Definition: BsDsKpipi_TD_ampFit.cpp:1183
AmpRatio::_im
FitParameter & _im
Definition: BsDsKpipi_TD_ampFit.cpp:159
AmpsPdfFlexiFastCPV_time_integrated::_delta
FitParameter & _delta
Definition: BsDsKpipi_TD_ampFit.cpp:1002
MINT::FitParameter::err
double err() const
Definition: FitParameter.cpp:331
AmpsPdfFlexiFast::un_normalised_noPs
double un_normalised_noPs(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:91
TimePdf::_eff_tag
FitParameter & _eff_tag
Definition: BsDsKpipi_TD_ampFit.cpp:1132
DecRateCoeff_Bd
cosh/sinh/cos/sin coefficients in decay rate equations
Definition: DecRateCoeff_Bd.h:21
FullAmpsPdfFlexiFastCPV::histoSet
virtual DalitzHistoSet histoSet()
Definition: BsDsKpipi_TD_ampFit.cpp:477
FullAmpsPdfFlexiFastCPV::getVal_withPs
virtual double getVal_withPs(IDalitzEvent *evt)
Definition: BsDsKpipi_TD_ampFit.cpp:468
TimePdf::_S_bar
FitParameter & _S_bar
Definition: BsDsKpipi_TD_ampFit.cpp:1126
TimePdf::getVal_noPs
virtual double getVal_noPs(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:1181
AmpsPdfFlexiFastCPV_mod::_cosh_coeff
DecRateCoeff_Bd * _cosh_coeff
Definition: BsDsKpipi_TD_ampFit.cpp:663
FullAmpsPdfFlexiFastCPV::getVal
virtual double getVal(IDalitzEvent *evt)
Definition: BsDsKpipi_TD_ampFit.cpp:464
AmpsPdfFlexiFastCPV_mod::_intA
double _intA
Definition: BsDsKpipi_TD_ampFit.cpp:640
DecRateCoeff_Bd::kSin
Definition: DecRateCoeff_Bd.h:24
AmpsPdfFlexiFast::ComplexVal_un_normalised_noPs
std::complex< double > ComplexVal_un_normalised_noPs(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:96
CPV_amp_polar::_r
FitParameter & _r
Definition: BsDsKpipi_TD_ampFit.cpp:186
AmpsPdfFlexiFastCPV::doFinalStatsAndSaveForAmp12
void doFinalStatsAndSaveForAmp12(MINT::Minimiser *min=0, const std::string &fname="FitAmpResults", const std::string &fnameROOT="fitFractions")
Definition: BsDsKpipi_TD_ampFit.cpp:978
AmpsPdfFlexiFastCPV_mod::_spline
RooCubicSplineFun * _spline
Definition: BsDsKpipi_TD_ampFit.cpp:657
TimePdf_mod::endFit
void endFit()
Definition: BsDsKpipi_TD_ampFit.cpp:1221
AmpsPdfFlexiFastCPV_mod::_r_mistag
RooRealVar * _r_mistag
Definition: BsDsKpipi_TD_ampFit.cpp:648
TimePdf::_w
FitParameter & _w
Definition: BsDsKpipi_TD_ampFit.cpp:1133
TimePdf_mod
Definition: BsDsKpipi_TD_ampFit.cpp:1202
calculateCoherence
void calculateCoherence(int step=0)
Definition: BsDsKpipi_TD_ampFit.cpp:2180
TimePdf_mod::_r
FitParameter & _r
Definition: BsDsKpipi_TD_ampFit.cpp:1205
DecRateCoeff_Bd::kCos
Definition: DecRateCoeff_Bd.h:25
DalitzEventList::saveAsNtuple
bool saveAsNtuple(const std::string &fname="DalitzEvents.root") const
Definition: MyDalitzEventList.h:324
FitAmpListBase::CConjugateFinalStateSameFitParameters
virtual bool CConjugateFinalStateSameFitParameters()
Definition: FitAmpListBase.cpp:157
TimePdf::getVal_noPs
virtual double getVal_noPs(IDalitzEvent *evt)
Definition: BsDsKpipi_TD_ampFit.cpp:1191
DalitzEventList::generatePhaseSpaceEvents
int generatePhaseSpaceEvents(int NumEvents, const DalitzEventPattern &pat, TRandom *rnd=0)
Definition: MyDalitzEventList.h:49
CPV_amp_polar::_sign
int _sign
Definition: BsDsKpipi_TD_ampFit.cpp:188
AmpsPdfFlexiFastCPV_time_integrated::getVal
virtual double getVal(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:1066
DalitzPdfBaseFlexiFastInteg::doFinalStatsAndSave
void doFinalStatsAndSave(MINT::Minimiser *min=0, const std::string &fname="FitAmpResults.txt", const std::string &fnameROOT="fitFractions.root")
Definition: DalitzPdfBaseFlexiFastInteg.cpp:318
CPV_amp_polar
Definition: BsDsKpipi_TD_ampFit.cpp:185
AmpsPdfFlexiFastCPV_time_integrated::_amps2
AmpsPdfFlexiFast * _amps2
Definition: BsDsKpipi_TD_ampFit.cpp:998
TimePdf::_dGamma
FitParameter & _dGamma
Definition: BsDsKpipi_TD_ampFit.cpp:1130
TimePdf::parametersChanged
void parametersChanged()
Definition: BsDsKpipi_TD_ampFit.cpp:1136
TimePdf_mod::getVal_noPs
virtual double getVal_noPs(IDalitzEvent *evt)
Definition: BsDsKpipi_TD_ampFit.cpp:1273
AmpsPdfFlexiFastCPV::histoSet
virtual DalitzHistoSet histoSet()
Definition: BsDsKpipi_TD_ampFit.cpp:976
TimePdf_mod::_w
FitParameter & _w
Definition: BsDsKpipi_TD_ampFit.cpp:1214
TimePdf::_tau
FitParameter & _tau
Definition: BsDsKpipi_TD_ampFit.cpp:1129
TimePdf_mod::parametersChanged
void parametersChanged()
Definition: BsDsKpipi_TD_ampFit.cpp:1217
FitAmpSum::RealVal
virtual double RealVal(IDalitzEvent &evt)
Definition: FitAmpSum.h:112
AmpsPdfFlexiFastCPV_time_integrated::un_normalised_noPs
double un_normalised_noPs(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:1038
TimePdf::_k
FitParameter & _k
Definition: BsDsKpipi_TD_ampFit.cpp:1127
AmpsPdfFlexiFastCPV_time_integrated::beginFit
void beginFit()
Definition: BsDsKpipi_TD_ampFit.cpp:1022
AmpsPdfFlexiFastCPV_mod::histoSet
virtual DalitzHistoSet histoSet()
Definition: BsDsKpipi_TD_ampFit.cpp:763
AmpsPdfFlexiFastCPV_time_integrated::parametersChanged
void parametersChanged()
Definition: BsDsKpipi_TD_ampFit.cpp:1016
AmpsPdfFlexiFastCPV_mod::_sinh_coeff
DecRateCoeff_Bd * _sinh_coeff
Definition: BsDsKpipi_TD_ampFit.cpp:665
TimePdf::getVal_withPs
virtual double getVal_withPs(IDalitzEvent *evt)
Definition: BsDsKpipi_TD_ampFit.cpp:1187
AmpsPdfFlexiFastCPV_mod::_min_TAU
NamedParameter< double > _min_TAU
Definition: BsDsKpipi_TD_ampFit.cpp:659
BaseGenerator::FillEventList
void FillEventList(DalitzEventList &evtList, int NEvents)
Definition: BaseGenerator.cpp:116
CPV_amp::ComplexVal
complex< double > ComplexVal()
Definition: BsDsKpipi_TD_ampFit.cpp:179
coshBasis
Definition: BsDsKpipi_TD_ampFit.cpp:246
AmpsPdfFlexiFastCPV_time_integrated::_intAbar
double _intAbar
Definition: BsDsKpipi_TD_ampFit.cpp:1012
AmpsPdfFlexiFastCPV::getVal_withPs
virtual double getVal_withPs(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:960
AmpsPdfFlexiFastCPV_mod::_pdf_omega_ss
RooAbsPdf * _pdf_omega_ss
Definition: BsDsKpipi_TD_ampFit.cpp:655
AmpsPdfFlexiFastCPV_mod::doFinalStatsAndSaveForAmp12
void doFinalStatsAndSaveForAmp12(MINT::Minimiser *min=0, const std::string &fname="FitAmpResults", const std::string &fnameROOT="fitFractions")
Definition: BsDsKpipi_TD_ampFit.cpp:765
FitAmpListBase::multiply
virtual void multiply(double r)
Definition: FitAmpListBase.cpp:404
FullAmpsPdfFlexiFastCPV::FullAmpsPdfFlexiFastCPV
FullAmpsPdfFlexiFastCPV(AmpsPdfFlexiFast *amps1, AmpsPdfFlexiFast *amps2, AmpsPdfFlexiFast *ampsSum, MINT::FitParameter &r, MINT::FitParameter &delta, MINT::FitParameter &gamma, MINT::FitParameter &tau, MINT::FitParameter &dGamma, MINT::FitParameter &dm, MINT::FitParameter &eff_tag, MINT::FitParameter &w)
Definition: BsDsKpipi_TD_ampFit.cpp:486
AmpsPdfFlexiFastCPV_mod::getVal
virtual double getVal(IDalitzEvent &evt)
Definition: BsDsKpipi_TD_ampFit.cpp:720
AmpsPdfFlexiFastCPV_time_integrated::histoSet
virtual DalitzHistoSet histoSet()
Definition: BsDsKpipi_TD_ampFit.cpp:1100
TimePdf_mod::_k
FitParameter & _k
Definition: BsDsKpipi_TD_ampFit.cpp:1208
TimePdf_mod::_delta
FitParameter & _delta
Definition: BsDsKpipi_TD_ampFit.cpp:1206
AmpRatio::_re
FitParameter & _re
Definition: BsDsKpipi_TD_ampFit.cpp:158
AmpsPdfFlexiFastCPV_mod::_r_f
RooCategory * _r_f
Definition: BsDsKpipi_TD_ampFit.cpp:649
AmpsPdfFlexiFastCPV_mod::parametersChanged
void parametersChanged()
Definition: BsDsKpipi_TD_ampFit.cpp:668
DecRateCoeff_Bd.h
Generated by   doxygen 1.8.15