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