MINT2
Classes | Enumerations | Functions
BsDsKpipi_TD_ampFit.cpp File Reference
#include "Mint/FitParameter.h"
#include "Mint/NamedParameter.h"
#include "Mint/Minimiser.h"
#include "Mint/Neg2LL.h"
#include "Mint/Neg2LLSum.h"
#include "Mint/DalitzEventList.h"
#include "Mint/NamedDecayTreeList.h"
#include "Mint/DecayTree.h"
#include "Mint/DiskResidentEventList.h"
#include "Mint/CLHEPPhysicalConstants.h"
#include "Mint/CLHEPSystemOfUnits.h"
#include "Mint/PdfBase.h"
#include "Mint/DalitzPdfBase.h"
#include "Mint/DalitzPdfBaseFastInteg.h"
#include "Mint/DalitzPdfBaseFlexiFastInteg.h"
#include "Mint/FitAmplitude.h"
#include "Mint/FitAmpSum.h"
#include "Mint/FitAmpIncoherentSum.h"
#include "Mint/DalitzEvent.h"
#include "Mint/AmpRatios.h"
#include "Mint/IEventGenerator.h"
#include "Mint/DalitzBWBoxSet.h"
#include "Mint/DalitzBoxSet.h"
#include "Mint/SignalGenerator.h"
#include "Mint/FromFileGenerator.h"
#include "Mint/DalitzSumPdf.h"
#include "Mint/cexp.h"
#include "Mint/DalitzPdfNormChecker.h"
#include "Mint/IFastAmplitudeIntegrable.h"
#include "Mint/DalitzPdfSaveInteg.h"
#include "Mint/Chi2Binning.h"
#include "Mint/FitAmpList.h"
#include "Mint/DalitzPdfBaseMCInteg.h"
#include "TGraph.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TNtupleD.h"
#include "TTree.h"
#include "TLegend.h"
#include <TStyle.h>
#include "TRandom2.h"
#include "TRandom3.h"
#include <ctime>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <TROOT.h>
#include "RooRealConstant.h"
#include "RooAbsReal.h"
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGlobalFunc.h"
#include "RooDecay.h"
#include "RooBDecay.h"
#include "RooPlot.h"
#include "RooEffProd.h"
#include "RooGenericPdf.h"
#include "RooGaussModel.h"
#include "RooProdPdf.h"
#include "RooConstVar.h"
#include "RooCategory.h"
#include "Mint/RooCubicSplineFun.h"
#include "Mint/RooCubicSplineKnot.h"
#include "Mint/RooGaussEfficiencyModel.h"
#include "Mint/DecRateCoeff_Bd.h"

Go to the source code of this file.

Classes

class  AmpsPdfFlexiFast
 
class  AmpRatio
 
class  CPV_amp
 
class  CPV_amp_polar
 
class  FullAmpsPdfFlexiFastCPV
 
class  AmpsPdfFlexiFastCPV_mod
 
class  AmpsPdfFlexiFastCPV
 
class  AmpsPdfFlexiFastCPV_time_integrated
 
class  TimePdf
 
class  TimePdf_mod
 

Enumerations

enum  basisType {
  noBasis =0, expBasis = 3, sinBasis =13, cosBasis =23,
  sinhBasis =63, coshBasis =53, noBasis =0, expBasis = 3,
  sinBasis =13, cosBasis =23, sinhBasis =63, coshBasis =53
}
 

Functions

void AddScaledAmpsToList (FitAmpSum &fas_tmp, FitAmpSum &fas, FitAmpSum &fasCC, std::string name, counted_ptr< IReturnComplex > &r_plus, counted_ptr< IReturnComplex > &r_minus)
 
std::vector< double > coherenceFactor (FitAmpSum &fas, FitAmpSum &fas_bar, double r, double delta, DalitzEventList &eventList)
 
int ampFit (int step=0)
 
void calculateCoherence (int step=0)
 
void coherence_plot ()
 
void makeIntegratorFile ()
 
int main (int argc, char **argv)
 

Enumeration Type Documentation

◆ basisType

enum basisType
Enumerator
noBasis 
expBasis 
sinBasis 
cosBasis 
sinhBasis 
coshBasis 
noBasis 
expBasis 
sinBasis 
cosBasis 
sinhBasis 
coshBasis 

Definition at line 244 of file BsDsKpipi_TD_ampFit.cpp.

