Overlaps corrected, new shape of sectors
[u/mrichter/AliRoot.git] / 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 }