StRoot  1
genRootDecayChain.cc
1 // Program to create ROOT files for EvtGen validation plots,
2 // generalising storage of daughter information for each decay chain step.
3 // There are three integers to keep track of the decay tree history:
4 // dVtx = present decay vertex integer number
5 // pVtx = parent decay vertex integer number
6 // daug = daughter order number for the given particle
7 //
8 // Consider the following decay chain history, where Vj is vertex number j:
9 // V1
10 // A -> B + C <---- First chain level
11 // | |
12 // +--> B1 + B2 +--> C1 C2 C3 <---- Second chain level
13 // V2 | V5 |
14 // +--> B1a B1b +--> C1a C1b C1c <---- Third chain level
15 // V3 | V6
16 // +--> d1 d2 <---- Fourth chain level
17 // V4
18 //
19 // Since the decay tree is stored recursively, the information order follows
20 // the immediate, sequential, sub-decay vertex order V1, V2, V3, V4, V5, V6:
21 // A, B, B1, B1a, B1b, d1, d2, B2, C, C1, C1a, C1b, C1c, C2, C3
22 //
23 // The unique set of integer "attributes" for each particle is as follows:
24 //
25 // Particle dVtx pVtx daug
26 // A 0 0 0 The mother particle
27 // B 1 0 1
28 // B1 2 1 1
29 // B1a 3 2 1
30 // B1b 3 2 2
31 // d1 4 3 1
32 // d2 4 3 2
33 // B2 2 1 2
34 // C 1 0 2
35 // C1 5 1 1
36 // C1a 6 5 1
37 // C1b 6 5 2
38 // C1c 6 5 3
39 // C2 5 1 2
40 // C3 5 1 3
41 //
42 // The decay tree will be defined by the decay file. Additional photons from PHOTOS will
43 // only appear as appended, extra daughters on the RHS of a given vertex, e.g. if a FSR
44 // photon is added to the RHS of V4, there will be a "d3" particle with attributes (4,3,3)
45
46 #include "genRootDecayChain.hh"
47
48 #include "EvtGen/EvtGen.hh"
49 #include "EvtGenBase/EvtRandomEngine.hh"
51 #include "EvtGenBase/EvtDecayBase.hh"
52 #include "EvtGenBase/EvtParticle.hh"
53 #include "EvtGenBase/EvtParticleFactory.hh"
54 #include "EvtGenBase/EvtPatches.hh"
55 #include "EvtGenBase/EvtPDL.hh"
56 #include "EvtGenBase/EvtSimpleRandomEngine.hh"
57 #include "EvtGenBase/EvtMTRandomEngine.hh"
58
59 #ifdef EVTGEN_EXTERNAL
60 #include "EvtGenExternal/EvtExternalGenList.hh"
61 #endif
62
63 #include "TStyle.h"
64 #include "TROOT.h"
65
66 #include <iostream>
67 #include <list>
68 #include <string>
69
70 using std::cout;
71 using std::endl;
72 using std::string;
73
74 genRootDecayChain::genRootDecayChain(const string& decayFileName,
75  const string& rootFileName,
76  const string& parentName,
77  int nEvents, bool storeMtmXYZ) :
78  _decayFileName(decayFileName),
79  _rootFileName(rootFileName),
80  _parentName(parentName),
81  _nEvents(nEvents),
82  _storeMtmXYZ(storeMtmXYZ),
83  _theFile(0),
84  _theTree(0),
85  _theCanvas(0)
86 {
87
88  _theFile = new TFile(rootFileName.c_str(), "recreate");
89  _theTree = new TTree("Data", "Data");
90  _theTree->SetDirectory(_theFile);
91
92  gROOT->SetStyle("Plain");
93  gStyle->SetOptStat(0);
94
95  _theCanvas = new TCanvas("theCanvas", "", 900, 700);
96  _theCanvas->Clear();
97  _theCanvas->UseCurrentStyle();
98
99 }
100
101 genRootDecayChain::~genRootDecayChain() {
102
103  _theFile->Close();
104  delete _theCanvas;
105
106 }
107
108 void genRootDecayChain::run() {
109
110  this->initTree();
111  this->generateEvents();
112  this->writeTree();
113 }
114
115 void genRootDecayChain::initTree() {
116
117  // Initialise internal variables
118  _eventId = 0;
119  _PDGId = 0;
120  _dVtx = 0, _pVtx = 0, _daug = 0;
121  _px = 0.0, _py = 0.0, _pz = 0.0;
122  _p = 0.0, _E = 0.0, _m = 0.0, _t = 0.0;
123
124  // Set-up the TTree
125  _theTree->Branch("eventId", &_eventId, "eventId/I");
126  _theTree->Branch("PDGId", &_PDGId, "PDGId/I");
127  _theTree->Branch("dVtx", &_dVtx, "dVtx/I");
128  _theTree->Branch("pVtx", &_pVtx, "pVtx/I");
129  _theTree->Branch("daug", &_daug, "daug/I");
130  _theTree->Branch("p", &_p, "p/D");
131  _theTree->Branch("E", &_E, "E/D");
132  _theTree->Branch("pL", &_pL, "pL/D");
133  _theTree->Branch("EL", &_EL, "EL/D");
134  _theTree->Branch("m", &_m, "m/D");
135
136  if (_storeMtmXYZ) {
137
138  // Store momentum components as well
139  _theTree->Branch("px", &_px, "px/D");
140  _theTree->Branch("py", &_py, "py/D");
141  _theTree->Branch("pz", &_pz, "pz/D");
142  _theTree->Branch("pxL", &_pxL, "pxL/D");
143  _theTree->Branch("pyL", &_pyL, "pyL/D");
144  _theTree->Branch("pzL", &_pzL, "pzL/D");
145
146  }
147
149  _theTree->Branch("t", &_t, "t/D");
150
151 }
152
153 void genRootDecayChain::writeTree() {
154
155  _theFile->cd();
156  _theTree->Write();
157
158 }
159
160 void genRootDecayChain::generateEvents() {
161
162  EvtRandomEngine* randomEngine = 0;
164  std::list<EvtDecayBase*> extraModels;
165
166  // Define the random number generator
167 #ifdef EVTGEN_CPP11
168  // Use the Mersenne-Twister generator (C++11 only)
169  randomEngine = new EvtMTRandomEngine();
170 #else
171  randomEngine = new EvtSimpleRandomEngine();
172 #endif
173
174  // Initialize the generator - read in the decay table and particle properties.
175  // For our validation purposes, we just want to read in one decay file
176
177 #ifdef EVTGEN_EXTERNAL
178  bool useEvtGenRandom(false);
179  EvtExternalGenList genList(true, "", "gamma", useEvtGenRandom);
181  extraModels = genList.getListOfModels();
182 #endif
183
184  int mixingType(1);
185  bool useXml(false);
186
187  EvtGen evtGen(_decayFileName.c_str(), "../evt.pdl", randomEngine,
189
190  EvtParticle* theParent(0);
191
192  EvtId theId = EvtPDL::getId(_parentName);
193  if (theId.getId() == -1 && theId.getAlias() == -1) {
194  cout<<"Error. Could not find valid EvtId for "<<_parentName<<endl;
195  return;
196  }
197
198  // Loop to create nEvents
199  int i;
200  for (i = 0; i < _nEvents; i++) {
201
202  // Parent particle 4-momentum
203  EvtVector4R pInit(EvtPDL::getMass(theId), 0.0, 0.0, 0.0);
204
205  _eventId = i;
206
207  if (i%1000 == 0) {cout<<"Event number = "<<i+1<<" out of "<<_nEvents<<std::endl;}
208
209  theParent = EvtParticleFactory::particleFactory(theId, pInit);
210  if (theParent->getSpinStates() == 3) {theParent->setVectorSpinDensity();}
211
212  // Generate the event
213  evtGen.generateDecay(theParent);
214
215  // Decay vertex index
216  theParent->setAttribute("dVtx", 0);
217  // Parent vertex index
218  theParent->setAttribute("pVtx", 0);
219  // Daughter index
220  theParent->setAttribute("daug", 0);
221
222  int nDaug = theParent->getNDaug();
223  int iDaug(0);
224
225  storeTreeInfo(theParent);
226
227  //cout<<"Event number = "<<i<<endl;
228  //theParent->printTree();
229
230  // Initialise the vertex number
231  _vertexNo = 1;
232
233  // Loop over the daughter tracks
234  for (iDaug = 0; iDaug < nDaug; iDaug++) {
235
236  EvtParticle* daug = theParent->getDaug(iDaug);
237  // Decay vertex index
238  daug->setAttribute("dVtx", 1);
239  // Parent decay vertex
240  daug->setAttribute("pVtx", 0);
241  // Daughter index: 1,2,..,nDaug
242  daug->setAttribute("daug", iDaug+1);
243
244  // Recursively store the daughter information
245  this->storeDaughterInfo(daug);
246
247  } // daughter loop
248
249  theParent->deleteTree();
250
251  } // event loop
252
253  delete randomEngine;
254
255 }
256
257 void genRootDecayChain::storeDaughterInfo(EvtParticle* theParticle) {
258
259  if (!theParticle) {return;}
260
261  // Store the particle information in the TTree
262  storeTreeInfo(theParticle);
263
264  // Loop over the particle's own daughter tracks and call this function
265  // recursively, keeping track of the decay chain order number
266  int nDaug = theParticle->getNDaug();
267  int iDaug(0);
268
269  // Increment the decay vertex integer if this particle has daughters
270  if (nDaug > 0) {_vertexNo++;}
271
272  // The parent vertex index for the daughters is equal to the
273  // current particle decay vertex index
274  int parentVtx = theParticle->getAttribute("dVtx");
275
276  // First, we need to loop over the given particle's daughters and set the
277  // attributes so we keep track of the decay tree history at this chain level
278  for (iDaug = 0; iDaug < nDaug; iDaug++) {
279
280  EvtParticle* daug = theParticle->getDaug(iDaug);
281
282  // Set the attributes for the daughters
283  daug->setAttribute("dVtx", _vertexNo);
284  daug->setAttribute("pVtx", parentVtx);
285  daug->setAttribute("daug", iDaug+1);
286
287  }
288
289  // Now loop over the daughters again and recursively call this function to
290  // follow the rest of the decay tree history
291  for (iDaug = 0; iDaug < nDaug; iDaug++) {
292
293  EvtParticle* daug = theParticle->getDaug(iDaug);
294  storeDaughterInfo(daug);
295
296  }
297
298 }
299
300 void genRootDecayChain::storeTreeInfo(EvtParticle* theParticle) {
301
302  if (!theParticle) {return;}
303
304  // Store the particle information in the TTree by first setting the internal
305  // variables and then calling Fill()
306  _dVtx = theParticle->getAttribute("dVtx");
307  _pVtx = theParticle->getAttribute("pVtx");
308  _daug = theParticle->getAttribute("daug");
309
310  //cout<<"Particle "<<theParticle->getName()<<", dVtx = "<<_dVtx
311  //<<", pVtx = "<<_pVtx<<", daug = "<<_daug<<endl;
312
313  EvtVector4R p4Lab = theParticle->getP4Lab(); // lab frame = mother frame
314  EvtVector4R p4 = theParticle->getP4(); // frame = parents frame
315
316  // PDG id
317  _PDGId = EvtPDL::getStdHep(theParticle->getId());
318
319  // 4-momenta in parent frame
320  _E = p4.get(0);
321  _px = p4.get(1);
322  _py = p4.get(2);
323  _pz = p4.get(3);
324  _p = sqrt(_px*_px + _py*_py + _pz*_pz);
325
326  // 4-momenta in lab frame
327  _EL = p4Lab.get(0);
328  _pxL = p4Lab.get(1);
329  _pyL = p4Lab.get(2);
330  _pzL = p4Lab.get(3);
331  _pL = sqrt(_pxL*_pxL + _pyL*_pyL + _pzL*_pzL);
332
333  // Rest mass and lifetime
334  _m = theParticle->mass();
336
337  // Store the information in the TTree
338  _theTree->Fill();
339
340 }
341
342 int main(int argc, char** argv) {
343
344  // Use the B0 -> K'*0 gamma, K'*0 -> K+ pi- decays as an example
345  string decayFileName("BKstarGamma.dec");
346  //string decayFileName("BuDst0rhop.dec");
347  if (argc > 1) {decayFileName = argv[1];}
348  cout<<"Decay file name is "<<decayFileName<<endl;
349
350  string rootFileName("BKstarGamma.root");
351  //string rootFileName("BuDst0rhop.root");
352  if (argc > 2) {rootFileName = argv[2];}
353  cout<<"Root file name is "<<rootFileName<<endl;
354
355  string parentName("B0");
356  //string parentName("B-");
357  if (argc > 3) {parentName = argv[3];}
358  cout<<"Parent name is "<<parentName<<endl;
359
360  int nEvents(100000);
361  if (argc > 4) {nEvents = atoi(argv[4]);}
362
363  bool storeMtmXYZ = true;
364
365  genRootDecayChain myGen(decayFileName, rootFileName, parentName, nEvents, storeMtmXYZ);
366  myGen.run();
367
368  return 0;
369
370 }
double mass() const
const EvtVector4R & getP4() const
EvtParticle * getDaug(int i)
EvtVector4R getP4Lab() const
size_t getNDaug() const
Definition: EvtId.hh:27
Definition: EvtGen.hh:46
EvtId getId() const