244  { noBasis=0 , expBasis= 3
245  , sinBasis=13, cosBasis=23
246  , sinhBasis=63, coshBasis=53 };
expBasis
Definition: BsDsKpipi_TD_ampFit.cpp:244
noBasis
Definition: BsDsKpipi_TD_ampFit.cpp:244
sinhBasis
Definition: BsDsKpipi_TD_ampFit.cpp:246
sinBasis
Definition: BsDsKpipi_TD_ampFit.cpp:245
cosBasis
Definition: BsDsKpipi_TD_ampFit.cpp:245
coshBasis
Definition: BsDsKpipi_TD_ampFit.cpp:246

Function Documentation

◆ AddScaledAmpsToList()

void AddScaledAmpsToList ( FitAmpSum &  fas_tmp,
FitAmpSum &  fas,
FitAmpSum &  fasCC,
std::string  name,
counted_ptr< IReturnComplex > &  r_plus,
counted_ptr< IReturnComplex > &  r_minus 
)

Definition at line 199 of file BsDsKpipi_TD_ampFit.cpp.

199  {
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 }
FitAmpSum::GetCloneOfSubsetSameFitParameters
virtual MINT::counted_ptr< FitAmpListBase > GetCloneOfSubsetSameFitParameters(std::string name) const
Definition: FitAmpSum.cpp:147
FitAmpListBase::addAsList
virtual int addAsList(const FitAmpListBase &other, double factor=1)
Definition: FitAmpListBase.cpp:69
MINT::counted_ptr
Definition: counted_ptr.h:39

◆ ampFit()

int ampFit ( int  step = 0)

Definition at line 1284 of file BsDsKpipi_TD_ampFit.cpp.

1284  {
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 }
MINT::EventList::Add
virtual bool Add(const EVENT_TYPE &evt)
Definition: EventList.h:63
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
FitAmplitude::AutogenerateFitFile
static void AutogenerateFitFile(const std::string &fname="protoFitAmplitudeFile.txt", const DalitzEventPattern &pat=DalitzEventPattern::NoPattern)
Definition: FitAmplitude.cpp:50
pi
static const double pi
Definition: CLHEPPhysicalConstants.h:44
s
static const double s
Definition: CLHEPSystemOfUnits.h:126
DalitzEvent::setValueInVector
virtual void setValueInVector(unsigned int i, double value)
Definition: DalitzEvent.cpp:415
coherenceFactor
std::vector< double > coherenceFactor(FitAmpSum &fas, FitAmpSum &fas_bar, double r, double delta, DalitzEventList &eventList)
Definition: BsDsKpipi_TD_ampFit.cpp:212
AmpsPdfFlexiFastCPV_time_integrated
Definition: BsDsKpipi_TD_ampFit.cpp:993
MINT::A_is_in_B
bool A_is_in_B(const std::string &a, const std::string &b)
Definition: Utils.cpp:34
DalitzEventList::fromNtuple
bool fromNtuple(TTree *ntp)
Definition: MyDalitzEventList.h:352
MINT::counted_ptr
Definition: counted_ptr.h:39
DalitzEvent::getGeneratorPdfRelativeToPhaseSpace
virtual double getGeneratorPdfRelativeToPhaseSpace() const
Definition: DalitzEvent.cpp:396
GeV
static const double GeV
Definition: CLHEPSystemOfUnits.h:152
MINT::EventList::size
virtual unsigned int size() const
Definition: EventList.h:59
MINT::anythingToString
std::string anythingToString(const T &anything)
Definition: Utils.h:62
TimePdf_mod
Definition: BsDsKpipi_TD_ampFit.cpp:1202
DalitzEventList::saveAsNtuple
bool saveAsNtuple(const std::string &fname="DalitzEvents.root") const
Definition: MyDalitzEventList.h:324
DalitzEventList::generatePhaseSpaceEvents
int generatePhaseSpaceEvents(int NumEvents, const DalitzEventPattern &pat, TRandom *rnd=0)
Definition: MyDalitzEventList.h:49
CPV_amp_polar
Definition: BsDsKpipi_TD_ampFit.cpp:185

◆ calculateCoherence()

void calculateCoherence ( int  step = 0)

Definition at line 2180 of file BsDsKpipi_TD_ampFit.cpp.

