Moving required CMake version from 2.8.4 to 2.8.8
[u/mrichter/AliRoot.git] / TEvtGen / THepMCParser.cxx
1 // module identifier line...
2 // Author: Brian Thorsbro, 24/6-2014
3
4 #include <iostream>
5 #include <sstream>
6 #include <string>
7 #include <list>
8 #include <set>
9 #include <time.h>
10
11 #include "THepMCParser.h"
12 #include "TObject.h"
13 #include "TTree.h"
14 #include "TClonesArray.h"
15 #include "TFile.h"
16 #include "TParticle.h"
17 #include "TDatabasePDG.h"
18 #include "HepMC/IO_GenEvent.h"
19
20
21 using namespace std;
22
23
24 THepMCParser::THepMCParser(const char * infile) : fTree(0)
25 {
26    HepMC::IO_BaseClass * events = new HepMC::IO_GenEvent(infile, std::ios::in);
27    init(events);
28    delete events;
29    events = 0; // nullptr
30 }
31 THepMCParser::THepMCParser(HepMC::IO_BaseClass * events) : fTree(0)
32 {
33    init(events);
34 }
35 void THepMCParser::init(HepMC::IO_BaseClass * events)
36 {
37    int particlecount = 0;
38    fTree = new TTree("treeEPOS","Tree EPOS");
39    TClonesArray * array = new TClonesArray("TParticle");
40    // array->BypassStreamer();
41    fTree->Branch("Particles",&array); // more flags?
42    THepMCParser::HeavyIonHeader_t heavyIonHeader;
43    fTree->Branch("HeavyIonInfo", &heavyIonHeader, THepMCParser::fgHeavyIonHeaderBranchString);
44    THepMCParser::PdfHeader_t pdfHeader;
45    fTree->Branch("PdfInfo", &pdfHeader, THepMCParser::fgPdfHeaderBranchString);
46    HepMC::GenEvent* evt =  0; // nullptr
47    while ((evt = events->read_next_event())) {
48       string errMsg1 = ParseGenEvent2TCloneArray(evt,array);
49       string errMsg2 = ParseGenEvent2HeaderStructs(evt,heavyIonHeader,pdfHeader);
50       if (errMsg1.length() == 0 && errMsg2.length() == 0) {
51          fTree->Fill();
52       } else {
53          if (errMsg1.length() != 0) cerr << errMsg1 << endl;
54          if (errMsg2.length() != 0) cerr << errMsg2 << endl;
55       }
56       particlecount += array->Capacity();
57    }
58 //   array->Clear();
59 //   delete array;
60 //   array = 0; // nullptr
61    cout << " parsed " << particlecount << " particles" << endl;
62 }
63
64
65 TTree * THepMCParser::GetTTree()
66 {
67    return fTree;
68 }
69 void THepMCParser::WriteTTreeToFile(const char *outfile)
70 {
71    TFile * f = new TFile(outfile, "recreate");
72    fTree->Write();
73    delete f;
74    f = 0; // nullptr
75 }
76
77
78
79 // Based on a validator written by Peter Hristov, CERN
80 bool THepMCParser::IsValidMotherDaughtersConsitency(bool useStdErr, bool requireSecondMotherBeforeDaughters)
81 {
82    bool valid = true;
83    TClonesArray * array = new TClonesArray("TParticle");
84    TBranch* branch = fTree->GetBranch("Particles");
85    branch->SetAddress(&array);
86    Int_t count = branch->GetEntries();
87    for (Int_t idx=0; idx<count; ++idx) {
88       array->Clear();
89       branch->GetEntry(idx); // "fill" the array
90       Int_t nkeep = array->GetEntriesFast();
91       for (Int_t i=0; i<nkeep; i++) {
92          TParticle * part = (TParticle*)array->AddrAt(i);
93          Int_t mum1 = part->GetFirstMother();
94          Int_t mum2 = part->GetSecondMother();
95          Int_t fd = part->GetFirstDaughter();
96          Int_t ld = part->GetLastDaughter();
97          if (mum1>-1 && i<mum1) {
98             valid = false;
99             if (useStdErr) cerr << "Problem: first_mother(" << mum1 << ") after daughter(" << i << ")" << endl;
100          }
101          if (mum2>-1 && i<mum2 && requireSecondMotherBeforeDaughters) {
102             valid = false;
103             if (useStdErr) cerr << "Problem: second_mother(" << mum2 << ") after daughter(" << i << ")" << endl;
104          }
105          if (fd > ld ) {
106             valid = false;
107             if (useStdErr) cerr << "Problem: first_daughter(" << fd << ") > last_daughter(" << ld << ")" << endl;
108          }
109          for (Int_t id=TMath::Max(fd,0); id<=ld; id++) {
110             TParticle * daughter = (TParticle*)array->AddrAt(id);
111             if (daughter->GetFirstMother() != i && daughter->GetSecondMother() != i) {
112                valid = false;
113                if (useStdErr) cerr << "Problem: mother("<< i << ") not registered as mother for either first_daughter("
114                      << daughter->GetFirstMother() << ") or second_daughter("
115                      << daughter->GetSecondMother() << ")" << endl;
116             }
117          }
118       }
119    }
120    delete array;
121    array = 0;
122    return valid;
123 }
124
125 bool THepMCParser::IsValidParticleInvariantMass(bool useStdErr, bool includeStatusCode2Particles)
126 {
127    bool valid = true;
128    TClonesArray *array = new TClonesArray("TParticle");
129    TBranch* branch = fTree->GetBranch("Particles");
130    branch->SetAddress(&array);
131    Int_t count = branch->GetEntries();
132    for (Int_t idx=0; idx<count; ++idx) {
133       array->Clear();
134       branch->GetEntry(idx);
135       Int_t nkeep = array->GetEntries();
136       for (Int_t i=0; i<nkeep; i++) {
137          TParticle * parton = (TParticle*)array->AddrAt(i);
138          if (parton->GetStatusCode()==2 && !includeStatusCode2Particles) {
139             continue;
140          }
141          TLorentzVector v;
142          parton->Momentum(v);
143          Double_t m1 = v.M(); // invariant mass from the particle in the HepMC event
144          TParticlePDG *dbParton = parton->GetPDG();
145          if (!dbParton) {
146             if (useStdErr) cerr << "Warning: could not look up particle with PDG Code: " << parton->GetPdgCode() << endl << endl;
147             continue;
148          }
149          Double_t m2 = dbParton->Mass();
150          bool checkok;
151          if (m2 == 0) {
152             checkok = abs(m1) < 0.0001; // no such thing as negative mass...
153          } else {
154             checkok = abs(1 - m1/m2) < 0.01;
155          }
156          if (!checkok && useStdErr) {
157             cerr << "Problem: " << GetParticleName(parton) << " HepMC:" << m1 << endl;
158             cerr << ListReactionChain(array,i);
159             cerr << endl;
160          }
161          if (!checkok)
162             valid = false;
163       }
164    }
165    delete array;
166    array = 0;
167    return valid;
168 }
169
170 bool THepMCParser::IsValidVertexInvariantMass(bool useStdErr, bool includeStatusCode2Particles)
171 {
172    bool valid = true;
173    TClonesArray * array = new TClonesArray("TParticle");
174    TBranch* branch = fTree->GetBranch("Particles");
175    branch->SetAddress(&array);
176    Int_t count = branch->GetEntries();
177    for (Int_t idx=0; idx<count; ++idx) {
178       array->Clear();
179       branch->GetEntry(idx); // "fill" the array
180       TLorentzVector v_st1;
181       TLorentzVector v_st4;
182       Int_t nkeep = array->GetEntriesFast();
183       for (Int_t i=0; i<nkeep; i++) {
184          TParticle * parton = (TParticle*)array->AddrAt(i);
185          TLorentzVector v_in;
186          parton->Momentum(v_in);
187          if (parton->GetStatusCode()==4) {
188             v_st4 += v_in;
189          } else if (parton->GetStatusCode()==1) {
190             v_st1 += v_in;
191          }
192          if (!includeStatusCode2Particles) { // only check beam particles vs final particles
193             continue;
194          }
195          Int_t fd = parton->GetFirstDaughter();
196          Int_t ld = parton->GetLastDaughter();
197          if (fd == -1) continue; // no daughters, continue loop
198          Int_t mother2 = -1;
199          TLorentzVector v_out;
200          bool oneok = false; bool allok = true; bool agreemother2 = true;
201          ostringstream daughterpdg;
202          ostringstream motherpdg;
203          for (Int_t id=TMath::Max(fd,0); id<=ld; id++) {
204             TParticle * daughter = (TParticle*)array->AddrAt(id);
205             if (fd==id) {
206                daughter->Momentum(v_out);
207                mother2 = daughter->GetSecondMother();
208             } else {
209                TLorentzVector d;
210                daughter->Momentum(d);
211                v_out += d;
212                if (daughter->GetSecondMother() != mother2) agreemother2 = false;
213             }
214             if (daughter->GetFirstMother() == i) {
215                oneok = true;
216             } else {
217                allok = false;
218             }
219             daughterpdg << " " << daughter->GetPdgCode();
220          }
221          motherpdg << " " << parton->GetPdgCode();
222          if (mother2 > -1 && agreemother2) {
223             TParticle * m2 = (TParticle*)array->AddrAt(mother2);
224             TLorentzVector m2v;
225             m2->Momentum(m2v);
226             v_in += m2v;
227             motherpdg << " " << m2->GetPdgCode();
228          }
229          if (oneok && allok && agreemother2) {
230             bool checkok = abs(1 - v_in.M()/v_out.M()) < 0.1;
231             if (!checkok) valid=false;
232             if (!checkok && useStdErr) {
233                //             cerr << "Problem: " << i << "[" << fd << "," << ld << "] PDG:" << motherpdg.str() << " ->" << daughterpdg.str() << "  Inv.mass: " << v_in.M() << " -> " << v_out.M() << endl;
234                cerr << ListReactionChain(array,i);
235                cerr << endl;
236                //             cerr << v_in.Px() << " " << v_in.Py() << " " << v_in.Pz() << " " << v_in.E() << endl;
237             } else if (useStdErr) {
238                //             cerr << "OK: " << i << "[" << fd << "," << ld << "] PDG:" << motherpdg.str() << " ->" << daughterpdg.str() << "  Inv.mass: " << v_in.M() << " -> " << v_out.M() << endl;
239             }
240          }
241       }
242       bool checkok = abs(1 - v_st4.M()/v_st1.M()) < 0.001;
243       if (!checkok) valid=false;
244       if (!checkok && useStdErr) {
245          cerr << " BEAM PARTICLES -> FINAL PARTICLES " << endl;
246          cerr << " status 4 (" << v_st4.M() << ") -> status 1 (" << v_st1.M() << ")" << endl << endl;
247       }
248    }
249    delete array;
250    array = 0;
251    return valid;
252 }
253
254 string THepMCParser::GetParticleName(TParticle * thisPart)
255 {
256    TParticlePDG *dbPart = thisPart->GetPDG();
257    ostringstream name;
258    if (dbPart) {
259       name << dbPart->GetName() << "(" << dbPart->PdgCode() << "," << dbPart->Mass() << ")";
260    } else {
261       name << thisPart->GetPdgCode() << "(NoDBinfo)";
262    }
263    return name.str();
264 }
265
266 string THepMCParser::ListReactionChain(TClonesArray * particles, Int_t particleId)
267 {
268    ostringstream output;
269
270    TParticle * part = (TParticle*)particles->AddrAt(particleId);
271    Int_t m1id = part->GetFirstMother();
272    Int_t m2id = part->GetSecondMother();
273    if (m1id > 1) { // ignore the initial collision with beam particles
274       ostringstream inStr;
275       ostringstream outStr;
276       TParticle * m1 = (TParticle*)particles->AddrAt(m1id);
277       TLorentzVector v_in;
278       m1->Momentum(v_in);
279       inStr << GetParticleName(m1) << "[" << v_in.M() << "]";
280       if (m2id > 1) {
281          TParticle * m2 = (TParticle*)particles->AddrAt(m2id);
282          TLorentzVector v_m2;
283          m2->Momentum(v_m2);
284          v_in += v_m2;
285          inStr << " " << GetParticleName(m2) << "[" << v_m2.M() << "]";
286       }
287       Int_t fd = m1->GetFirstDaughter();
288       Int_t ld = m1->GetLastDaughter();
289       TLorentzVector v_out;
290       part->Momentum(v_out);
291       outStr << GetParticleName(part) << "[" << v_out.M() << "]";
292       for (Int_t i=fd; i<=ld; ++i) {
293          if (i!=particleId) {
294             TParticle * d = (TParticle*)particles->AddrAt(i);
295             TLorentzVector v_d;
296             d->Momentum(v_d);
297             v_out += v_d;
298             outStr << " " << GetParticleName(d) << "[" << v_d.M() << "]";
299          }
300       }
301       output << "Parent reaction, inv mass: " << v_in.M() << " -> " << v_out.M() << endl;
302       output << " - partons: " << inStr.str() << " -> " << outStr.str() << endl;
303    }
304    Int_t fd = part->GetFirstDaughter();
305    Int_t ld = part->GetLastDaughter();
306    if (fd > -1) {
307       ostringstream inStr;
308       ostringstream outStr;
309       TLorentzVector v_in;
310       part->Momentum(v_in);
311       inStr << GetParticleName(part) << "[" << v_in.M() << "]";
312
313       TParticle * f = (TParticle*)particles->AddrAt(fd);
314       m2id = f->GetSecondMother();
315       if (m2id == particleId) {
316          m2id = f->GetFirstMother();
317       }
318       if (m2id > -1) {
319          TParticle * m2 = (TParticle*)particles->AddrAt(m2id);
320          TLorentzVector v_m2;
321          m2->Momentum(v_m2);
322          v_in += v_m2;
323          inStr << " " << GetParticleName(m2) << "[" << v_m2.M() << "]";
324       }
325       TLorentzVector v_out;
326       f->Momentum(v_out);
327       outStr << GetParticleName(f) << "[" << v_out.M() << "]";
328       for (Int_t i=fd+1; i<=ld; ++i) {
329          TParticle * d = (TParticle*)particles->AddrAt(i);
330          TLorentzVector v_d;
331          d->Momentum(v_d);
332          v_out += v_d;
333          outStr << " " << GetParticleName(d) << "[" << v_d.M() << "]";
334       }
335       output << "Child reaction, inv mass: " << v_in.M() << " -> " << v_out.M() << endl;
336       output << " - partons: " << inStr.str() << " -> " << outStr.str() << endl;
337    } else {
338       output << "Child reaction" << endl << " - none" << endl;
339    }
340
341    return output.str();
342 }
343
344
345 string THepMCParser::ParseGenEvent2TCloneArray(HepMC::GenEvent * genEvent, TClonesArray * array, string momUnit, string lenUnit, bool requireSecondMotherBeforeDaughters)
346 {
347    if (requireSecondMotherBeforeDaughters) {
348       return "requireSecondMotherBeforeDaughters not implemented yet!";
349    }
350    genEvent->use_units(momUnit, lenUnit);
351    array->Clear();
352    ostringstream errMsgStream;
353    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
354
355    // Check event with HepMC's internal validation algorithm
356    if (!genEvent->is_valid()) {
357       errMsgStream << "Error with event id: " << genEvent->event_number() << ", event is not valid!";
358       return errMsgStream.str();
359    }
360
361    // Pull out the beam particles from the event
362    const pair<HepMC::GenParticle *,HepMC::GenParticle *> beamparts = genEvent->beam_particles();
363
364    // Four sanity checks:
365    // - Beam particles exists and are not the same
366    // - Both beam particles should have no production vertices, they come from the beams
367    // - Both beam particles should have defined end vertices, as they both should contribute
368    // - Both beam particles should have the exact same end vertex
369    if (!beamparts.first || !beamparts.second || beamparts.first->barcode()==beamparts.second->barcode()) {
370       errMsgStream << "Error with event id: " << genEvent->event_number() << ", beam particles doesn't exists or are the same";
371       return errMsgStream.str();
372    }
373    if (beamparts.first->production_vertex() || beamparts.second->production_vertex()) {
374       errMsgStream << "Error with event id: " << genEvent->event_number() << ", beam particles have production vertex/vertices...";
375       return errMsgStream.str();
376    }
377    if (!beamparts.first->end_vertex() || !beamparts.second->end_vertex()) {
378       errMsgStream << "Error with event id: " << genEvent->event_number() << ", beam particles have undefined end vertex/vertices...";
379       return errMsgStream.str();
380    }
381    if (beamparts.first->end_vertex() != beamparts.second->end_vertex()) {
382       errMsgStream << "Error with event id: " << genEvent->event_number() << ", beam particles do not collide in the same end vertex.";
383       return errMsgStream.str();
384    }
385
386    // Set the array to hold the number of particles in the event
387    array->Expand(genEvent->particles_size());
388
389    // Create a TParticle for each beam particle
390    new((*array)[0]) TParticle(
391          beamparts.first->pdg_id(),
392          beamparts.first->status(), // check if status has the same meaning
393          -1, // no mother1
394          -1, // no mother2
395          -1, // first daughter not known yet
396          -1, // last daughter not known yet
397          beamparts.first->momentum().px(),
398          beamparts.first->momentum().py(),
399          beamparts.first->momentum().pz(),
400          beamparts.first->momentum().e(),
401          0, // no production vertex, so zero?
402          0,
403          0,
404          0
405    );
406    partonMemory[beamparts.first->barcode()] = 0;
407    new((*array)[1]) TParticle(
408          beamparts.second->pdg_id(),
409          beamparts.second->status(),
410          -1, // no mother1
411          -1, // no mother2
412          -1, // first daughter not known yet
413          -1, // last daughter not known yet
414          beamparts.second->momentum().px(),
415          beamparts.second->momentum().py(),
416          beamparts.second->momentum().pz(),
417          beamparts.second->momentum().e(),
418          0, // no production vertex, so zero?
419          0,
420          0,
421          0
422    );
423    partonMemory[beamparts.second->barcode()] = 1;
424    Int_t arrayID = 2; // start counting IDs after the beam particles
425
426    // Do first vertex as a special case for performance, since its often has the most daughters and both mothers are known
427    Int_t firstDaughter = arrayID;
428    for (HepMC::GenVertex::particles_out_const_iterator iter = beamparts.first->end_vertex()->particles_out_const_begin();
429          iter != beamparts.first->end_vertex()->particles_out_const_end();
430          ++iter) {
431       new((*array)[arrayID]) TParticle(
432             (*iter)->pdg_id(),
433             (*iter)->status(),
434             0, // beam particle 1
435             1, // beam particle 2
436             -1, // first daughter not known yet
437             -1, // last daughter not known yet
438             (*iter)->momentum().px(),
439             (*iter)->momentum().py(),
440             (*iter)->momentum().pz(),
441             (*iter)->momentum().e(),
442             beamparts.first->end_vertex()->position().x(),
443             beamparts.first->end_vertex()->position().y(),
444             beamparts.first->end_vertex()->position().z(),
445             beamparts.first->end_vertex()->position().t()
446       );
447       partonMemory[(*iter)->barcode()] = arrayID;
448       ++arrayID;
449    }
450    Int_t lastDaughter = arrayID-1;
451    ((TParticle*)array->AddrAt(0))->SetFirstDaughter(firstDaughter); // beam particle 1
452    ((TParticle*)array->AddrAt(0))->SetLastDaughter(lastDaughter);
453    ((TParticle*)array->AddrAt(1))->SetFirstDaughter(firstDaughter); // beam particle 2
454    ((TParticle*)array->AddrAt(1))->SetLastDaughter(lastDaughter);
455
456    // Build vertex list by exploring tree and sorting such that daughters comes after mothers
457    list<HepMC::GenVertex*> vertexList;
458    set<int> vertexSearchSet;
459    ExploreVertex(beamparts.first->end_vertex(),vertexList,vertexSearchSet,requireSecondMotherBeforeDaughters);
460
461    // Analyze each vertex
462    for (list<HepMC::GenVertex*>::iterator i = vertexList.begin(); i != vertexList.end(); ++i) {
463       HepMC::GenVertex * vertex = (*i);
464
465       // first establish mother relations
466       HepMC::GenVertex::particles_in_const_iterator iter = vertex->particles_in_const_begin();
467       if (iter == vertex->particles_in_const_end()) {
468          return "Particle without a mother, and its not a beam particle!";
469       }
470       int motherA = partonMemory[(*iter)->barcode()];
471       if (((TParticle*)array->AddrAt(motherA))->GetFirstDaughter() > -1) {
472          return "Trying to assign new daughters to a particle that already has daughters defined!";
473       }
474       ++iter;
475       int motherB = -1;
476       if (iter != vertex->particles_in_const_end()) {
477          motherB = partonMemory[(*iter)->barcode()];
478          if (((TParticle*)array->AddrAt(motherB))->GetFirstDaughter() > -1) {
479             return "Trying to assign new daughters to a particle that already has daughters defined!";
480          }
481          ++iter;
482          if (iter != vertex->particles_in_const_end()) {
483             return "Particle with more than two mothers!";
484          }
485       }
486       if (motherB > -1 && motherB < motherA) {
487          int swap = motherA; motherA = motherB; motherB = swap;
488       }
489
490       // add the particles to the array, important that they are add in succession with respect to arrayID
491       firstDaughter = arrayID;
492       for (iter = vertex->particles_out_const_begin();
493             iter != vertex->particles_out_const_end();
494             ++iter) {
495          new((*array)[arrayID]) TParticle(
496                (*iter)->pdg_id(),
497                (*iter)->status(),
498                motherA, // mother 1
499                motherB, // mother 2
500                -1, // first daughter, if applicable, not known yet
501                -1, // last daughter, if applicable, not known yet
502                (*iter)->momentum().px(),
503                (*iter)->momentum().py(),
504                (*iter)->momentum().pz(),
505                (*iter)->momentum().e(),
506                vertex->position().x(),
507                vertex->position().y(),
508                vertex->position().z(),
509                vertex->position().t()
510          );
511          partonMemory[(*iter)->barcode()] = arrayID;
512          ++arrayID;
513       }
514       lastDaughter = arrayID-1;
515       if (lastDaughter < firstDaughter) {
516          return "Vertex with no out particles, should not be possible!";
517       }
518       // update mother with daughter interval
519       ((TParticle*)array->AddrAt(motherA))->SetFirstDaughter(firstDaughter);
520       ((TParticle*)array->AddrAt(motherA))->SetLastDaughter(lastDaughter);
521       if (motherB > -1) {
522          ((TParticle*)array->AddrAt(motherB))->SetFirstDaughter(firstDaughter);
523          ((TParticle*)array->AddrAt(motherB))->SetLastDaughter(lastDaughter);
524       }
525    }
526
527    return "";
528 }
529
530
531 void THepMCParser::ExploreVertex(HepMC::GenVertex * vertex, list<HepMC::GenVertex*> & vertexList, set<int> & vertexSearchSet, bool requireSecondMotherBeforeDaughters)
532 {
533    for (HepMC::GenVertex::particles_out_const_iterator iter = vertex->particles_out_const_begin();
534          iter != vertex->particles_out_const_end();
535          ++iter) {
536       HepMC::GenVertex * testVertex = (*iter)->end_vertex();
537       if (testVertex) {
538          bool foundVertex = vertexSearchSet.find((*iter)->end_vertex()->barcode()) != vertexSearchSet.end();
539          if (requireSecondMotherBeforeDaughters) {
540             // redo this algorithem to move subtree instead of node....
541             // its not completely error proof in its current implementation even though the error is extremely rare
542
543             if (foundVertex) for (list<HepMC::GenVertex*>::iterator i = vertexList.begin(); i != vertexList.end(); ++i) {
544                if ((*i)->barcode() == testVertex->barcode()) {
545                   vertexList.erase(i);
546                   cout << " it happened, the vertex parsing order had to be changed " << endl;
547                   break;
548                }
549             } else {
550                vertexSearchSet.insert((*iter)->end_vertex()->barcode());
551             }
552             vertexList.push_back(testVertex);
553             if (!foundVertex) ExploreVertex(testVertex,vertexList,vertexSearchSet,requireSecondMotherBeforeDaughters);
554
555          } else {
556             if (!foundVertex) {
557                vertexSearchSet.insert((*iter)->end_vertex()->barcode());
558                vertexList.push_back(testVertex);
559                ExploreVertex(testVertex,vertexList,vertexSearchSet,requireSecondMotherBeforeDaughters);
560             }
561          }
562       }
563    }
564 }
565
566
567
568 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";
569 const char * THepMCParser::fgPdfHeaderBranchString = "id1/I:id2:pdf_id1:pdf_id2:x1/D:x2:scalePDF:pdf1:pdf2";
570
571 string THepMCParser::ParseGenEvent2HeaderStructs(HepMC::GenEvent * genEvent, HeavyIonHeader_t & heavyIonHeader, PdfHeader_t & pdfHeader, bool fillZeroOnMissingHeavyIon, bool fillZeroOnMissingPdf)
572 {
573    HepMC::HeavyIon * heavyIon = genEvent->heavy_ion();
574    HepMC::PdfInfo * pdfInfo = genEvent->pdf_info();
575    if ((!heavyIon && !fillZeroOnMissingHeavyIon) || (!pdfInfo && !fillZeroOnMissingPdf)) {
576       return "HeavyIonInfo and/or PdfInfo not defined for this event, did you read it with IO_GenEvent?";
577    }
578    if (heavyIon) {
579       heavyIonHeader.Ncoll_hard = heavyIon->Ncoll_hard();
580       heavyIonHeader.Npart_proj = heavyIon->Npart_proj();
581       heavyIonHeader.Npart_targ = heavyIon->Npart_targ();
582       heavyIonHeader.Ncoll = heavyIon->Ncoll();
583       heavyIonHeader.spectator_neutrons = heavyIon->spectator_neutrons();
584       heavyIonHeader.spectator_protons = heavyIon->spectator_protons();
585       heavyIonHeader.N_Nwounded_collisions = heavyIon->N_Nwounded_collisions();
586       heavyIonHeader.Nwounded_N_collisions = heavyIon->Nwounded_N_collisions();
587       heavyIonHeader.Nwounded_Nwounded_collisions = heavyIon->Nwounded_Nwounded_collisions();
588       heavyIonHeader.impact_parameter = heavyIon->impact_parameter();
589       heavyIonHeader.event_plane_angle = heavyIon->event_plane_angle();
590       heavyIonHeader.eccentricity = heavyIon->eccentricity();
591       heavyIonHeader.sigma_inel_NN = heavyIon->sigma_inel_NN();
592    } else {
593       heavyIonHeader.Ncoll_hard = 0;
594       heavyIonHeader.Npart_proj = 0;
595       heavyIonHeader.Npart_targ = 0;
596       heavyIonHeader.Ncoll = 0;
597       heavyIonHeader.spectator_neutrons = 0;
598       heavyIonHeader.spectator_protons = 0;
599       heavyIonHeader.N_Nwounded_collisions = 0;
600       heavyIonHeader.Nwounded_N_collisions = 0;
601       heavyIonHeader.Nwounded_Nwounded_collisions = 0;
602       heavyIonHeader.impact_parameter = 0.0;
603       heavyIonHeader.event_plane_angle = 0.0;
604       heavyIonHeader.eccentricity = 0.0;
605       heavyIonHeader.sigma_inel_NN = 0.0;
606    }
607    if (pdfInfo) {
608       pdfHeader.id1 = pdfInfo->id1();
609       pdfHeader.id2 = pdfInfo->id2();
610       pdfHeader.pdf_id1 = pdfInfo->pdf_id1();
611       pdfHeader.pdf_id2 = pdfInfo->pdf_id2();
612       pdfHeader.x1 = pdfInfo->x1();
613       pdfHeader.x2 = pdfInfo->x2();
614       pdfHeader.scalePDF = pdfInfo->scalePDF();
615       pdfHeader.pdf1 = pdfInfo->pdf1();
616       pdfHeader.pdf2 = pdfInfo->pdf2();
617    } else {
618       pdfHeader.id1 = 0;
619       pdfHeader.id2 = 0;
620       pdfHeader.pdf_id1 = 0;
621       pdfHeader.pdf_id2 = 0;
622       pdfHeader.x1 = 0.0;
623       pdfHeader.x2 = 0.0;
624       pdfHeader.scalePDF = 0.0;
625       pdfHeader.pdf1 = 0.0;
626       pdfHeader.pdf2 = 0.0;
627    }
628    return "";
629 }
630
631
632
633
634
635
636