MINT2
Functions
genTimeDependent.cpp File Reference
#include "Mint/NamedParameter.h"
#include "Mint/DiskResidentEventList.h"
#include "Mint/DalitzEvent.h"
#include "TFile.h"
#include "TRandom3.h"
#include <ctime>
#include <iostream>
#include <memory>
#include <Mint/SplineGenerator.h>
#include <Mint/TimeDependentGenerator.h>
#include <Mint/Eff3piSymmetric.h>

Go to the source code of this file.

Functions

int genTimeDependent ()
 
int main ()
 

Function Documentation

◆ genTimeDependent()

int genTimeDependent ( )

Definition at line 23 of file genTimeDependent.cpp.

23  {
24  time_t startTime = time(0);
25 
26  TRandom3 ranLux;
27  NamedParameter<int> RandomSeed("RandomSeed", 0);
28  ranLux.SetSeed((int)RandomSeed);
29  gRandom = &ranLux;
30 
32 
33  NamedParameter<int> Nevents("Nevents", 10000);
34 
35  NamedParameter<int> EventPattern("Event Pattern", 421, -321, 211, 211, -211);
36  DalitzEventPattern pat(EventPattern.getVector());
37  DalitzEventPattern cpPat(pat) ;
38  cpPat[0].antiThis() ;
39 
40  SignalGenerator genD0(pat);
41  SignalGenerator genD0bar(cpPat);
42 
43  NamedParameter<int> saveEvents("SaveEvents", 1);
44  NamedParameter<int> genTimeDependent("genTimeDependent", 0);
45  NamedParameter<double> lifetime("lifetime", 0.4101) ;
46  double width = 1./lifetime ;
47  NamedParameter<double> x("x", 0.0039) ;
48  double deltam = x * width ;
49  NamedParameter<double> y("y", 0.0065) ;
50  double deltagamma = y * width * 2. ;
51  NamedParameter<double> qoverp("qoverp", 1.) ;
52  NamedParameter<double> phi("phi", 0.) ;
53 
54  NamedParameter<double> s13("s13", 0., NamedParameterBase::QUIET) ;
55  NamedParameter<double> s23("s23", 0., NamedParameterBase::QUIET) ;
56 
57  NamedParameter<string> efficiencyFile("efficiencyFile", string("/home/ppe/n/nmchugh/SummerProject/DaVinciDev_v44r10p1/AGammaD0Tohhpi0/scripts/mint/h_efficiency.root")) ;
58  NamedParameter<string> h_efficiencyName( "h_efficiencyName", string("h_efficiency") ) ;
59 
60  NamedParameter<int> addExpEffects("addExpEffects", 0) ;
61  NamedParameter <float> resWidth("resWidth", 0.05) ;
62 
63  NamedParameter<string> outputfname("outputFileName", string("pipipi0_1.root"), (char*)0) ;
64 
65  TH1F* h_efficiency = NULL ;
66  Eff3piSymmetric* sEfficiency = NULL;
67  if((bool)addExpEffects){
68  TFile* eff_infile = TFile::Open( ((string) efficiencyFile).c_str() ) ;
69  eff_infile->GetObject(((string) h_efficiencyName).c_str(), h_efficiency) ;
70  sEfficiency = new Eff3piSymmetric();
71  }
72 
73 
74  cout << " got event pattern: " << pat << endl;
75 
76  unique_ptr<TimeDependentGenerator> timedepgen ;
77  if(genTimeDependent){
78  int startinit(time(0)) ;
79  timedepgen.reset(new TimeDependentGenerator(pat,
80  width, deltam, deltagamma, qoverp, phi, &ranLux,
81  h_efficiency, resWidth, (bool)addExpEffects, sEfficiency)) ;
82  cout << "Initialise TimeDependentGenerator took " << time(0) - startinit << " s" << endl ;
83  }
84 
85  if(!(bool)saveEvents)
86  return 0 ;
87 
88  DiskResidentEventList eventList1(pat, outputfname) ;
89  TNtupleD* evttuple = (TNtupleD*)gDirectory->Get("DalitzEventList") ;
90  TTree* ntuple = new TTree("TimeEventList", "TimeEventList") ;
91  evttuple->AddFriend(ntuple) ;
92  int tag(0) ;
93  double decaytime(0.) ;
94  double smeareddecaytime(0.) ;
95  TBranch* tagbranch = ntuple->Branch("tag", &tag, "tag/I") ;
96  TBranch* decaytimebranch = ntuple->Branch("decaytime", &decaytime, "decaytime/D") ;
97  TBranch* smeareddecaytimebranch = ntuple->Branch("smeareddecaytime", &smeareddecaytime, "smeareddecaytime/D") ;
98 
99  cout << "Generating " << Nevents << " signal events (1)." << endl;
100  int startTimeGen(time(0)) ;
101  for(int i = 0 ; i < Nevents ; ++i){
102 
103  if(i%1000 == 0)
104  cout << "Generating candidate " << i << " (" << (time(0)-startTimeGen)/float(i)
105  << " s per candidate)" << endl ;
106 
108  if(genTimeDependent){
109  if(s13 && s23)
110  evt = timedepgen->generate_event(s13, s23) ;
111  else
112  evt = timedepgen->generate_event() ;
113  }
114  else {
115  // Decide if it's a D0 or D0bar that's being generated.
116  SignalGenerator* gen(&genD0) ;
117  int _tag = 1 ;
118  if(ranLux.Rndm() > 0.5){
119  gen = &genD0bar ;
120  _tag = -1 ;
121  }
122 
123  // Generate the decay time of the candidate.
124  double _decaytime = ranLux.Exp(lifetime) ;
125 
126  evt = gen->newEvent() ;
127  evt->setValueInVector(TimeDependentGenerator::GenTimeEvent::ITAG, _tag) ;
128  evt->setValueInVector(TimeDependentGenerator::GenTimeEvent::IDECAYTIME, _decaytime) ;
129  evt->setValueInVector(TimeDependentGenerator::GenTimeEvent::ISMEAREDDECAYTIME, _decaytime) ;
130  }
131  eventList1.Add(*evt) ;
132  tag = evt->getValueFromVector(TimeDependentGenerator::GenTimeEvent::ITAG) ;
133  decaytime = evt->getValueFromVector(TimeDependentGenerator::GenTimeEvent::IDECAYTIME) ;
134  smeareddecaytime = evt->getValueFromVector(TimeDependentGenerator::GenTimeEvent::ISMEAREDDECAYTIME) ;
135  ntuple->Fill() ;
136  }
137 
138  if(genTimeDependent){
139  cout << "Generator efficiency: " << timedepgen->get_gen_efficiency() << endl ;
140  cout << "Generator scale: " << timedepgen->get_scale() << endl ;
141  }
142 
143  cout << "Save data" << endl ;
144  ntuple->Write() ;
145  eventList1.Close() ;
146 
147  cout << "Save Dalitz plots." << endl ;
148  DalitzHistoSet datH = eventList1.histoSet();
149  datH.save("plotsFromEventList.root");
150 
151  cout << " genTimeDependent done. Took " << (time(0) - startTime)/60.
152  << " min. Returning 0." << endl;
153 
154  // For debugging, plot PDF vs time at given point in phase space for each tag.
155  if(s13 && s23){
156  TFile fplots("plotsVsTime.root", "recreate") ;
157  unsigned nbins = 1000 ;
158  float tmin = 0. ;
159  float tmax = 8.2 ;
160  DalitzEvent evt(cpPat, s13, s23) ;
161  for(int tag = -1 ; tag < 2 ; tag += 2){
162  evt.setValueInVector(TimeDependentGenerator::GenTimeEvent::ITAG, tag) ;
163  std::ostringstream tagstr ;
164  tagstr << tag ;
165  TH1F pdf = timedepgen->draw_pdf_vs_time(evt, nbins, tmin, tmax,
166  "pdf_vs_time_tag" + tagstr.str()) ;
167  TH1F env = timedepgen->draw_envelope_vs_time(evt, nbins, tmin, tmax,
168  "env_vs_time_tag" + tagstr.str()) ;
169  pdf.Write() ;
170  env.Write() ;
171  evt = DalitzEvent(pat, s13, s23) ;
172  }
173  fplots.Close() ;
174  }
175 
176  return 0;
177 }
static void AutogenerateFitFile(const std::string &fname="protoFitAmplitudeFile.txt", const DalitzEventPattern &pat=DalitzEventPattern::NoPattern)
int genTimeDependent()
bool save(const std::string &filename="DalitzHistos.root") const

◆ main()

int main ( )

Definition at line 180 of file genTimeDependent.cpp.

180  {
181 
182  return genTimeDependent();
183 
184 }
int genTimeDependent()