Little changes to run for pdc04
[u/mrichter/AliRoot.git] / JETAN / AliJetParticlesReaderKine.cxx
1 // $Id$
2
3 //_______________________________________________________________________
4 /////////////////////////////////////////////////////////////////////////
5 //
6 // class AliJetParticlesReaderKine
7 //
8 // Reader for Kinematics
9 //
10 // loizides@ikf.uni-frankfurt.de
11 //
12 /////////////////////////////////////////////////////////////////////////
13
14 #include <Riostream.h>
15 #include <TString.h>
16 #include <TParticle.h>
17 #include <TLorentzVector.h>
18 #include <AliRunLoader.h>
19 #include <AliStack.h>
20 #include <AliHeader.h>
21 #include <AliGenEventHeader.h>
22 #include <AliGenPythiaEventHeader.h>
23 #include <AliGenHijingEventHeader.h>
24 #include "AliJetParticle.h"
25 #include "AliJetEventParticles.h"
26 #include "AliJetParticlesReaderKine.h"
27
28 ClassImp(AliJetParticlesReaderKine)
29
30 AliJetParticlesReaderKine::AliJetParticlesReaderKine() :
31   AliJetParticlesReader(),
32   fFileName("galice.root"),
33   fRunLoader(0),
34   fNeutral(kTRUE),fCharged(kTRUE),fEM(kTRUE),
35   fUseTracks(kFALSE)
36 {
37   //constructor
38 }
39
40 AliJetParticlesReaderKine::AliJetParticlesReaderKine(TString& fname) :
41   AliJetParticlesReader(),
42   fFileName(fname),
43   fRunLoader(0),
44   fNeutral(kTRUE),fCharged(kTRUE),fEM(kTRUE),
45   fUseTracks(kFALSE)
46 {
47   //constructor
48 }
49
50 AliJetParticlesReaderKine::AliJetParticlesReaderKine(TObjArray* dirs, const Char_t *filename):
51   AliJetParticlesReader(dirs),
52   fFileName(filename),
53   fRunLoader(0),
54   fNeutral(kTRUE),fCharged(kTRUE),fEM(kTRUE),
55   fUseTracks(kFALSE)
56 {
57   //constructor
58 }
59
60 AliJetParticlesReaderKine::~AliJetParticlesReaderKine()
61 {
62   //destructor
63   if(fRunLoader) delete fRunLoader;
64 }
65
66 void AliJetParticlesReaderKine::Rewind()
67 {
68   //Rewinds to the beginning
69   if(fRunLoader) delete fRunLoader;
70   fRunLoader   = 0;
71   fCurrentDir  = 0;
72   fNEventsRead = 0;  
73 }
74
75 Int_t AliJetParticlesReaderKine::ReadNext()
76 {
77   //Reads Kinematics Tree
78
79   if((!fOwner) || (fEventParticles == 0)) 
80     fEventParticles = new AliJetEventParticles();
81
82   while(fCurrentDir < GetNumberOfDirs())
83     { 
84
85       if (!OpenFile(fCurrentDir)) 
86       { 
87         delete fRunLoader; //close current session
88         fRunLoader = 0;    //assure pointer is null
89         fCurrentDir++;
90         continue;
91       }
92     
93       if (fCurrentEvent == fRunLoader->GetNumberOfEvents())
94         {
95           //read next directory
96           delete fRunLoader; //close current session
97           fRunLoader = 0;    //assure pointer is null
98           fCurrentDir++;     //go to next dir
99           continue; 
100         }
101      
102       Info("ReadNext","Reading Event %d",fCurrentEvent);
103       fRunLoader->GetEvent(fCurrentEvent);
104       AliStack* stack = fRunLoader->Stack();
105       if (!stack)
106         {
107           Error("ReadNext","Can't get stack for event %d",fCurrentEvent);
108           continue;
109         }
110
111       //clear old event
112       Int_t nprim = stack->GetNprimary();
113       Int_t npart = nprim;
114       if(fUseTracks)
115         npart = stack->GetNtrack();
116       fEventParticles->Reset(npart);
117
118       TString headdesc="";
119       Int_t gentype=-1;
120       AliHeader *header=fRunLoader->GetHeader();
121       if(!header) {
122         Warning("ReadNext","Header not found in event %d",fCurrentEvent);
123       } else {
124         AliGenEventHeader *genheader=header->GenEventHeader();
125         TString hname=genheader->GetName();
126         if(hname.CompareTo("HIJINGparam")==0) {
127           gentype=1;      ;
128           headdesc+="HIJINGparam";
129         } else if(hname.CompareTo("Hijing")==0) {
130           gentype=2;
131           AliGenHijingEventHeader *hheader=(AliGenHijingEventHeader*)header->GenEventHeader();
132           headdesc+="Hijing";
133           if(!hheader) {
134             Warning("ReadNext","Hijing-Header not found in event %d",fCurrentEvent);
135           } else {
136             Int_t ntrials=hheader->Trials();
137             fEventParticles->SetTrials(ntrials);
138             fEventParticles->SetImpact(hheader->ImpactParameter());
139             fEventParticles->SetNhard(hheader->HardScatters());
140             fEventParticles->SetNpart(hheader->NwNw());
141
142             TLorentzVector  jet1; 
143             TLorentzVector  jet2; 
144             TLorentzVector  jetFsr1;
145             TLorentzVector  jetFsr2;
146             hheader->GetJets(jet1,jet2,jetFsr1,jetFsr2);
147             fEventParticles->AddHard(0,jet1.X(),jet1.Y(),jet1.Z(),jet1.T(),0);
148             fEventParticles->AddHard(1,jet2.X(),jet2.Y(),jet2.Z(),jet2.T(),0);
149           }
150         } else if(hname.CompareTo("Pythia")==0){
151           gentype=3;
152           headdesc+="Pythia";
153           AliGenPythiaEventHeader *pheader=(AliGenPythiaEventHeader*)header->GenEventHeader();
154           if(!pheader) {
155             Warning("ReadNext","Pythia-Header not found in event %d",fCurrentEvent);
156           } else {
157             Int_t ntruq=0;
158 #ifndef NOUQHEADERINFO
159             ntruq=pheader->NUQTriggerJets();
160             if(ntruq){
161               Double_t x0=pheader->GetXJet();
162               Double_t y0=pheader->GetYJet();
163               if(x0==y0==-1){
164                 x0=y0=0.;
165               }
166               Double_t zquench[4];
167               pheader->GetZQuench(zquench);
168               fEventParticles->SetXYJet(x0,y0);
169               fEventParticles->SetZQuench(zquench);
170               for(Int_t j=0;j<ntruq;j++){
171                 Float_t pjet[4];
172                 pheader->UQJet(j,pjet);
173                 fEventParticles->AddUQJet(pjet);
174               }
175             }
176 #endif
177             //Int_t ptyp=pheader->ProcessType();
178             Int_t ntrials=pheader->Trials();
179             fEventParticles->SetTrials(ntrials);
180
181             Int_t ntr=pheader->NTriggerJets();
182             if(ntr){
183               for(Int_t j=0;j<ntr;j++){
184                 Float_t pjet[4];
185                 pheader->TriggerJet(j,pjet);
186                 fEventParticles->AddJet(pjet);
187                 if(!ntruq) fEventParticles->AddUQJet(pjet);
188               }
189             }
190             for(Int_t i=6;i<=7;i++){
191               TParticle *MP = stack->Particle(i);
192               if(!MP) break;
193               Int_t type=0;
194               if(MP->GetPdgCode()==21) type=2;
195               else type=1;
196               fEventParticles->AddHard(i-6,MP->Px(),MP->Py(),MP->Pz(),MP->Energy(),type);
197             }
198           }
199         }
200         headdesc+=" Run ";
201         headdesc+=header->GetRun();
202         headdesc+=": Ev ";
203         headdesc+=header->GetEventNrInRun();
204       }
205       fEventParticles->SetHeader(headdesc);
206
207       //get vertex
208       const TParticle *kv = stack->Particle(0);
209       if(kv) {
210         fEventParticles->SetVertex(kv->Vx(),kv->Vy(),kv->Vz());
211       }
212
213       //loop over particles
214       for (Int_t i = 0;i<npart; i++)
215         {
216           TParticle *p = stack->Particle(i);
217           if(!p) continue;
218           Int_t child1 = p->GetFirstDaughter();
219           //Int_t child2 = p->GetLastDaughter();        
220           //Int_t mother = p->GetFirstMother();    
221           //cout << child1 << " " << child2 << " " << mother << endl;
222           if((child1>=0) && (child1<nprim)) continue; 
223
224           //check status code depending on gentype
225           if(gentype==1){ // Hijing Param
226             if(p->GetStatusCode()!=0) continue;
227           } else if(gentype==2){ //Hijing
228             //if( p->GetStatusCode()!=0){
229             //  cout <<p->GetStatusCode() << " ";p->Print(); 
230             //}
231           } else if(gentype==3){ //Pythia
232             if(p->GetStatusCode()!=1) continue;
233             p->SetWeight(-123); //mark particle
234           }
235
236           //kinematic cuts
237           if(IsAcceptedParticle(p)){ //put particle in event
238             //if(p->Pt()>20){cout <<p->GetStatusCode() << " ";p->Print(); }
239             fEventParticles->AddParticle(p,i); 
240           }
241         }
242       fCurrentEvent++;
243       fNEventsRead++;
244       return kTRUE;
245     }
246
247   //end of loop over directories specified in fDirs ObjArray
248   return kFALSE;
249 }
250
251 Int_t AliJetParticlesReaderKine::OpenFile(Int_t n)
252 {
253   //opens file with kine tree
254
255   if(fRunLoader){
256     if(fCurrentEvent < fRunLoader->GetNumberOfEvents()) return kTRUE;
257     else return kFALSE;
258   }
259
260   const TString& dirname = GetDirName(n);
261   if (dirname == "")
262     { 
263       Error("OpenNextFile","Can't get directory name with index %d",n);
264       return kFALSE;
265     }
266
267   TString filename = dirname +"/"+ fFileName;
268   fRunLoader = AliRunLoader::Open(filename.Data()); 
269
270   if (fRunLoader == 0)
271     {
272       Error("OpenNextFile","Can't open session from file %s",filename.Data());
273       return kFALSE;
274     }
275   
276   if (fRunLoader->GetNumberOfEvents() <= 0)
277     {
278       Error("OpenNextFile","There is no events in this directory.");
279       delete fRunLoader;
280       fRunLoader = 0;
281       return kFALSE;
282     }
283
284   if (fRunLoader->LoadKinematics())
285     {
286       Error("OpenNextFile","Error occured while loading kinematics.");
287       return kFALSE;
288     }
289
290   fCurrentEvent = 0;
291   return kTRUE;
292 }