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    // Read the next event
  57    if ((fGenEvent = fEventsHandle->read_next_event())) {
  58       THepMCParser::ParseGenEvent2TCloneArray(fGenEvent,fParticleArray,false);
  59       fParticleIterator->Reset();
  60       THepMCParser::HeavyIonHeader_t heavyIonHeader;
  61       THepMCParser::PdfHeader_t pdfHeader;
  62       THepMCParser::ParseGenEvent2HeaderStructs(fGenEvent,heavyIonHeader,pdfHeader,true,true);
  63       fGenEventHeader = new AliGenHepMCEventHeader(
  64             heavyIonHeader.Ncoll_hard,
  65             heavyIonHeader.Npart_proj,
  66             heavyIonHeader.Npart_targ,
  67             heavyIonHeader.Ncoll,
  68             heavyIonHeader.spectator_neutrons,
  69             heavyIonHeader.spectator_protons,
  70             heavyIonHeader.N_Nwounded_collisions,
  71             heavyIonHeader.Nwounded_N_collisions,
  72             heavyIonHeader.Nwounded_Nwounded_collisions,
  73             heavyIonHeader.impact_parameter,
  74             heavyIonHeader.event_plane_angle,
  75             heavyIonHeader.eccentricity,
  76             heavyIonHeader.sigma_inel_NN,
  77             pdfHeader.id1,
  78             pdfHeader.id2,
  79             pdfHeader.pdf_id1,
  80             pdfHeader.pdf_id2,
  81             pdfHeader.x1,
  82             pdfHeader.x2,
  83             pdfHeader.scalePDF,
  84             pdfHeader.pdf1,
  85             pdfHeader.pdf2
  86       );
  87       printf("Parsed event with %d particles.\n", fGenEvent->particles_size());
  88       return fGenEvent->particles_size();
  89    }
  90    printf("No more events in the file.\n");
  91    return 0;
  92 }
  93
  94 TParticle* AliGenReaderHepMC::NextParticle()
  95 {
  96    // Read next particle
  97    TParticle * particle = (TParticle*)fParticleIterator->Next();
  98    if (particle && particle->GetStatusCode()==1) {
  99       particle->SetBit(kTransportBit);
 100    }
 101    return particle;
 102 }
 103
 104 void AliGenReaderHepMC::RewindEvent()
 105 {
 106    fParticleIterator->Reset();
 107 }