Changes for alien analysis.
[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 <AliRunLoader.h>
18 #include <AliStack.h>
19
20 #include "AliJetParticle.h"
21 #include "AliJetEventParticles.h"
22 #include "AliJetParticlesReaderKine.h"
23 #include "AliHeader.h"
24 #include "AliGenPythiaEventHeader.h"
25
26 ClassImp(AliJetParticlesReaderKine)
27
28 AliJetParticlesReaderKine::AliJetParticlesReaderKine() :
29   AliJetParticlesReader(),
30   fFileName("galice.root"),
31   fRunLoader(0),
32   fNeutral(kTRUE),fCharged(kTRUE),fEM(kTRUE),
33   fUseTracks(kFALSE)
34 {
35   //constructor
36 }
37
38 AliJetParticlesReaderKine::AliJetParticlesReaderKine(TString& fname) :
39   AliJetParticlesReader(),
40   fFileName(fname),
41   fRunLoader(0),
42   fNeutral(kTRUE),fCharged(kTRUE),fEM(kTRUE),
43   fUseTracks(kFALSE)
44 {
45   //constructor
46 }
47
48 AliJetParticlesReaderKine::AliJetParticlesReaderKine(TObjArray* dirs, const Char_t *filename):
49   AliJetParticlesReader(dirs),
50   fFileName(filename),
51   fRunLoader(0),
52   fNeutral(kTRUE),fCharged(kTRUE),fEM(kTRUE),
53   fUseTracks(kFALSE)
54 {
55   //constructor
56 }
57
58 AliJetParticlesReaderKine::~AliJetParticlesReaderKine()
59 {
60   //destructor
61   if(fRunLoader) delete fRunLoader;
62 }
63
64 void AliJetParticlesReaderKine::Rewind()
65 {
66   //Rewinds to the beginning
67   if(fRunLoader) delete fRunLoader;
68   fRunLoader   = 0;
69   fCurrentDir  = 0;
70   fNEventsRead = 0;  
71 }
72
73 Int_t AliJetParticlesReaderKine::ReadNext()
74 {
75   //Reads Kinematics Tree
76
77   if((!fOwner) || (fEventParticles == 0)) 
78     fEventParticles = new AliJetEventParticles();
79
80   while(fCurrentDir < GetNumberOfDirs())
81     { 
82
83       if (!OpenFile(fCurrentDir)) 
84       { 
85         delete fRunLoader; //close current session
86         fRunLoader = 0;    //assure pointer is null
87         fCurrentDir++;
88         continue;
89       }
90     
91       if (fCurrentEvent == fRunLoader->GetNumberOfEvents())
92         {
93           //read next directory
94           delete fRunLoader; //close current session
95           fRunLoader = 0;    //assure pointer is null
96           fCurrentDir++;     //go to next dir
97           continue; 
98         }
99      
100       Info("ReadNext","Reading Event %d",fCurrentEvent);
101       fRunLoader->GetEvent(fCurrentEvent);
102       AliStack* stack = fRunLoader->Stack();
103       if (!stack)
104         {
105           Error("ReadNext","Can't get stack for event %d",fCurrentEvent);
106           continue;
107         }
108
109       //clear old event
110       Int_t nprim = stack->GetNprimary();
111       Int_t npart = nprim;
112       if(fUseTracks)
113         npart = stack->GetNtrack();
114       fEventParticles->Reset(npart);
115
116       TString headdesc="";
117       AliHeader *header=fRunLoader->GetHeader();
118       if(!header) {
119         Warning("ReadNext","Header not found in event %d",fCurrentEvent);
120       } else {
121         headdesc+="Run ";
122         headdesc+=header->GetRun();
123         headdesc+=": Ev ";
124         headdesc+=header->GetEventNrInRun();
125         AliGenPythiaEventHeader *pheader=(AliGenPythiaEventHeader*)header->GenEventHeader();
126         if(!pheader) {
127           Warning("ReadNext","Pythia-Header not found in event %d",fCurrentEvent);
128         } else {
129           Int_t ntruq=pheader->NUQTriggerJets();
130           if(ntruq){
131             Double_t x0=pheader->GetXJet();
132             Double_t y0=pheader->GetYJet();
133             Double_t zquench[4];
134             pheader->GetZQuench(zquench);
135             fEventParticles->SetXYJet(x0,y0);
136             fEventParticles->SetZQuench(zquench);
137             for(Int_t j=0;j<ntruq;j++){
138               Float_t pjet[4];
139               pheader->UQJet(j,pjet);
140               fEventParticles->AddUQJet(pjet);
141             }
142           }
143           //Int_t ptyp=pheader->ProcessType();
144           Int_t ntrials=pheader->Trials();
145           headdesc+=": Tr ";
146           headdesc+=ntrials;
147           Int_t ntr=pheader->NTriggerJets();
148           if(ntr){
149             for(Int_t j=0;j<ntr;j++){
150               Float_t pjet[4];
151               pheader->TriggerJet(j,pjet);
152               fEventParticles->AddJet(pjet);
153               if(!ntruq) fEventParticles->AddUQJet(pjet);
154             }
155           }
156         }
157       }
158       fEventParticles->SetHeader(headdesc);
159
160       //get vertex
161       const TParticle *kv = stack->Particle(0);
162       if(kv) {
163         fEventParticles->SetVertex(kv->Vx(),kv->Vy(),kv->Vz());
164       }
165
166       //loop over particles
167       for (Int_t i = 0;i<npart; i++)
168         {
169           TParticle *p = stack->Particle(i);
170           if(!p) continue;
171           Int_t child1 = p->GetFirstDaughter();
172           //Int_t child2 = p->GetLastDaughter();        
173           //Int_t mother = p->GetFirstMother();    
174           //cout << child1 << " " << child2 << " " << mother << endl;
175           if((child1>=0) && (child1<nprim)) continue; 
176           if(p->GetStatusCode()!=1){
177             //p->Print();
178             continue;
179           }
180           if(IsAcceptedParticle(p)) //put particle in event
181             fEventParticles->AddParticle(p,i); 
182         }
183       fCurrentEvent++;
184       fNEventsRead++;
185       return kTRUE;
186     }
187
188   //end of loop over directories specified in fDirs ObjArray
189   return kFALSE;
190 }
191
192 Int_t AliJetParticlesReaderKine::OpenFile(Int_t n)
193 {
194   //opens file with kine tree
195
196   if(fRunLoader){
197     if(fCurrentEvent < fRunLoader->GetNumberOfEvents()) return kTRUE;
198     else return kFALSE;
199   }
200
201   const TString& dirname = GetDirName(n);
202   if (dirname == "")
203     { 
204       Error("OpenNextFile","Can't get directory name with index %d",n);
205       return kFALSE;
206     }
207
208   TString filename = dirname +"/"+ fFileName;
209   fRunLoader = AliRunLoader::Open(filename.Data()); 
210
211   if (fRunLoader == 0)
212     {
213       Error("OpenNextFile","Can't open session from file %s",filename.Data());
214       return kFALSE;
215     }
216   
217   if (fRunLoader->GetNumberOfEvents() <= 0)
218     {
219       Error("OpenNextFile","There is no events in this directory.");
220       delete fRunLoader;
221       fRunLoader = 0;
222       return kFALSE;
223     }
224
225   if (fRunLoader->LoadKinematics())
226     {
227       Error("OpenNextFile","Error occured while loading kinematics.");
228       return kFALSE;
229     }
230
231   fCurrentEvent = 0;
232   return kTRUE;
233 }