1 // module identifier line...
2 // Author: Brian Thorsbro, 24/6-2014
11 #include "THepMCParser.h"
12 #include "HepMC/IO_GenEvent.h"
15 #include "TClonesArray.h"
17 #include "TParticle.h"
18 #include "TDatabasePDG.h"
23 THepMCParser::THepMCParser(const char * infile) : fTree(0)
25 HepMC::IO_BaseClass * events = new HepMC::IO_GenEvent(infile, std::ios::in);
28 events = 0; // nullptr
30 THepMCParser::THepMCParser(HepMC::IO_BaseClass * events) : fTree(0)
34 void THepMCParser::init(HepMC::IO_BaseClass * events)
36 int particlecount = 0;
37 fTree = new TTree("treeEPOS","Tree EPOS");
38 TClonesArray * array = new TClonesArray("TParticle");
39 // array->BypassStreamer();
40 fTree->Branch("Particles",&array); // more flags?
41 THepMCParser::HeavyIonHeader_t heavyIonHeader;
42 fTree->Branch("HeavyIonInfo", &heavyIonHeader, THepMCParser::fgHeavyIonHeaderBranchString);
43 THepMCParser::PdfHeader_t pdfHeader;
44 fTree->Branch("PdfInfo", &pdfHeader, THepMCParser::fgPdfHeaderBranchString);
45 HepMC::GenEvent* evt = 0; // nullptr
46 while ((evt = events->read_next_event())) {
47 string errMsg1 = ParseGenEvent2TCloneArray(evt,array);
48 string errMsg2 = ParseGenEvent2HeaderStructs(evt,heavyIonHeader,pdfHeader);
49 if (errMsg1.length() == 0 && errMsg2.length() == 0) {
52 if (errMsg1.length() != 0) cerr << errMsg1 << endl;
53 if (errMsg2.length() != 0) cerr << errMsg2 << endl;
55 particlecount += array->Capacity();
59 // array = 0; // nullptr
60 cout << " parsed " << particlecount << " particles" << endl;
64 TTree * THepMCParser::GetTTree()
68 void THepMCParser::WriteTTreeToFile(const char *outfile)
70 TFile * f = new TFile(outfile, "recreate");
78 // Based on a validator written by Peter Hristov, CERN
79 bool THepMCParser::IsValidMotherDaughtersConsitency(bool useStdErr, bool requireSecondMotherBeforeDaughters)
82 TClonesArray * array = new TClonesArray("TParticle");
83 TBranch* branch = fTree->GetBranch("Particles");
84 branch->SetAddress(&array);
85 Int_t count = branch->GetEntries();
86 for (Int_t idx=0; idx<count; ++idx) {
88 branch->GetEntry(idx); // "fill" the array
89 Int_t nkeep = array->GetEntriesFast();
90 for (Int_t i=0; i<nkeep; i++) {
91 TParticle * part = (TParticle*)array->AddrAt(i);
92 Int_t mum1 = part->GetFirstMother();
93 Int_t mum2 = part->GetSecondMother();
94 Int_t fd = part->GetFirstDaughter();
95 Int_t ld = part->GetLastDaughter();
96 if (mum1>-1 && i<mum1) {
98 if (useStdErr) cerr << "Problem: first_mother(" << mum1 << ") after daughter(" << i << ")" << endl;
100 if (mum2>-1 && i<mum2 && requireSecondMotherBeforeDaughters) {
102 if (useStdErr) cerr << "Problem: second_mother(" << mum2 << ") after daughter(" << i << ")" << endl;
106 if (useStdErr) cerr << "Problem: first_daughter(" << fd << ") > last_daughter(" << ld << ")" << endl;
108 for (Int_t id=TMath::Max(fd,0); id<=ld; id++) {
109 TParticle * daughter = (TParticle*)array->AddrAt(id);
110 if (daughter->GetFirstMother() != i && daughter->GetSecondMother() != i) {
112 if (useStdErr) cerr << "Problem: mother("<< i << ") not registered as mother for either first_daughter("
113 << daughter->GetFirstMother() << ") or second_daughter("
114 << daughter->GetSecondMother() << ")" << endl;
124 bool THepMCParser::IsValidParticleInvariantMass(bool useStdErr, bool includeStatusCode2Particles)
127 TClonesArray *array = new TClonesArray("TParticle");
128 TBranch* branch = fTree->GetBranch("Particles");
129 branch->SetAddress(&array);
130 Int_t count = branch->GetEntries();
131 for (Int_t idx=0; idx<count; ++idx) {
133 branch->GetEntry(idx);
134 Int_t nkeep = array->GetEntries();
135 for (Int_t i=0; i<nkeep; i++) {
136 TParticle * parton = (TParticle*)array->AddrAt(i);
137 if (parton->GetStatusCode()==2 && !includeStatusCode2Particles) {
142 Double_t m1 = v.M(); // invariant mass from the particle in the HepMC event
143 TParticlePDG *dbParton = parton->GetPDG();
145 if (useStdErr) cerr << "Warning: could not look up particle with PDG Code: " << parton->GetPdgCode() << endl << endl;
148 Double_t m2 = dbParton->Mass();
151 checkok = abs(m1) < 0.0001; // no such thing as negative mass...
153 checkok = abs(1 - m1/m2) < 0.01;
155 if (!checkok && useStdErr) {
156 cerr << "Problem: " << GetParticleName(parton) << " HepMC:" << m1 << endl;
157 cerr << ListReactionChain(array,i);
169 bool THepMCParser::IsValidVertexInvariantMass(bool useStdErr, bool includeStatusCode2Particles)
172 TClonesArray * array = new TClonesArray("TParticle");
173 TBranch* branch = fTree->GetBranch("Particles");
174 branch->SetAddress(&array);
175 Int_t count = branch->GetEntries();
176 for (Int_t idx=0; idx<count; ++idx) {
178 branch->GetEntry(idx); // "fill" the array
179 TLorentzVector v_st1;
180 TLorentzVector v_st4;
181 Int_t nkeep = array->GetEntriesFast();
182 for (Int_t i=0; i<nkeep; i++) {
183 TParticle * parton = (TParticle*)array->AddrAt(i);
185 parton->Momentum(v_in);
186 if (parton->GetStatusCode()==4) {
188 } else if (parton->GetStatusCode()==1) {
191 if (!includeStatusCode2Particles) { // only check beam particles vs final particles
194 Int_t fd = parton->GetFirstDaughter();
195 Int_t ld = parton->GetLastDaughter();
196 if (fd == -1) continue; // no daughters, continue loop
198 TLorentzVector v_out;
199 bool oneok = false; bool allok = true; bool agreemother2 = true;
200 ostringstream daughterpdg;
201 ostringstream motherpdg;
202 for (Int_t id=TMath::Max(fd,0); id<=ld; id++) {
203 TParticle * daughter = (TParticle*)array->AddrAt(id);
205 daughter->Momentum(v_out);
206 mother2 = daughter->GetSecondMother();
209 daughter->Momentum(d);
211 if (daughter->GetSecondMother() != mother2) agreemother2 = false;
213 if (daughter->GetFirstMother() == i) {
218 daughterpdg << " " << daughter->GetPdgCode();
220 motherpdg << " " << parton->GetPdgCode();
221 if (mother2 > -1 && agreemother2) {
222 TParticle * m2 = (TParticle*)array->AddrAt(mother2);
226 motherpdg << " " << m2->GetPdgCode();
228 if (oneok && allok && agreemother2) {
229 bool checkok = abs(1 - v_in.M()/v_out.M()) < 0.1;
230 if (!checkok) valid=false;
231 if (!checkok && useStdErr) {
232 // cerr << "Problem: " << i << "[" << fd << "," << ld << "] PDG:" << motherpdg.str() << " ->" << daughterpdg.str() << " Inv.mass: " << v_in.M() << " -> " << v_out.M() << endl;
233 cerr << ListReactionChain(array,i);
235 // cerr << v_in.Px() << " " << v_in.Py() << " " << v_in.Pz() << " " << v_in.E() << endl;
236 } else if (useStdErr) {
237 // cerr << "OK: " << i << "[" << fd << "," << ld << "] PDG:" << motherpdg.str() << " ->" << daughterpdg.str() << " Inv.mass: " << v_in.M() << " -> " << v_out.M() << endl;
241 bool checkok = abs(1 - v_st4.M()/v_st1.M()) < 0.001;
242 if (!checkok) valid=false;
243 if (!checkok && useStdErr) {
244 cerr << " BEAM PARTICLES -> FINAL PARTICLES " << endl;
245 cerr << " status 4 (" << v_st4.M() << ") -> status 1 (" << v_st1.M() << ")" << endl << endl;
253 string THepMCParser::GetParticleName(TParticle * thisPart)
255 TParticlePDG *dbPart = thisPart->GetPDG();
258 name << dbPart->GetName() << "(" << dbPart->PdgCode() << "," << dbPart->Mass() << ")";
260 name << thisPart->GetPdgCode() << "(NoDBinfo)";
265 string THepMCParser::ListReactionChain(TClonesArray * particles, Int_t particleId)
267 ostringstream output;
269 TParticle * part = (TParticle*)particles->AddrAt(particleId);
270 Int_t m1id = part->GetFirstMother();
271 Int_t m2id = part->GetSecondMother();
272 if (m1id > 1) { // ignore the initial collision with beam particles
274 ostringstream outStr;
275 TParticle * m1 = (TParticle*)particles->AddrAt(m1id);
278 inStr << GetParticleName(m1) << "[" << v_in.M() << "]";
280 TParticle * m2 = (TParticle*)particles->AddrAt(m2id);
284 inStr << " " << GetParticleName(m2) << "[" << v_m2.M() << "]";
286 Int_t fd = m1->GetFirstDaughter();
287 Int_t ld = m1->GetLastDaughter();
288 TLorentzVector v_out;
289 part->Momentum(v_out);
290 outStr << GetParticleName(part) << "[" << v_out.M() << "]";
291 for (Int_t i=fd; i<=ld; ++i) {
293 TParticle * d = (TParticle*)particles->AddrAt(i);
297 outStr << " " << GetParticleName(d) << "[" << v_d.M() << "]";
300 output << "Parent reaction, inv mass: " << v_in.M() << " -> " << v_out.M() << endl;
301 output << " - partons: " << inStr.str() << " -> " << outStr.str() << endl;
303 Int_t fd = part->GetFirstDaughter();
304 Int_t ld = part->GetLastDaughter();
307 ostringstream outStr;
309 part->Momentum(v_in);
310 inStr << GetParticleName(part) << "[" << v_in.M() << "]";
312 TParticle * f = (TParticle*)particles->AddrAt(fd);
313 m2id = f->GetSecondMother();
314 if (m2id == particleId) {
315 m2id = f->GetFirstMother();
318 TParticle * m2 = (TParticle*)particles->AddrAt(m2id);
322 inStr << " " << GetParticleName(m2) << "[" << v_m2.M() << "]";
324 TLorentzVector v_out;
326 outStr << GetParticleName(f) << "[" << v_out.M() << "]";
327 for (Int_t i=fd+1; i<=ld; ++i) {
328 TParticle * d = (TParticle*)particles->AddrAt(i);
332 outStr << " " << GetParticleName(d) << "[" << v_d.M() << "]";
334 output << "Child reaction, inv mass: " << v_in.M() << " -> " << v_out.M() << endl;
335 output << " - partons: " << inStr.str() << " -> " << outStr.str() << endl;
337 output << "Child reaction" << endl << " - none" << endl;
344 string THepMCParser::ParseGenEvent2TCloneArray(HepMC::GenEvent * genEvent, TClonesArray * array, bool requireSecondMotherBeforeDaughters)
346 if (requireSecondMotherBeforeDaughters) {
347 return "requireSecondMotherBeforeDaughters not implemented yet!";
350 ostringstream errMsgStream;
351 map<int,Int_t> partonMemory; // unordered_map in c++11 - but probably not much performance gain from that: log(n) vs log(1) where constant can be high
353 // Check event with HepMC's internal validation algorithm
354 if (!genEvent->is_valid()) {
355 errMsgStream << "Error with event id: " << genEvent->event_number() << ", event is not valid!";
356 return errMsgStream.str();
359 // Pull out the beam particles from the event
360 const pair<HepMC::GenParticle *,HepMC::GenParticle *> beamparts = genEvent->beam_particles();
362 // Four sanity checks:
363 // - Beam particles exists and are not the same
364 // - Both beam particles should have no production vertices, they come from the beams
365 // - Both beam particles should have defined end vertices, as they both should contribute
366 // - Both beam particles should have the exact same end vertex
367 if (!beamparts.first || !beamparts.second || beamparts.first->barcode()==beamparts.second->barcode()) {
368 errMsgStream << "Error with event id: " << genEvent->event_number() << ", beam particles doesn't exists or are the same";
369 return errMsgStream.str();
371 if (beamparts.first->production_vertex() || beamparts.second->production_vertex()) {
372 errMsgStream << "Error with event id: " << genEvent->event_number() << ", beam particles have production vertex/vertices...";
373 return errMsgStream.str();
375 if (!beamparts.first->end_vertex() || !beamparts.second->end_vertex()) {
376 errMsgStream << "Error with event id: " << genEvent->event_number() << ", beam particles have undefined end vertex/vertices...";
377 return errMsgStream.str();
379 if (beamparts.first->end_vertex() != beamparts.second->end_vertex()) {
380 errMsgStream << "Error with event id: " << genEvent->event_number() << ", beam particles do not collide in the same end vertex.";
381 return errMsgStream.str();
384 // Set the array to hold the number of particles in the event
385 array->Expand(genEvent->particles_size());
387 // Create a TParticle for each beam particle
388 new((*array)[0]) TParticle(
389 beamparts.first->pdg_id(),
390 beamparts.first->status(), // check if status has the same meaning
393 -1, // first daughter not known yet
394 -1, // last daughter not known yet
395 beamparts.first->momentum().px(),
396 beamparts.first->momentum().py(),
397 beamparts.first->momentum().pz(),
398 beamparts.first->momentum().e(),
399 0, // no production vertex, so zero?
404 partonMemory[beamparts.first->barcode()] = 0;
405 new((*array)[1]) TParticle(
406 beamparts.second->pdg_id(),
407 beamparts.second->status(),
410 -1, // first daughter not known yet
411 -1, // last daughter not known yet
412 beamparts.second->momentum().px(),
413 beamparts.second->momentum().py(),
414 beamparts.second->momentum().pz(),
415 beamparts.second->momentum().e(),
416 0, // no production vertex, so zero?
421 partonMemory[beamparts.second->barcode()] = 1;
422 Int_t arrayID = 2; // start counting IDs after the beam particles
424 // Do first vertex as a special case for performance, since its often has the most daughters and both mothers are known
425 Int_t firstDaughter = arrayID;
426 for (HepMC::GenVertex::particles_out_const_iterator iter = beamparts.first->end_vertex()->particles_out_const_begin();
427 iter != beamparts.first->end_vertex()->particles_out_const_end();
429 new((*array)[arrayID]) TParticle(
432 0, // beam particle 1
433 1, // beam particle 2
434 -1, // first daughter not known yet
435 -1, // last daughter not known yet
436 (*iter)->momentum().px(),
437 (*iter)->momentum().py(),
438 (*iter)->momentum().pz(),
439 (*iter)->momentum().e(),
440 beamparts.first->end_vertex()->position().x(),
441 beamparts.first->end_vertex()->position().y(),
442 beamparts.first->end_vertex()->position().z(),
443 beamparts.first->end_vertex()->position().t()
445 partonMemory[(*iter)->barcode()] = arrayID;
448 Int_t lastDaughter = arrayID-1;
449 ((TParticle*)array->AddrAt(0))->SetFirstDaughter(firstDaughter); // beam particle 1
450 ((TParticle*)array->AddrAt(0))->SetLastDaughter(lastDaughter);
451 ((TParticle*)array->AddrAt(1))->SetFirstDaughter(firstDaughter); // beam particle 2
452 ((TParticle*)array->AddrAt(1))->SetLastDaughter(lastDaughter);
454 // Build vertex list by exploring tree and sorting such that daughters comes after mothers
455 list<HepMC::GenVertex*> vertexList;
456 set<int> vertexSearchSet;
457 ExploreVertex(beamparts.first->end_vertex(),vertexList,vertexSearchSet,requireSecondMotherBeforeDaughters);
459 // Analyze each vertex
460 for (list<HepMC::GenVertex*>::iterator i = vertexList.begin(); i != vertexList.end(); ++i) {
461 HepMC::GenVertex * vertex = (*i);
463 // first establish mother relations
464 HepMC::GenVertex::particles_in_const_iterator iter = vertex->particles_in_const_begin();
465 if (iter == vertex->particles_in_const_end()) {
466 return "Particle without a mother, and its not a beam particle!";
468 int motherA = partonMemory[(*iter)->barcode()];
469 if (((TParticle*)array->AddrAt(motherA))->GetFirstDaughter() > -1) {
470 return "Trying to assign new daughters to a particle that already has daughters defined!";
474 if (iter != vertex->particles_in_const_end()) {
475 motherB = partonMemory[(*iter)->barcode()];
476 if (((TParticle*)array->AddrAt(motherB))->GetFirstDaughter() > -1) {
477 return "Trying to assign new daughters to a particle that already has daughters defined!";
480 if (iter != vertex->particles_in_const_end()) {
481 return "Particle with more than two mothers!";
484 if (motherB > -1 && motherB < motherA) {
485 int swap = motherA; motherA = motherB; motherB = swap;
488 // add the particles to the array, important that they are add in succession with respect to arrayID
489 firstDaughter = arrayID;
490 for (iter = vertex->particles_out_const_begin();
491 iter != vertex->particles_out_const_end();
493 new((*array)[arrayID]) TParticle(
498 -1, // first daughter, if applicable, not known yet
499 -1, // last daughter, if applicable, not known yet
500 (*iter)->momentum().px(),
501 (*iter)->momentum().py(),
502 (*iter)->momentum().pz(),
503 (*iter)->momentum().e(),
504 vertex->position().x(),
505 vertex->position().y(),
506 vertex->position().z(),
507 vertex->position().t()
509 partonMemory[(*iter)->barcode()] = arrayID;
512 lastDaughter = arrayID-1;
513 if (lastDaughter < firstDaughter) {
514 return "Vertex with no out particles, should not be possible!";
516 // update mother with daughter interval
517 ((TParticle*)array->AddrAt(motherA))->SetFirstDaughter(firstDaughter);
518 ((TParticle*)array->AddrAt(motherA))->SetLastDaughter(lastDaughter);
520 ((TParticle*)array->AddrAt(motherB))->SetFirstDaughter(firstDaughter);
521 ((TParticle*)array->AddrAt(motherB))->SetLastDaughter(lastDaughter);
529 void THepMCParser::ExploreVertex(HepMC::GenVertex * vertex, list<HepMC::GenVertex*> & vertexList, set<int> & vertexSearchSet, bool requireSecondMotherBeforeDaughters)
531 for (HepMC::GenVertex::particles_out_const_iterator iter = vertex->particles_out_const_begin();
532 iter != vertex->particles_out_const_end();
534 HepMC::GenVertex * testVertex = (*iter)->end_vertex();
536 bool foundVertex = vertexSearchSet.find((*iter)->end_vertex()->barcode()) != vertexSearchSet.end();
537 if (requireSecondMotherBeforeDaughters) {
538 // redo this algorithem to move subtree instead of node....
539 // its not completely error proof in its current implementation even though the error is extremely rare
541 if (foundVertex) for (list<HepMC::GenVertex*>::iterator i = vertexList.begin(); i != vertexList.end(); ++i) {
542 if ((*i)->barcode() == testVertex->barcode()) {
544 cout << " it happened, the vertex parsing order had to be changed " << endl;
548 vertexSearchSet.insert((*iter)->end_vertex()->barcode());
550 vertexList.push_back(testVertex);
551 if (!foundVertex) ExploreVertex(testVertex,vertexList,vertexSearchSet,requireSecondMotherBeforeDaughters);
555 vertexSearchSet.insert((*iter)->end_vertex()->barcode());
556 vertexList.push_back(testVertex);
557 ExploreVertex(testVertex,vertexList,vertexSearchSet,requireSecondMotherBeforeDaughters);
566 const char * THepMCParser::fgHeavyIonHeaderBranchString = "Ncoll_hard/I:Npart_proj:Npart_targ:Ncoll:spectator_neutrons:spectator_protons:N_Nwounded_collisions:Nwounded_N_collisions:Nwounded_Nwounded_collisions:impact_parameter/F:event_plane_angle:eccentricity:sigma_inel_NN";
567 const char * THepMCParser::fgPdfHeaderBranchString = "id1/I:id2:pdf_id1:pdf_id2:x1/D:x2:scalePDF:pdf1:pdf2";
569 string THepMCParser::ParseGenEvent2HeaderStructs(HepMC::GenEvent * genEvent, HeavyIonHeader_t & heavyIonHeader, PdfHeader_t & pdfHeader, bool fillZeroOnMissingHeavyIon, bool fillZeroOnMissingPdf)
571 HepMC::HeavyIon * heavyIon = genEvent->heavy_ion();
572 HepMC::PdfInfo * pdfInfo = genEvent->pdf_info();
573 if ((!heavyIon && !fillZeroOnMissingHeavyIon) || (!pdfInfo && !fillZeroOnMissingPdf)) {
574 return "HeavyIonInfo and/or PdfInfo not defined for this event, did you read it with IO_GenEvent?";
577 heavyIonHeader.Ncoll_hard = heavyIon->Ncoll_hard();
578 heavyIonHeader.Npart_proj = heavyIon->Npart_proj();
579 heavyIonHeader.Npart_targ = heavyIon->Npart_targ();
580 heavyIonHeader.Ncoll = heavyIon->Ncoll();
581 heavyIonHeader.spectator_neutrons = heavyIon->spectator_neutrons();
582 heavyIonHeader.spectator_protons = heavyIon->spectator_protons();
583 heavyIonHeader.N_Nwounded_collisions = heavyIon->N_Nwounded_collisions();
584 heavyIonHeader.Nwounded_N_collisions = heavyIon->Nwounded_N_collisions();
585 heavyIonHeader.Nwounded_Nwounded_collisions = heavyIon->Nwounded_Nwounded_collisions();
586 heavyIonHeader.impact_parameter = heavyIon->impact_parameter();
587 heavyIonHeader.event_plane_angle = heavyIon->event_plane_angle();
588 heavyIonHeader.eccentricity = heavyIon->eccentricity();
589 heavyIonHeader.sigma_inel_NN = heavyIon->sigma_inel_NN();
591 heavyIonHeader.Ncoll_hard = 0;
592 heavyIonHeader.Npart_proj = 0;
593 heavyIonHeader.Npart_targ = 0;
594 heavyIonHeader.Ncoll = 0;
595 heavyIonHeader.spectator_neutrons = 0;
596 heavyIonHeader.spectator_protons = 0;
597 heavyIonHeader.N_Nwounded_collisions = 0;
598 heavyIonHeader.Nwounded_N_collisions = 0;
599 heavyIonHeader.Nwounded_Nwounded_collisions = 0;
600 heavyIonHeader.impact_parameter = 0.0;
601 heavyIonHeader.event_plane_angle = 0.0;
602 heavyIonHeader.eccentricity = 0.0;
603 heavyIonHeader.sigma_inel_NN = 0.0;
606 pdfHeader.id1 = pdfInfo->id1();
607 pdfHeader.id2 = pdfInfo->id2();
608 pdfHeader.pdf_id1 = pdfInfo->pdf_id1();
609 pdfHeader.pdf_id2 = pdfInfo->pdf_id2();
610 pdfHeader.x1 = pdfInfo->x1();
611 pdfHeader.x2 = pdfInfo->x2();
612 pdfHeader.scalePDF = pdfInfo->scalePDF();
613 pdfHeader.pdf1 = pdfInfo->pdf1();
614 pdfHeader.pdf2 = pdfInfo->pdf2();
618 pdfHeader.pdf_id1 = 0;
619 pdfHeader.pdf_id2 = 0;
622 pdfHeader.scalePDF = 0.0;
623 pdfHeader.pdf1 = 0.0;
624 pdfHeader.pdf2 = 0.0;