Line data    Source code

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