Added HepMC parser. Build it to to a library independent of HepMC. Updated HepMC...
authorbthorsbr <brian.peter.thorsbro@cern.ch>
Mon, 14 Jul 2014 12:21:30 +0000 (14:21 +0200)
committerhristov <Peter.Hristov@cern.ch>
Fri, 8 Aug 2014 14:37:53 +0000 (16:37 +0200)
TEvtGen/CMakelibHepMC.pkg
TEvtGen/CMakelibHepMCParser.pkg [new file with mode: 0644]
TEvtGen/HepMC/HEPEVT_Wrapper.h
TEvtGen/THepMCParser.cxx [new file with mode: 0644]
TEvtGen/THepMCParser.h [new file with mode: 0644]
TEvtGen/THepMCParserLinkDef.h [new file with mode: 0644]

index 4fb0887..79426df 100644 (file)
@@ -44,6 +44,7 @@ HepMC/IO_HERWIG.cc
 HepMC/PdfInfo.cc
 HepMC/Polarization.cc
 HepMC/SearchVector.cc
+HepMC/SimpleVector.icc
 HepMC/StreamHelpers.cc
 HepMC/StreamInfo.cc
 HepMC/WeightContainer.cc
@@ -79,7 +80,6 @@ HepMC/PythiaWrapper6_4_WIN32.h
 HepMC/PythiaWrapper.h
 HepMC/SearchVector.h
 HepMC/SimpleVector.h
-HepMC/SimpleVector.icc
 HepMC/StreamHelpers.h
 HepMC/StreamInfo.h
 HepMC/TempParticleMap.h
diff --git a/TEvtGen/CMakelibHepMCParser.pkg b/TEvtGen/CMakelibHepMCParser.pkg
new file mode 100644 (file)
index 0000000..aefff43
--- /dev/null
@@ -0,0 +1,38 @@
+#--------------------------------------------------------------------------------#
+# Package File for EvtGen                                                        #
+# Author : Johny Jose (johny.jose@cern.ch)                                       #
+# Variables Defined :                                                            #
+#                                                                                #
+# SRCS - C++ source files                                                        #
+# HDRS - C++ header files                                                        #
+# DHDR - ROOT Dictionary Linkdef header file                                     #
+# CSRCS - C source files                                                         #
+# CHDRS - C header files                                                         #
+# EINCLUDE - Include directories                                                 #
+# EDEFINE - Compiler definitions                                                 #
+# ELIBS - Extra libraries to link                                                #
+# ELIBSDIR - Extra library directories                                           #
+# PACKFFLAGS - Fortran compiler flags for package                                #
+# PACKCXXFLAGS - C++ compiler flags for package                                  #
+# PACKCFLAGS - C compiler flags for package                                      #
+# PACKSOFLAGS - Shared library linking flags                                     #
+# PACKLDFLAGS - Module linker flags                                              #
+# PACKBLIBS - Libraries to link (Executables only)                               #
+# EXPORT - Header files to be exported                                           #
+# CINTHDRS - Dictionary header files                                             #
+# CINTAUTOLINK - Set automatic dictionary generation                             #
+# ARLIBS - Archive Libraries and objects for linking (Executables only)          #
+# SHLIBS - Shared Libraries and objects for linking (Executables only)           #
+#--------------------------------------------------------------------------------#
+
+set ( SRCS  THepMCParser.cxx)
+
+set ( HDRS THepMCParser.h)
+
+set ( DHDR  THepMCParserLinkDef.h)
+
+set ( EINCLUDE TEvtGen TEvtGen/HepMC) 
+
+
+
+
index cc56250..99e8d36 100644 (file)
@@ -1,7 +1,7 @@
 //--------------------------------------------------------------------------
 
 #ifndef HEPEVT_EntriesAllocation
-#define HEPEVT_EntriesAllocation 10000
+#define HEPEVT_EntriesAllocation 20000
 #endif  // HEPEVT_EntriesAllocation
 
 //--------------------------------------------------------------------------
