Bugfix.
[u/mrichter/AliRoot.git] / JETAN / AliJetParticlesReaderKine.cxx
CommitLineData
d7c6ab14 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>
d7073a83 17#include <TLorentzVector.h>
d7c6ab14 18#include <AliRunLoader.h>
19#include <AliStack.h>
d7073a83 20#include <AliHeader.h>
21#include <AliGenEventHeader.h>
22#include <AliGenPythiaEventHeader.h>
23#include <AliGenHijingEventHeader.h>
d7c6ab14 24#include "AliJetParticle.h"
25#include "AliJetEventParticles.h"
26#include "AliJetParticlesReaderKine.h"
27
28ClassImp(AliJetParticlesReaderKine)
29
d7c6ab14 30AliJetParticlesReaderKine::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
d7c6ab14 40AliJetParticlesReaderKine::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
d7c6ab14 50AliJetParticlesReaderKine::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
d7c6ab14 60AliJetParticlesReaderKine::~AliJetParticlesReaderKine()
61{
62 //destructor
63 if(fRunLoader) delete fRunLoader;
64}
65
d7c6ab14 66void AliJetParticlesReaderKine::Rewind()
67{
68 //Rewinds to the beginning
69 if(fRunLoader) delete fRunLoader;
04a02430 70 fRunLoader = 0;
71 fCurrentDir = 0;
72 fNEventsRead = 0;
d7c6ab14 73}
74
d7c6ab14 75Int_t AliJetParticlesReaderKine::ReadNext()
76{
77 //Reads Kinematics Tree
78
79 if((!fOwner) || (fEventParticles == 0))
80 fEventParticles = new AliJetEventParticles();
81
b2760c9e 82 while(fCurrentDir < GetNumberOfDirs())
d7c6ab14 83 {
04a02430 84
85 if (!OpenFile(fCurrentDir))
d7c6ab14 86 {
04a02430 87 delete fRunLoader; //close current session
88 fRunLoader = 0; //assure pointer is null
d7c6ab14 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
04a02430 102 Info("ReadNext","Reading Event %d",fCurrentEvent);
d7c6ab14 103 fRunLoader->GetEvent(fCurrentEvent);
104 AliStack* stack = fRunLoader->Stack();
105 if (!stack)
106 {
04a02430 107 Error("ReadNext","Can't get stack for event %d",fCurrentEvent);
d7c6ab14 108 continue;
109 }
110
04a02430 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="";
d7073a83 119 Int_t gentype=-1;
04a02430 120 AliHeader *header=fRunLoader->GetHeader();
121 if(!header) {
122 Warning("ReadNext","Header not found in event %d",fCurrentEvent);
123 } else {
d7073a83 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;
042d3031 158#ifndef NOUQHEADERINFO
d7073a83 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 }
04a02430 175 }
042d3031 176#endif
d7073a83 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);
04a02430 197 }
198 }
199 }
301a24f1 200 headdesc+=" Run ";
d7073a83 201 headdesc+=header->GetRun();
202 headdesc+=": Ev ";
203 headdesc+=header->GetEventNrInRun();
04a02430 204 }
205 fEventParticles->SetHeader(headdesc);
206
d7c6ab14 207 //get vertex
208 const TParticle *kv = stack->Particle(0);
209 if(kv) {
210 fEventParticles->SetVertex(kv->Vx(),kv->Vy(),kv->Vz());
211 }
212
04a02430 213 //loop over particles
d7c6ab14 214 for (Int_t i = 0;i<npart; i++)
215 {
216 TParticle *p = stack->Particle(i);
217 if(!p) continue;
b2760c9e 218 Int_t child1 = p->GetFirstDaughter();
219 //Int_t child2 = p->GetLastDaughter();
220 //Int_t mother = p->GetFirstMother();
b2760c9e 221 //cout << child1 << " " << child2 << " " << mother << endl;
04a02430 222 if((child1>=0) && (child1<nprim)) continue;
d7073a83 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;
301a24f1 233 p->SetWeight(-123); //mark particle
04a02430 234 }
d7073a83 235
236 //kinematic cuts
237 if(IsAcceptedParticle(p)){ //put particle in event
238 //if(p->Pt()>20){cout <<p->GetStatusCode() << " ";p->Print(); }
d7c6ab14 239 fEventParticles->AddParticle(p,i);
d7073a83 240 }
d7c6ab14 241 }
d7c6ab14 242 fCurrentEvent++;
243 fNEventsRead++;
244 return kTRUE;
b2760c9e 245 }
04a02430 246
247 //end of loop over directories specified in fDirs ObjArray
d7c6ab14 248 return kFALSE;
249}
250
d7c6ab14 251Int_t AliJetParticlesReaderKine::OpenFile(Int_t n)
252{
253 //opens file with kine tree
254
04a02430 255 if(fRunLoader){
256 if(fCurrentEvent < fRunLoader->GetNumberOfEvents()) return kTRUE;
257 else return kFALSE;
258 }
259
d7c6ab14 260 const TString& dirname = GetDirName(n);
261 if (dirname == "")
262 {
04a02430 263 Error("OpenNextFile","Can't get directory name with index %d",n);
d7c6ab14 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 }
04a02430 283
d7c6ab14 284 if (fRunLoader->LoadKinematics())
285 {
286 Error("OpenNextFile","Error occured while loading kinematics.");
287 return kFALSE;
288 }
04a02430 289
d7c6ab14 290 fCurrentEvent = 0;
291 return kTRUE;
292}