2180  {
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 }
MINT::EventList::Add
virtual bool Add(const EVENT_TYPE &evt)
Definition: EventList.h:63
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
FitAmplitude::AutogenerateFitFile
static void AutogenerateFitFile(const std::string &fname="protoFitAmplitudeFile.txt", const DalitzEventPattern &pat=DalitzEventPattern::NoPattern)
Definition: FitAmplitude.cpp:50
s
static const double s
Definition: CLHEPSystemOfUnits.h:126
coherenceFactor
std::vector< double > coherenceFactor(FitAmpSum &fas, FitAmpSum &fas_bar, double r, double delta, DalitzEventList &eventList)
Definition: BsDsKpipi_TD_ampFit.cpp:212
DalitzEventList::fromNtuple
bool fromNtuple(TTree *ntp)
Definition: MyDalitzEventList.h:352
MINT::counted_ptr
Definition: counted_ptr.h:39
GeV
static const double GeV
Definition: CLHEPSystemOfUnits.h:152
MINT::EventList::size
virtual unsigned int size() const
Definition: EventList.h:59
MINT::anythingToString
std::string anythingToString(const T &anything)
Definition: Utils.h:62
DalitzEventList::generatePhaseSpaceEvents
int generatePhaseSpaceEvents(int NumEvents, const DalitzEventPattern &pat, TRandom *rnd=0)
Definition: MyDalitzEventList.h:49
CPV_amp_polar
Definition: BsDsKpipi_TD_ampFit.cpp:185

◆ coherence_plot()

void coherence_plot ( )

Definition at line 2403 of file BsDsKpipi_TD_ampFit.cpp.

2403  {
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 }
MINT::EventList::Add
virtual bool Add(const EVENT_TYPE &evt)
Definition: EventList.h:63
MINT::PolymorphVector::push_back
void push_back(const T &c)
Definition: PolymorphVector.h:50
s
static const double s
Definition: CLHEPSystemOfUnits.h:126
DalitzEventList::fromNtuple
bool fromNtuple(TTree *ntp)
Definition: MyDalitzEventList.h:352
GeV
static const double GeV
Definition: CLHEPSystemOfUnits.h:152
MINT::EventList::size
virtual unsigned int size() const
Definition: EventList.h:59

◆ coherenceFactor()

std::vector<double> coherenceFactor ( FitAmpSum &  fas,
FitAmpSum &  fas_bar,
double  r,
double  delta,
DalitzEventList &  eventList 
)

Definition at line 212 of file BsDsKpipi_TD_ampFit.cpp.

212  {
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 }
MINT::PolymorphVector::push_back
void push_back(const T &c)
Definition: PolymorphVector.h:50
pi
static const double pi
Definition: CLHEPPhysicalConstants.h:44
FitAmpSum::getVal
virtual std::complex< double > getVal(IDalitzEvent &evt)
Definition: FitAmpSum.cpp:182
MINT::EventList::size
virtual unsigned int size() const
Definition: EventList.h:59

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 2527 of file BsDsKpipi_TD_ampFit.cpp.

2527  {
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 }
makeIntegratorFile
void makeIntegratorFile()
Definition: BsDsKpipi_TD_ampFit.cpp:2488
ampFit
int ampFit(int step=0)
Definition: BsDsKpipi_TD_ampFit.cpp:1284

◆ makeIntegratorFile()

void makeIntegratorFile ( )

Definition at line 2488 of file BsDsKpipi_TD_ampFit.cpp.

2488  {
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 }
MINT::EventList::Add
virtual bool Add(const EVENT_TYPE &evt)
Definition: EventList.h:63
FitAmplitude::AutogenerateFitFile
static void AutogenerateFitFile(const std::string &fname="protoFitAmplitudeFile.txt", const DalitzEventPattern &pat=DalitzEventPattern::NoPattern)
Definition: FitAmplitude.cpp:50
s
static const double s
Definition: CLHEPSystemOfUnits.h:126
GeV
static const double GeV
Definition: CLHEPSystemOfUnits.h:152
MINT::EventList::size
virtual unsigned int size() const
Definition: EventList.h:59
DalitzEventList::saveAsNtuple
bool saveAsNtuple(const std::string &fname="DalitzEvents.root") const
Definition: MyDalitzEventList.h:324
DalitzEventList::generatePhaseSpaceEvents
int generatePhaseSpaceEvents(int NumEvents, const DalitzEventPattern &pat, TRandom *rnd=0)
Definition: MyDalitzEventList.h:49
Generated by   doxygen 1.8.15