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