diff --git a/TEvtGen/THepMCParser.cxx b/TEvtGen/THepMCParser.cxx
new file mode 100644 (file)
index 0000000..2d18d72
--- /dev/null
@@ -0,0 +1,634 @@
+// module identifier line...
+// Author: Brian Thorsbro, 24/6-2014
+
+#include <iostream>
+#include <sstream>
+#include <string>
+#include <list>
+#include <set>
+#include <time.h>
+
+#include "THepMCParser.h"
+#include "HepMC/IO_GenEvent.h"
+#include "TObject.h"
+#include "TTree.h"
+#include "TClonesArray.h"
+#include "TFile.h"
+#include "TParticle.h"
+#include "TDatabasePDG.h"
+
+using namespace std;
+
+
+THepMCParser::THepMCParser(const char * infile) : fTree(0)
+{
+   HepMC::IO_BaseClass * events = new HepMC::IO_GenEvent(infile, std::ios::in);
+   init(events);
+   delete events;
+   events = 0; // nullptr
+}
+THepMCParser::THepMCParser(HepMC::IO_BaseClass * events) : fTree(0)
+{
+   init(events);
+}
+void THepMCParser::init(HepMC::IO_BaseClass * events)
+{
+   int particlecount = 0;
+   fTree = new TTree("treeEPOS","Tree EPOS");
+   TClonesArray * array = new TClonesArray("TParticle");
+   // array->BypassStreamer();
+   fTree->Branch("Particles",&array); // more flags?
+   THepMCParser::HeavyIonHeader_t heavyIonHeader;
+   fTree->Branch("HeavyIonInfo", &heavyIonHeader, THepMCParser::fgHeavyIonHeaderBranchString);
+   THepMCParser::PdfHeader_t pdfHeader;
+   fTree->Branch("PdfInfo", &pdfHeader, THepMCParser::fgPdfHeaderBranchString);
+   HepMC::GenEvent* evt =  0; // nullptr
+   while ((evt = events->read_next_event())) {
+      string errMsg1 = ParseGenEvent2TCloneArray(evt,array);
+      string errMsg2 = ParseGenEvent2HeaderStructs(evt,heavyIonHeader,pdfHeader);
+      if (errMsg1.length() == 0 && errMsg2.length() == 0) {
+         fTree->Fill();
+      } else {
+         if (errMsg1.length() != 0) cerr << errMsg1 << endl;
+         if (errMsg2.length() != 0) cerr << errMsg2 << endl;
+      }
+      particlecount += array->Capacity();
+   }
+//   array->Clear();
+//   delete array;
+//   array = 0; // nullptr
+   cout << " parsed " << particlecount << " particles" << endl;
+}
+
+
+TTree * THepMCParser::GetTTree()
+{
+   return fTree;
+}
+void THepMCParser::WriteTTreeToFile(const char *outfile)
+{
+   TFile * f = new TFile(outfile, "recreate");
+   fTree->Write();
+   delete f;
+   f = 0; // nullptr
+}
+
+
+
+// Based on a validator written by Peter Hristov, CERN
+bool THepMCParser::IsValidMotherDaughtersConsitency(bool useStdErr, bool requireSecondMotherBeforeDaughters)
+{
+   bool valid = true;
+   TClonesArray * array = new TClonesArray("TParticle");
+   TBranch* branch = fTree->GetBranch("Particles");
+   branch->SetAddress(&array);
+   Int_t count = branch->GetEntries();
+   for (Int_t idx=0; idx<count; ++idx) {
+      array->Clear();
+      branch->GetEntry(idx); // "fill" the array
+      Int_t nkeep = array->GetEntriesFast();
+      for (Int_t i=0; i<nkeep; i++) {
+         TParticle * part = (TParticle*)array->AddrAt(i);
+         Int_t mum1 = part->GetFirstMother();
+         Int_t mum2 = part->GetSecondMother();
+         Int_t fd = part->GetFirstDaughter();
+         Int_t ld = part->GetLastDaughter();
+         if (mum1>-1 && i<mum1) {
+            valid = false;
+            if (useStdErr) cerr << "Problem: first_mother(" << mum1 << ") after daughter(" << i << ")" << endl;
+         }
+         if (mum2>-1 && i<mum2 && requireSecondMotherBeforeDaughters) {
+            valid = false;
+            if (useStdErr) cerr << "Problem: second_mother(" << mum2 << ") after daughter(" << i << ")" << endl;
+         }
+         if (fd > ld ) {
+            valid = false;
+            if (useStdErr) cerr << "Problem: first_daughter(" << fd << ") > last_daughter(" << ld << ")" << endl;
+         }
+         for (Int_t id=TMath::Max(fd,0); id<=ld; id++) {
+            TParticle * daughter = (TParticle*)array->AddrAt(id);
+            if (daughter->GetFirstMother() != i && daughter->GetSecondMother() != i) {
+               valid = false;
+               if (useStdErr) cerr << "Problem: mother("<< i << ") not registered as mother for either first_daughter("
+                     << daughter->GetFirstMother() << ") or second_daughter("
+                     << daughter->GetSecondMother() << ")" << endl;
+            }
+         }
+      }
+   }
+   delete array;
+   array = 0;
+   return valid;
+}
+
+bool THepMCParser::IsValidParticleInvariantMass(bool useStdErr, bool includeStatusCode2Particles)
+{
+   bool valid = true;
+   TClonesArray *array = new TClonesArray("TParticle");
+   TBranch* branch = fTree->GetBranch("Particles");
+   branch->SetAddress(&array);
+   Int_t count = branch->GetEntries();
+   for (Int_t idx=0; idx<count; ++idx) {
+      array->Clear();
+      branch->GetEntry(idx);
+      Int_t nkeep = array->GetEntries();
+      for (Int_t i=0; i<nkeep; i++) {
+         TParticle * parton = (TParticle*)array->AddrAt(i);
+         if (parton->GetStatusCode()==2 && !includeStatusCode2Particles) {
+            continue;
+         }
+         TLorentzVector v;
+         parton->Momentum(v);
+         Double_t m1 = v.M(); // invariant mass from the particle in the HepMC event
+         TParticlePDG *dbParton = parton->GetPDG();
+         if (!dbParton) {
+            if (useStdErr) cerr << "Warning: could not look up particle with PDG Code: " << parton->GetPdgCode() << endl << endl;
+            continue;
+         }
+         Double_t m2 = dbParton->Mass();
+         bool checkok;
+         if (m2 == 0) {
+            checkok = abs(m1) < 0.0001; // no such thing as negative mass...
+         } else {
+            checkok = abs(1 - m1/m2) < 0.01;
+         }
+         if (!checkok && useStdErr) {
+            cerr << "Problem: " << GetParticleName(parton) << " HepMC:" << m1 << endl;
+            cerr << ListReactionChain(array,i);
+            cerr << endl;
+         }
+         if (!checkok)
+            valid = false;
+      }
+   }
+   delete array;
+   array = 0;
+   return valid;
+}
+
+bool THepMCParser::IsValidVertexInvariantMass(bool useStdErr, bool includeStatusCode2Particles)
+{
+   bool valid = true;
+   TClonesArray * array = new TClonesArray("TParticle");
+   TBranch* branch = fTree->GetBranch("Particles");
+   branch->SetAddress(&array);
+   Int_t count = branch->GetEntries();
+   for (Int_t idx=0; idx<count; ++idx) {
+      array->Clear();
+      branch->GetEntry(idx); // "fill" the array
+      TLorentzVector v_st1;
+      TLorentzVector v_st4;
+      Int_t nkeep = array->GetEntriesFast();
+      for (Int_t i=0; i<nkeep; i++) {
+         TParticle * parton = (TParticle*)array->AddrAt(i);
+         TLorentzVector v_in;
+         parton->Momentum(v_in);
+         if (parton->GetStatusCode()==4) {
+            v_st4 += v_in;
+         } else if (parton->GetStatusCode()==1) {
+            v_st1 += v_in;
+         }
+         if (!includeStatusCode2Particles) { // only check beam particles vs final particles
+            continue;
+         }
+         Int_t fd = parton->GetFirstDaughter();
+         Int_t ld = parton->GetLastDaughter();
+         if (fd == -1) continue; // no daughters, continue loop
+         Int_t mother2 = -1;
+         TLorentzVector v_out;
+         bool oneok = false; bool allok = true; bool agreemother2 = true;
+         ostringstream daughterpdg;
+         ostringstream motherpdg;
+         for (Int_t id=TMath::Max(fd,0); id<=ld; id++) {
+            TParticle * daughter = (TParticle*)array->AddrAt(id);
+            if (fd==id) {
+               daughter->Momentum(v_out);
+               mother2 = daughter->GetSecondMother();
+            } else {
+               TLorentzVector d;
+               daughter->Momentum(d);
+               v_out += d;
+               if (daughter->GetSecondMother() != mother2) agreemother2 = false;
+            }
+            if (daughter->GetFirstMother() == i) {
+               oneok = true;
+            } else {
+               allok = false;
+            }
+            daughterpdg << " " << daughter->GetPdgCode();
+         }
+         motherpdg << " " << parton->GetPdgCode();
+         if (mother2 > -1 && agreemother2) {
+            TParticle * m2 = (TParticle*)array->AddrAt(mother2);
+            TLorentzVector m2v;
+            m2->Momentum(m2v);
+            v_in += m2v;
+            motherpdg << " " << m2->GetPdgCode();
+         }
+         if (oneok && allok && agreemother2) {
+            bool checkok = abs(1 - v_in.M()/v_out.M()) < 0.1;
+            if (!checkok) valid=false;
+            if (!checkok && useStdErr) {
+               //             cerr << "Problem: " << i << "[" << fd << "," << ld << "] PDG:" << motherpdg.str() << " ->" << daughterpdg.str() << "  Inv.mass: " << v_in.M() << " -> " << v_out.M() << endl;
+               cerr << ListReactionChain(array,i);
+               cerr << endl;
+               //             cerr << v_in.Px() << " " << v_in.Py() << " " << v_in.Pz() << " " << v_in.E() << endl;
+            } else if (useStdErr) {
+               //             cerr << "OK: " << i << "[" << fd << "," << ld << "] PDG:" << motherpdg.str() << " ->" << daughterpdg.str() << "  Inv.mass: " << v_in.M() << " -> " << v_out.M() << endl;
+            }
+         }
+      }
+      bool checkok = abs(1 - v_st4.M()/v_st1.M()) < 0.001;
+      if (!checkok) valid=false;
+      if (!checkok && useStdErr) {
+         cerr << " BEAM PARTICLES -> FINAL PARTICLES " << endl;
+         cerr << " status 4 (" << v_st4.M() << ") -> status 1 (" << v_st1.M() << ")" << endl << endl;
+      }
+   }
+   delete array;
+   array = 0;
+   return valid;
+}
+
+string THepMCParser::GetParticleName(TParticle * thisPart)
+{
+   TParticlePDG *dbPart = thisPart->GetPDG();
+   ostringstream name;
+   if (dbPart) {
+      name << dbPart->GetName() << "(" << dbPart->PdgCode() << "," << dbPart->Mass() << ")";
+   } else {
+      name << thisPart->GetPdgCode() << "(NoDBinfo)";
+   }
+   return name.str();
+}
+
+string THepMCParser::ListReactionChain(TClonesArray * particles, Int_t particleId)
+{
+   ostringstream output;
+
+   TParticle * part = (TParticle*)particles->AddrAt(particleId);
+   Int_t m1id = part->GetFirstMother();
+   Int_t m2id = part->GetSecondMother();
+   if (m1id > 1) { // ignore the initial collision with beam particles
+      ostringstream inStr;
+      ostringstream outStr;
+      TParticle * m1 = (TParticle*)particles->AddrAt(m1id);
+      TLorentzVector v_in;
+      m1->Momentum(v_in);
+      inStr << GetParticleName(m1) << "[" << v_in.M() << "]";
+      if (m2id > 1) {
+         TParticle * m2 = (TParticle*)particles->AddrAt(m2id);
+         TLorentzVector v_m2;
+         m2->Momentum(v_m2);
+         v_in += v_m2;
+         inStr << " " << GetParticleName(m2) << "[" << v_m2.M() << "]";
+      }
+      Int_t fd = m1->GetFirstDaughter();
+      Int_t ld = m1->GetLastDaughter();
+      TLorentzVector v_out;
+      part->Momentum(v_out);
+      outStr << GetParticleName(part) << "[" << v_out.M() << "]";
+      for (Int_t i=fd; i<=ld; ++i) {
+         if (i!=particleId) {
+            TParticle * d = (TParticle*)particles->AddrAt(i);
+            TLorentzVector v_d;
+            d->Momentum(v_d);
+            v_out += v_d;
+            outStr << " " << GetParticleName(d) << "[" << v_d.M() << "]";
+         }
+      }
+      output << "Parent reaction, inv mass: " << v_in.M() << " -> " << v_out.M() << endl;
+      output << " - partons: " << inStr.str() << " -> " << outStr.str() << endl;
+   }
+   Int_t fd = part->GetFirstDaughter();
+   Int_t ld = part->GetLastDaughter();
+   if (fd > -1) {
+      ostringstream inStr;
+      ostringstream outStr;
+      TLorentzVector v_in;
+      part->Momentum(v_in);
+      inStr << GetParticleName(part) << "[" << v_in.M() << "]";
+
+      TParticle * f = (TParticle*)particles->AddrAt(fd);
+      m2id = f->GetSecondMother();
+      if (m2id == particleId) {
+         m2id = f->GetFirstMother();
+      }
+      if (m2id > -1) {
+         TParticle * m2 = (TParticle*)particles->AddrAt(m2id);
+         TLorentzVector v_m2;
+         m2->Momentum(v_m2);
+         v_in += v_m2;
+         inStr << " " << GetParticleName(m2) << "[" << v_m2.M() << "]";
+      }
+      TLorentzVector v_out;
+      f->Momentum(v_out);
+      outStr << GetParticleName(f) << "[" << v_out.M() << "]";
+      for (Int_t i=fd+1; i<=ld; ++i) {
+         TParticle * d = (TParticle*)particles->AddrAt(i);
+         TLorentzVector v_d;
+         d->Momentum(v_d);
+         v_out += v_d;
+         outStr << " " << GetParticleName(d) << "[" << v_d.M() << "]";
+      }
+      output << "Child reaction, inv mass: " << v_in.M() << " -> " << v_out.M() << endl;
+      output << " - partons: " << inStr.str() << " -> " << outStr.str() << endl;
+   } else {
+      output << "Child reaction" << endl << " - none" << endl;
+   }
+
+   return output.str();
+}
+
+
+string THepMCParser::ParseGenEvent2TCloneArray(HepMC::GenEvent * genEvent, TClonesArray * array, bool requireSecondMotherBeforeDaughters)
+{
+   if (requireSecondMotherBeforeDaughters) {
+      return "requireSecondMotherBeforeDaughters not implemented yet!";
+   }
+   array->Clear();
+   ostringstream errMsgStream;
+   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
+
+   // Check event with HepMC's internal validation algorithm
+   if (!genEvent->is_valid()) {
+      errMsgStream << "Error with event id: " << genEvent->event_number() << ", event is not valid!";
+      return errMsgStream.str();
+   }
+
+   // Pull out the beam particles from the event
+   const pair<HepMC::GenParticle *,HepMC::GenParticle *> beamparts = genEvent->beam_particles();
+
+   // Four sanity checks:
+   // - Beam particles exists and are not the same
+   // - Both beam particles should have no production vertices, they come from the beams
+   // - Both beam particles should have defined end vertices, as they both should contribute
+   // - Both beam particles should have the exact same end vertex
+   if (!beamparts.first || !beamparts.second || beamparts.first->barcode()==beamparts.second->barcode()) {
+      errMsgStream << "Error with event id: " << genEvent->event_number() << ", beam particles doesn't exists or are the same";
+      return errMsgStream.str();
+   }
+   if (beamparts.first->production_vertex() || beamparts.second->production_vertex()) {
+      errMsgStream << "Error with event id: " << genEvent->event_number() << ", beam particles have production vertex/vertices...";
+      return errMsgStream.str();
+   }
+   if (!beamparts.first->end_vertex() || !beamparts.second->end_vertex()) {
+      errMsgStream << "Error with event id: " << genEvent->event_number() << ", beam particles have undefined end vertex/vertices...";
+      return errMsgStream.str();
+   }
+   if (beamparts.first->end_vertex() != beamparts.second->end_vertex()) {
+      errMsgStream << "Error with event id: " << genEvent->event_number() << ", beam particles do not collide in the same end vertex.";
+      return errMsgStream.str();
+   }
+
+   // Set the array to hold the number of particles in the event
+   array->Expand(genEvent->particles_size());
+
+   // Create a TParticle for each beam particle
+   new((*array)[0]) TParticle(
+         beamparts.first->pdg_id(),
+         beamparts.first->status(), // check if status has the same meaning
+         -1, // no mother1
+         -1, // no mother2
+         -1, // first daughter not known yet
+         -1, // last daughter not known yet
+         beamparts.first->momentum().px(),
+         beamparts.first->momentum().py(),
+         beamparts.first->momentum().pz(),
+         beamparts.first->momentum().e(),
+         0, // no production vertex, so zero?
+         0,
+         0,
+         0
+   );
+   partonMemory[beamparts.first->barcode()] = 0;
+   new((*array)[1]) TParticle(
+         beamparts.second->pdg_id(),
+         beamparts.second->status(),
+         -1, // no mother1
+         -1, // no mother2
+         -1, // first daughter not known yet
+         -1, // last daughter not known yet
+         beamparts.second->momentum().px(),
+         beamparts.second->momentum().py(),
+         beamparts.second->momentum().pz(),
+         beamparts.second->momentum().e(),
+         0, // no production vertex, so zero?
+         0,
+         0,
+         0
+   );
+   partonMemory[beamparts.second->barcode()] = 1;
+   Int_t arrayID = 2; // start counting IDs after the beam particles
+
+   // Do first vertex as a special case for performance, since its often has the most daughters and both mothers are known
+   Int_t firstDaughter = arrayID;
+   for (HepMC::GenVertex::particles_out_const_iterator iter = beamparts.first->end_vertex()->particles_out_const_begin();
+         iter != beamparts.first->end_vertex()->particles_out_const_end();
+         ++iter) {
+      new((*array)[arrayID]) TParticle(
+            (*iter)->pdg_id(),
+            (*iter)->status(),
+            0, // beam particle 1
+            1, // beam particle 2
+            -1, // first daughter not known yet
+            -1, // last daughter not known yet
+            (*iter)->momentum().px(),
+            (*iter)->momentum().py(),
+            (*iter)->momentum().pz(),
+            (*iter)->momentum().e(),
+            beamparts.first->end_vertex()->position().x(),
+            beamparts.first->end_vertex()->position().y(),
+            beamparts.first->end_vertex()->position().z(),
+            beamparts.first->end_vertex()->position().t()
+      );
+      partonMemory[(*iter)->barcode()] = arrayID;
+      ++arrayID;
+   }
+   Int_t lastDaughter = arrayID-1;
+   ((TParticle*)array->AddrAt(0))->SetFirstDaughter(firstDaughter); // beam particle 1
+   ((TParticle*)array->AddrAt(0))->SetLastDaughter(lastDaughter);
+   ((TParticle*)array->AddrAt(1))->SetFirstDaughter(firstDaughter); // beam particle 2
+   ((TParticle*)array->AddrAt(1))->SetLastDaughter(lastDaughter);
+
+   // Build vertex list by exploring tree and sorting such that daughters comes after mothers
+   list<HepMC::GenVertex*> vertexList;
+   set<int> vertexSearchSet;
+   ExploreVertex(beamparts.first->end_vertex(),vertexList,vertexSearchSet,requireSecondMotherBeforeDaughters);
+
+   // Analyze each vertex
+   for (list<HepMC::GenVertex*>::iterator i = vertexList.begin(); i != vertexList.end(); ++i) {
+      HepMC::GenVertex * vertex = (*i);
+
+      // first establish mother relations
+      HepMC::GenVertex::particles_in_const_iterator iter = vertex->particles_in_const_begin();
+      if (iter == vertex->particles_in_const_end()) {
+         return "Particle without a mother, and its not a beam particle!";
+      }
+      int motherA = partonMemory[(*iter)->barcode()];
+      if (((TParticle*)array->AddrAt(motherA))->GetFirstDaughter() > -1) {
+         return "Trying to assign new daughters to a particle that already has daughters defined!";
+      }
+      ++iter;
+      int motherB = -1;
+      if (iter != vertex->particles_in_const_end()) {
+         motherB = partonMemory[(*iter)->barcode()];
+         if (((TParticle*)array->AddrAt(motherB))->GetFirstDaughter() > -1) {
+            return "Trying to assign new daughters to a particle that already has daughters defined!";
+         }
+         ++iter;
+         if (iter != vertex->particles_in_const_end()) {
+            return "Particle with more than two mothers!";
+         }
+      }
+      if (motherB > -1 && motherB < motherA) {
+         int swap = motherA; motherA = motherB; motherB = swap;
+      }
+
+      // add the particles to the array, important that they are add in succession with respect to arrayID
+      firstDaughter = arrayID;
+      for (iter = vertex->particles_out_const_begin();
+            iter != vertex->particles_out_const_end();
+            ++iter) {
+         new((*array)[arrayID]) TParticle(
+               (*iter)->pdg_id(),
+               (*iter)->status(),
+               motherA, // mother 1
+               motherB, // mother 2
+               -1, // first daughter, if applicable, not known yet
+               -1, // last daughter, if applicable, not known yet
+               (*iter)->momentum().px(),
+               (*iter)->momentum().py(),
+               (*iter)->momentum().pz(),
+               (*iter)->momentum().e(),
+               vertex->position().x(),
+               vertex->position().y(),
+               vertex->position().z(),
+               vertex->position().t()
+         );
+         partonMemory[(*iter)->barcode()] = arrayID;
+         ++arrayID;
+      }
+      lastDaughter = arrayID-1;
+      if (lastDaughter < firstDaughter) {
+         return "Vertex with no out particles, should not be possible!";
+      }
+      // update mother with daughter interval
+      ((TParticle*)array->AddrAt(motherA))->SetFirstDaughter(firstDaughter);
+      ((TParticle*)array->AddrAt(motherA))->SetLastDaughter(lastDaughter);
+      if (motherB > -1) {
+         ((TParticle*)array->AddrAt(motherB))->SetFirstDaughter(firstDaughter);
+         ((TParticle*)array->AddrAt(motherB))->SetLastDaughter(lastDaughter);
+      }
+   }
+
+   return "";
+}
+
+
+void THepMCParser::ExploreVertex(HepMC::GenVertex * vertex, list<HepMC::GenVertex*> & vertexList, set<int> & vertexSearchSet, bool requireSecondMotherBeforeDaughters)
+{
+   for (HepMC::GenVertex::particles_out_const_iterator iter = vertex->particles_out_const_begin();
+         iter != vertex->particles_out_const_end();
+         ++iter) {
+      HepMC::GenVertex * testVertex = (*iter)->end_vertex();
+      if (testVertex) {
+         bool foundVertex = vertexSearchSet.find((*iter)->end_vertex()->barcode()) != vertexSearchSet.end();
+         if (requireSecondMotherBeforeDaughters) {
+            // redo this algorithem to move subtree instead of node....
+            // its not completely error proof in its current implementation even though the error is extremely rare
+
+            if (foundVertex) for (list<HepMC::GenVertex*>::iterator i = vertexList.begin(); i != vertexList.end(); ++i) {
+               if ((*i)->barcode() == testVertex->barcode()) {
+                  vertexList.erase(i);
+                  cout << " it happened, the vertex parsing order had to be changed " << endl;
+                  break;
+               }
+            } else {
+               vertexSearchSet.insert((*iter)->end_vertex()->barcode());
+            }
+            vertexList.push_back(testVertex);
+            if (!foundVertex) ExploreVertex(testVertex,vertexList,vertexSearchSet,requireSecondMotherBeforeDaughters);
+
+         } else {
+            if (!foundVertex) {
+               vertexSearchSet.insert((*iter)->end_vertex()->barcode());
+               vertexList.push_back(testVertex);
+               ExploreVertex(testVertex,vertexList,vertexSearchSet,requireSecondMotherBeforeDaughters);
+            }
+         }
+      }
+   }
+}
+
+
+
+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";
+const char * THepMCParser::fgPdfHeaderBranchString = "id1/I:id2:pdf_id1:pdf_id2:x1/D:x2:scalePDF:pdf1:pdf2";
+
+string THepMCParser::ParseGenEvent2HeaderStructs(HepMC::GenEvent * genEvent, HeavyIonHeader_t & heavyIonHeader, PdfHeader_t & pdfHeader, bool fillZeroOnMissingHeavyIon, bool fillZeroOnMissingPdf)
+{
+   HepMC::HeavyIon * heavyIon = genEvent->heavy_ion();
+   HepMC::PdfInfo * pdfInfo = genEvent->pdf_info();
+   if ((!heavyIon && !fillZeroOnMissingHeavyIon) || (!pdfInfo && !fillZeroOnMissingPdf)) {
+      return "HeavyIonInfo and/or PdfInfo not defined for this event, did you read it with IO_GenEvent?";
+   }
+   if (heavyIon) {
+      heavyIonHeader.Ncoll_hard = heavyIon->Ncoll_hard();
+      heavyIonHeader.Npart_proj = heavyIon->Npart_proj();
+      heavyIonHeader.Npart_targ = heavyIon->Npart_targ();
+      heavyIonHeader.Ncoll = heavyIon->Ncoll();
+      heavyIonHeader.spectator_neutrons = heavyIon->spectator_neutrons();
+      heavyIonHeader.spectator_protons = heavyIon->spectator_protons();
+      heavyIonHeader.N_Nwounded_collisions = heavyIon->N_Nwounded_collisions();
+      heavyIonHeader.Nwounded_N_collisions = heavyIon->Nwounded_N_collisions();
+      heavyIonHeader.Nwounded_Nwounded_collisions = heavyIon->Nwounded_Nwounded_collisions();
+      heavyIonHeader.impact_parameter = heavyIon->impact_parameter();
+      heavyIonHeader.event_plane_angle = heavyIon->event_plane_angle();
+      heavyIonHeader.eccentricity = heavyIon->eccentricity();
+      heavyIonHeader.sigma_inel_NN = heavyIon->sigma_inel_NN();
+   } else {
+      heavyIonHeader.Ncoll_hard = 0;
+      heavyIonHeader.Npart_proj = 0;
+      heavyIonHeader.Npart_targ = 0;
+      heavyIonHeader.Ncoll = 0;
+      heavyIonHeader.spectator_neutrons = 0;
+      heavyIonHeader.spectator_protons = 0;
+      heavyIonHeader.N_Nwounded_collisions = 0;
+      heavyIonHeader.Nwounded_N_collisions = 0;
+      heavyIonHeader.Nwounded_Nwounded_collisions = 0;
+      heavyIonHeader.impact_parameter = 0.0;
+      heavyIonHeader.event_plane_angle = 0.0;
+      heavyIonHeader.eccentricity = 0.0;
+      heavyIonHeader.sigma_inel_NN = 0.0;
+   }
+   if (pdfInfo) {
+      pdfHeader.id1 = pdfInfo->id1();
+      pdfHeader.id2 = pdfInfo->id2();
+      pdfHeader.pdf_id1 = pdfInfo->pdf_id1();
+      pdfHeader.pdf_id2 = pdfInfo->pdf_id2();
+      pdfHeader.x1 = pdfInfo->x1();
+      pdfHeader.x2 = pdfInfo->x2();
+      pdfHeader.scalePDF = pdfInfo->scalePDF();
+      pdfHeader.pdf1 = pdfInfo->pdf1();
+      pdfHeader.pdf2 = pdfInfo->pdf2();
+   } else {
+      pdfHeader.id1 = 0;
+      pdfHeader.id2 = 0;
+      pdfHeader.pdf_id1 = 0;
+      pdfHeader.pdf_id2 = 0;
+      pdfHeader.x1 = 0.0;
+      pdfHeader.x2 = 0.0;
+      pdfHeader.scalePDF = 0.0;
+      pdfHeader.pdf1 = 0.0;
+      pdfHeader.pdf2 = 0.0;
+   }
+   return "";
+}
+
+
+
+
+
+
+
diff --git a/TEvtGen/THepMCParser.h b/TEvtGen/THepMCParser.h
new file mode 100644 (file)
index 0000000..03d73c5
--- /dev/null
@@ -0,0 +1,115 @@
+// module identifier line...
+// Author: Brian Thorsbro, 24/6-2014
+
+
+#ifndef THEPMCPARSER_H
+#define THEPMCPARSER_H
+
+#include <string>
+#include <list>
+#include <set>
+#include "HepMC/IO_GenEvent.h"
+#include "TTree.h"
+#include "TClonesArray.h"
+#include "TParticle.h"
+
+
+class THepMCParser {
+
+private:
+
+   TTree * fTree;
+   void init(HepMC::IO_BaseClass *);
+
+   // 'name(pdgcode,invmass)' if db entry or 'pdgcode()' if no db entry
+   static std::string GetParticleName(TParticle *);
+
+   // Fold out the vertex tree in a list, such that daughters are always last
+   static void ExploreVertex(HepMC::GenVertex *, std::list<HepMC::GenVertex*> &, std::set<int> &, bool);
+
+
+public:
+
+   // Header struct
+   static const char * 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";
+   struct HeavyIonHeader_t {
+      Int_t   Ncoll_hard;                    // Number of hard scatterings
+      Int_t   Npart_proj;                    // Number of projectile participants
+      Int_t   Npart_targ;                    // Number of target participants
+      Int_t   Ncoll;                         // Number of NN (nucleon-nucleon) collisions
+      Int_t   spectator_neutrons;            // Number of spectator neutrons
+      Int_t   spectator_protons;             // Number of spectator protons
+      Int_t   N_Nwounded_collisions;         // Number of N-Nwounded collisions
+      Int_t   Nwounded_N_collisions;         // Number of Nwounded-N collisons
+      Int_t   Nwounded_Nwounded_collisions;  // Number of Nwounded-Nwounded collisions
+      Float_t impact_parameter;              // Impact Parameter(in fm) of collision
+      Float_t event_plane_angle;             // Azimuthal angle of event plane
+      Float_t eccentricity;                  // eccentricity of participating nucleons in the transverse plane (as in phobos nucl-ex/0510031)
+      Float_t sigma_inel_NN;                 // nucleon-nucleon inelastic (including diffractive) cross-section
+   };
+   static const char * fgPdfHeaderBranchString; // "id1/I,id2,pdf_id1,pdf_id2,x1/D,x2,scalePDF,pdf1,pdf2";
+   struct PdfHeader_t {
+      Int_t    id1;        // flavour code of first parton
+      Int_t    id2;        // flavour code of second parton
+      Int_t    pdf_id1;    // LHAPDF set id of first parton
+      Int_t    pdf_id2;    // LHAPDF set id of second parton
+      Double_t x1;         // fraction of beam momentum carried by first parton ("beam side")
+      Double_t x2;         // fraction of beam momentum carried by second parton ("target side")
+      Double_t scalePDF;   // Q-scale used in evaluation of PDF's   (in GeV)
+      Double_t pdf1;       // PDF (id1, x1, Q) - x*f(x)
+      Double_t pdf2;       // PDF (id2, x2, Q) - x*f(x)
+   };
+
+   // Default constructor/destructor stuff, don't inherit from this class unless you handle the tree pointer
+   inline THepMCParser() : fTree(0) {;} // nullptr in c++11
+   inline virtual ~THepMCParser() {;} // should be a memory management of the TTree...
+
+   // The actual useful constructors, either take:
+   //  - a file name for a file with HepMC data or
+   //  - a HepMC event data structure
+   THepMCParser(const char *);
+   THepMCParser(HepMC::IO_BaseClass *);
+
+   // Optional validators, set the argument to true for verbose output to STDERR
+   // WARNING: including status code 2 may produce invalid flag when it is in fact valid, not recommended to use
+   bool IsValidMotherDaughtersConsitency(bool useStdErr = false, bool requireSecondMotherBeforeDaughters = false);
+   bool IsValidParticleInvariantMass(bool useStdErr = false, bool includeStatusCode2Particles = false);
+   bool IsValidVertexInvariantMass(bool useStdErr = false, bool includeStatusCode2Particles = false);
+
+   // Show the decay chain by point of view of some particle
+   static std::string ListReactionChain(TClonesArray *, Int_t);
+
+   // Access the TTree generated or write it to a file
+   TTree * GetTTree();
+   void WriteTTreeToFile(const char *);
+
+   // ***** The constructors will parse all the events with this function for you *****
+   // The work horse of the parser, takes one event and generates the corresponding TClonesArray of TParticles,
+   // requireSecondMotherBeforeDaughters defaults to false
+   //  - Success if length of return string is 0, the provided TClonesArray has been filled
+   //  - Failure otherwise and the return string then contains the error message
+   //
+   // Note:
+   // The TClonesArray must be initialized before hand
+   // the array will be cleared and filled with TParticles
+   // The capacity of the array will be set to the number of particles parsed
+   // First mother will be before the daughters, but second mother may be after
+   // The two first particles in the array will be the beam particles
+   // The status code set on the particle is copied from HepMC, i.e.
+   //  - 1 = final particle
+   //  - 2 = transitory particle
+   //  - 4 = beam particle
+   // The function is static to enable customized wrappers around it
+   static std::string ParseGenEvent2TCloneArray(HepMC::GenEvent *, TClonesArray *, bool requireSecondMotherBeforeDaughters = false);
+
+
+   // Depending on the implementation of HepMC::IO_BaseClass there may be information on
+   // heavy ions or parton distribution functions available. This function will pull them
+   // out if available.
+   // The caller must supply allocated structures which will then be filled.
+   static std::string ParseGenEvent2HeaderStructs(HepMC::GenEvent *, HeavyIonHeader_t &, PdfHeader_t &, bool fillZeroOnMissingHeavyIon = true, bool fillZeroOnMissingPdf = true);
+
+
+};
+
+#endif
diff --git a/TEvtGen/THepMCParserLinkDef.h b/TEvtGen/THepMCParserLinkDef.h
new file mode 100644 (file)
index 0000000..c140fea
--- /dev/null
@@ -0,0 +1,9 @@
+#ifdef __CINT__
+
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+
+#pragma link C++ class THepMCParser-;
+
+#endif