]>
Commit | Line | Data |
---|---|---|
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",fCurrentDir*1000+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 | } |