EVGEN/AliGenReaderHepMC.cxx

   1 #include <TVirtualMC.h>
   2 #include <TDatabasePDG.h>
   3 #include <TParticle.h>
   4
   5 #include "AliGenReaderHepMC.h"
   6 #include "AliRun.h"
   7 #include "AliStack.h"
   8 #include "AliGenHepMCEventHeader.h"
   9
  10 #include "HepMC/IO_BaseClass.h"
  11 #include "HepMC/GenEvent.h"
  12 #include "HepMC/IO_GenEvent.h"
  13
  14 ClassImp(AliGenReaderHepMC)
  15
  16 AliGenReaderHepMC::AliGenReaderHepMC():fEventsHandle(0), fGenEvent(0), fParticleArray(0), fParticleIterator(0), fGenEventHeader(0) {;}
  17
  18 AliGenReaderHepMC::AliGenReaderHepMC(const AliGenReaderHepMC &reader)
  19    :AliGenReader(reader), fEventsHandle(0), fGenEvent(0), fParticleArray(0), fParticleIterator(0), fGenEventHeader(0) {reader.Copy(*this);}
  20
  21
  22 AliGenReaderHepMC& AliGenReaderHepMC::operator=(const  AliGenReaderHepMC& rhs)
  23 {
  24    // Assignment operator
  25    rhs.Copy(*this);
  26    return *this;
  27 }
  28
  29 AliGenReaderHepMC::~AliGenReaderHepMC(){ delete fEventsHandle; delete fGenEvent; delete fParticleArray; delete fParticleIterator;} // not deleting fGenEventHeader as it is returned out
  30
  31 void AliGenReaderHepMC::Copy(TObject&) const
  32 {
  33    //
  34    // Copy
  35    //
  36    Fatal("Copy","Not implemented!\n");
  37 }
  38
  39 void AliGenReaderHepMC::Init()
  40 {
  41    // check if file exists, using FILE to avoid (the otherwise faster) POSIX dependencies
  42    if (FILE *file = fopen(fFileName,"r"))  {
  43       printf("File %s opened \n", fFileName);
  44       fclose(file);
  45    } else {
  46       printf("Couldn't open input file: %s \n", fFileName);
  47    }
  48    // Initialisation
  49    fEventsHandle = new HepMC::IO_GenEvent(fFileName, std::ios::in);
  50    fParticleArray = new TClonesArray("TParticle");
  51    fParticleIterator = new TIter(fParticleArray);
  52 }
  53
  54 Int_t AliGenReaderHepMC::NextEvent()
  55 {
  56    // Clean memory
  57    if (fGenEvent) delete fGenEvent;
  58    // Read the next event
  59    if ((fGenEvent = fEventsHandle->read_next_event())) {
  60       THepMCParser::ParseGenEvent2TCloneArray(fGenEvent,fParticleArray,"GEV","CM",false);
  61       fParticleIterator->Reset();
  62       THepMCParser::HeavyIonHeader_t heavyIonHeader;
  63       THepMCParser::PdfHeader_t pdfHeader;
  64       THepMCParser::ParseGenEvent2HeaderStructs(fGenEvent,heavyIonHeader,pdfHeader,true,true);
  65       fGenEventHeader = new AliGenHepMCEventHeader(
  66             heavyIonHeader.Ncoll_hard,
  67             heavyIonHeader.Npart_proj,
  68             heavyIonHeader.Npart_targ,
  69             heavyIonHeader.Ncoll,
  70             heavyIonHeader.spectator_neutrons,
  71             heavyIonHeader.spectator_protons,
  72             heavyIonHeader.N_Nwounded_collisions,
  73             heavyIonHeader.Nwounded_N_collisions,
  74             heavyIonHeader.Nwounded_Nwounded_collisions,
  75             heavyIonHeader.impact_parameter,
  76             heavyIonHeader.event_plane_angle,
  77             heavyIonHeader.eccentricity,
  78             heavyIonHeader.sigma_inel_NN,
  79             pdfHeader.id1,
  80             pdfHeader.id2,
  81             pdfHeader.pdf_id1,
  82             pdfHeader.pdf_id2,
  83             pdfHeader.x1,
  84             pdfHeader.x2,
  85             pdfHeader.scalePDF,
  86             pdfHeader.pdf1,
  87             pdfHeader.pdf2
  88       );
  89       printf("Parsed event %d with %d particles.\n", fGenEvent->event_number(), fGenEvent->particles_size());
  90       return fGenEvent->particles_size();
  91    }
  92    printf("No more events in the file.\n");
  93    return 0;
  94 }
  95
  96 TParticle* AliGenReaderHepMC::NextParticle()
  97 {
  98    // Read next particle
  99    TParticle * particle = (TParticle*)fParticleIterator->Next();
 100    if (particle && particle->GetStatusCode()==1) {
 101       particle->SetBit(kTransportBit);
 102    }
 103    return particle;
 104 }
 105
 106 void AliGenReaderHepMC::RewindEvent()
 107 {
 108    fParticleIterator->Reset();
 109 }