]>
Commit | Line | Data |
---|---|---|
a5556ea5 | 1 | #include "AliReaderESD.h" |
2 | //____________________________________________________________________ | |
3 | ////////////////////////////////////////////////////////////////////// | |
4 | // // | |
5 | // class AliReaderESD // | |
6 | // // | |
7 | // reader for ALICE Event Summary Data (ESD). // | |
8 | // // | |
9 | // Piotr.Skowronski@cern.ch // | |
10 | // // | |
11 | ////////////////////////////////////////////////////////////////////// | |
12 | ||
13 | #include <TPDGCode.h> | |
14 | #include <TString.h> | |
15 | #include <TObjString.h> | |
16 | #include <TTree.h> | |
17 | #include <TFile.h> | |
18 | #include <TKey.h> | |
19 | #include <TParticle.h> | |
20 | #include <TH1.h> | |
21 | ||
22 | #include <AliRun.h> | |
23 | #include <AliRunLoader.h> | |
24 | #include <AliStack.h> | |
25 | #include <AliESDtrack.h> | |
7814b246 | 26 | #include <AliKalmanTrack.h> |
a5556ea5 | 27 | #include <AliESD.h> |
28 | ||
29 | #include "AliAnalysis.h" | |
30 | #include "AliAODRun.h" | |
31 | #include "AliAOD.h" | |
afa8b37b | 32 | #include "AliAODParticle.h" |
a5556ea5 | 33 | #include "AliAODParticleCut.h" |
34 | #include "AliTrackPoints.h" | |
35 | #include "AliClusterMap.h" | |
36 | ||
37 | ClassImp(AliReaderESD) | |
38 | ||
39 | AliReaderESD::AliReaderESD(const Char_t* esdfilename, const Char_t* galfilename): | |
40 | fESDFileName(esdfilename), | |
41 | fGAlFileName(galfilename), | |
42 | fFile(0x0), | |
43 | fRunLoader(0x0), | |
44 | fKeyIterator(0x0), | |
45 | fReadSim(kFALSE), | |
46 | fCheckParticlePID(kFALSE), | |
47 | fReadMostProbableOnly(kFALSE), | |
48 | fNTrackPoints(0), | |
49 | fdR(0.0), | |
50 | fClusterMap(kFALSE), | |
7814b246 | 51 | fITSTrackPoints(kFALSE), |
52 | fMustTPC(kFALSE), | |
e1a64564 | 53 | fReadCentralBarrel(kTRUE), |
beb1c41d | 54 | fReadMuon(kFALSE), |
55 | fReadPHOS(kFALSE), | |
a5556ea5 | 56 | fNTPCClustMin(0), |
57 | fNTPCClustMax(150), | |
58 | fTPCChi2PerClustMin(0.0), | |
59 | fTPCChi2PerClustMax(10e5), | |
60 | fChi2Min(0.0), | |
61 | fChi2Max(10e5), | |
62 | fC00Min(0.0), | |
63 | fC00Max(10e5), | |
64 | fC11Min(0.0), | |
65 | fC11Max(10e5), | |
66 | fC22Min(0.0), | |
67 | fC22Max(10e5), | |
68 | fC33Min(0.0), | |
69 | fC33Max(10e5), | |
70 | fC44Min(0.0), | |
71 | fC44Max(10e5), | |
72 | fTPCC00Min(0.0), | |
73 | fTPCC00Max(10e5), | |
74 | fTPCC11Min(0.0), | |
75 | fTPCC11Max(10e5), | |
76 | fTPCC22Min(0.0), | |
77 | fTPCC22Max(10e5), | |
78 | fTPCC33Min(0.0), | |
79 | fTPCC33Max(10e5), | |
80 | fTPCC44Min(0.0), | |
81 | fTPCC44Max(10e5) | |
82 | ||
83 | { | |
84 | //cosntructor | |
85 | if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES)) | |
86 | Fatal("AliReaderESD","ESD defintions probobly changed. Ask Youra."); | |
87 | } | |
88 | /********************************************************************/ | |
89 | ||
90 | AliReaderESD::AliReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char_t* galfilename): | |
91 | AliReader(dirs), | |
92 | fESDFileName(esdfilename), | |
93 | fGAlFileName(galfilename), | |
94 | fFile(0x0), | |
95 | fRunLoader(0x0), | |
96 | fKeyIterator(0x0), | |
97 | fReadSim(kFALSE), | |
98 | fCheckParticlePID(kFALSE), | |
99 | fReadMostProbableOnly(kFALSE), | |
100 | fNTrackPoints(0), | |
101 | fdR(0.0), | |
102 | fClusterMap(kFALSE), | |
7814b246 | 103 | fITSTrackPoints(kFALSE), |
104 | fMustTPC(kFALSE), | |
e1a64564 | 105 | fReadCentralBarrel(kTRUE), |
beb1c41d | 106 | fReadMuon(kFALSE), |
107 | fReadPHOS(kFALSE), | |
a5556ea5 | 108 | fNTPCClustMin(0), |
109 | fNTPCClustMax(150), | |
110 | fTPCChi2PerClustMin(0.0), | |
111 | fTPCChi2PerClustMax(10e5), | |
112 | fChi2Min(0.0), | |
113 | fChi2Max(10e5), | |
114 | fC00Min(0.0), | |
115 | fC00Max(10e5), | |
116 | fC11Min(0.0), | |
117 | fC11Max(10e5), | |
118 | fC22Min(0.0), | |
119 | fC22Max(10e5), | |
120 | fC33Min(0.0), | |
121 | fC33Max(10e5), | |
122 | fC44Min(0.0), | |
123 | fC44Max(10e5), | |
124 | fTPCC00Min(0.0), | |
125 | fTPCC00Max(10e5), | |
126 | fTPCC11Min(0.0), | |
127 | fTPCC11Max(10e5), | |
128 | fTPCC22Min(0.0), | |
129 | fTPCC22Max(10e5), | |
130 | fTPCC33Min(0.0), | |
131 | fTPCC33Max(10e5), | |
132 | fTPCC44Min(0.0), | |
133 | fTPCC44Max(10e5) | |
134 | { | |
135 | //cosntructor | |
136 | if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES)) | |
137 | Fatal("AliReaderESD","ESD defintions probobly changed. Ask Youra."); | |
138 | } | |
139 | /********************************************************************/ | |
140 | ||
141 | AliReaderESD::~AliReaderESD() | |
142 | { | |
143 | //desctructor | |
beb1c41d | 144 | delete fRunLoader; |
145 | delete fKeyIterator; | |
146 | delete fFile; | |
a5556ea5 | 147 | } |
beb1c41d | 148 | |
a5556ea5 | 149 | /**********************************************************/ |
150 | Int_t AliReaderESD::ReadNext() | |
151 | { | |
152 | //reads next event from fFile | |
153 | //fRunLoader is for reading Kine | |
154 | ||
afa8b37b | 155 | if (AliVAODParticle::GetDebug()) |
a5556ea5 | 156 | Info("ReadNext","Entered"); |
157 | ||
158 | if (fEventSim == 0x0) fEventSim = new AliAOD(); | |
159 | if (fEventRec == 0x0) fEventRec = new AliAOD(); | |
160 | ||
161 | fEventSim->Reset(); | |
162 | fEventRec->Reset(); | |
163 | ||
164 | do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./" | |
beb1c41d | 165 | { |
166 | if (fFile == 0x0) | |
e1a64564 | 167 | { |
168 | fFile = OpenFile(fCurrentDir);//rl is opened here | |
169 | if (fFile == 0x0) | |
170 | { | |
171 | Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir); | |
172 | fCurrentDir++; | |
173 | continue; | |
174 | } | |
175 | fCurrentEvent = 0; | |
176 | } | |
a5556ea5 | 177 | TString esdname = "ESD"; |
178 | esdname+=fCurrentEvent; | |
179 | AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname)); | |
180 | if (esd == 0x0) | |
181 | { | |
afa8b37b | 182 | if (AliVAODParticle::GetDebug() > 2 ) |
a5556ea5 | 183 | { |
184 | Info("ReadNext","Can not find AliESD object named %s",esdname.Data()); | |
185 | } | |
186 | fCurrentDir++; | |
a5556ea5 | 187 | delete fFile;//we have to assume there is no more ESD objects in the fFile |
188 | fFile = 0x0; | |
189 | delete fRunLoader; | |
190 | fRunLoader = 0x0; | |
191 | continue; | |
192 | } | |
beb1c41d | 193 | ReadESD(esd); |
a5556ea5 | 194 | |
beb1c41d | 195 | fCurrentEvent++; |
196 | fNEventsRead++; | |
197 | delete esd; | |
198 | return 0;//success -> read one event | |
199 | }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array | |
a5556ea5 | 200 | |
201 | return 1; //no more directories to read | |
202 | } | |
a5556ea5 | 203 | |
beb1c41d | 204 | /**********************************************************/ |
a5556ea5 | 205 | Int_t AliReaderESD::ReadESD(AliESD* esd) |
beb1c41d | 206 | { |
207 | //Reads esd data | |
208 | if (esd == 0x0) | |
209 | { | |
210 | Error("ReadESD","ESD is NULL"); | |
211 | return 1; | |
212 | } | |
213 | ||
214 | // seperate each method | |
215 | if (fReadCentralBarrel) ReadESDCentral(esd); | |
216 | ||
217 | if (fReadMuon) ReadESDMuon(esd); | |
218 | ||
219 | if (fReadPHOS) ReadESDPHOS(esd); | |
220 | ||
221 | return 1; | |
222 | } | |
223 | ||
224 | /**********************************************************/ | |
225 | Int_t AliReaderESD::ReadESDCentral(AliESD* esd) | |
a5556ea5 | 226 | { |
227 | //****** Tentative particle type "concentrations" | |
228 | static const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05}; | |
229 | ||
230 | Double_t pidtable[kNSpecies];//array used for reading pid probabilities from ESD track | |
231 | Double_t w[kNSpecies]; | |
232 | Double_t mom[3];//momentum | |
233 | Double_t pos[3];//position | |
234 | Double_t vertexpos[3];//vertex position | |
235 | //Reads one ESD | |
beb1c41d | 236 | |
a5556ea5 | 237 | TDatabasePDG* pdgdb = TDatabasePDG::Instance(); |
238 | if (pdgdb == 0x0) | |
239 | { | |
240 | Error("ReadESD","Can not get PDG Database Instance."); | |
241 | return 1; | |
242 | } | |
243 | ||
7257b3b6 | 244 | Float_t mf = esd->GetMagneticField()/10.; //AliESD::GetMagnField returns mf in kG |
a5556ea5 | 245 | |
7814b246 | 246 | if ( (mf == 0.0) && ((fNTrackPoints > 0) || fITSTrackPoints) ) |
a5556ea5 | 247 | { |
248 | Error("ReadESD","Magnetic Field is 0 and Track Points Demended. Skipping to next event."); | |
249 | return 1; | |
250 | } | |
251 | ||
7814b246 | 252 | if (fITSTrackPoints) |
253 | { | |
7257b3b6 | 254 | Info("ReadESD","Magnetic Field is %f",mf); |
255 | AliKalmanTrack::SetMagneticField(mf); | |
7814b246 | 256 | } |
257 | ||
a5556ea5 | 258 | AliStack* stack = 0x0; |
259 | if (fReadSim && fRunLoader) | |
260 | { | |
261 | fRunLoader->GetEvent(fCurrentEvent); | |
262 | stack = fRunLoader->Stack(); | |
263 | } | |
264 | ||
265 | const AliESDVertex* vertex = esd->GetVertex(); | |
266 | if (vertex == 0x0) | |
267 | { | |
268 | Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)"); | |
269 | vertexpos[0] = 0.0; | |
270 | vertexpos[1] = 0.0; | |
271 | vertexpos[2] = 0.0; | |
272 | } | |
273 | else | |
274 | { | |
275 | vertex->GetXYZ(vertexpos); | |
276 | } | |
277 | ||
afa8b37b | 278 | if (AliVAODParticle::GetDebug() > 0) |
a5556ea5 | 279 | { |
280 | Info("ReadESD","Primary Vertex is (%f,%f,%f)",vertexpos[0],vertexpos[1],vertexpos[2]); | |
281 | } | |
282 | ||
283 | Info("ReadESD","Reading Event %d",fCurrentEvent); | |
284 | ||
285 | Int_t ntr = esd->GetNumberOfTracks(); | |
286 | Info("ReadESD","Found %d tracks.",ntr); | |
287 | for (Int_t i = 0;i<ntr; i++) | |
288 | { | |
289 | AliESDtrack *esdtrack = esd->GetTrack(i); | |
290 | if (esdtrack == 0x0) | |
291 | { | |
292 | Error("Next","Can not get track %d", i); | |
293 | continue; | |
294 | } | |
295 | ||
296 | //if (esdtrack->HasVertexParameters() == kFALSE) | |
297 | if ((esdtrack->GetStatus() & AliESDtrack::kITSrefit) == kFALSE) | |
298 | { | |
afa8b37b | 299 | if (AliVAODParticle::GetDebug() > 2) |
a5556ea5 | 300 | Info("ReadNext","Particle skipped: Data at vertex not available."); |
301 | continue; | |
302 | } | |
303 | ||
7814b246 | 304 | if (fMustTPC) |
305 | { | |
306 | if ((esdtrack->GetStatus() & AliESDtrack::kTPCin) == kFALSE) | |
307 | { | |
308 | if (AliVAODParticle::GetDebug() > 2) | |
309 | Info("ReadNext","Particle skipped: Was not reconstructed in TPC."); | |
310 | continue; | |
311 | } | |
312 | } | |
313 | ||
a5556ea5 | 314 | if ((esdtrack->GetStatus() & AliESDtrack::kESDpid) == kFALSE) |
315 | { | |
afa8b37b | 316 | if (AliVAODParticle::GetDebug() > 2) |
a5556ea5 | 317 | Info("ReadNext","Particle skipped: PID BIT is not set."); |
318 | continue; | |
319 | } | |
320 | ||
321 | ||
322 | Double_t extx; | |
323 | Double_t extp[5]; | |
324 | esdtrack->GetConstrainedExternalParameters(extx,extp); | |
325 | if (extp[4] == 0.0) | |
326 | { | |
afa8b37b | 327 | if (AliVAODParticle::GetDebug() > 2) |
a5556ea5 | 328 | Info("ReadNext","Track has 0 contrianed curvature -> Probobly parameters never updated. Skipping."); |
329 | continue; | |
330 | } | |
331 | esdtrack->GetESDpid(pidtable); | |
332 | esdtrack->GetConstrainedPxPyPz(mom); | |
333 | esdtrack->GetConstrainedXYZ(pos); | |
334 | pos[0] -= vertexpos[0];//we are interested only in relative position to Primary vertex at this point | |
335 | pos[1] -= vertexpos[1]; | |
336 | pos[2] -= vertexpos[2]; | |
337 | ||
338 | Int_t charge = (extp[4] > 0)?1:-1;//if curvature=charg/Pt is positive charge is positive | |
339 | ||
340 | //Particle from kinematics | |
afa8b37b | 341 | AliAODParticle* particle = 0; |
a5556ea5 | 342 | Bool_t keeppart = kFALSE; |
343 | if ( fReadSim && stack ) | |
344 | { | |
345 | if (esdtrack->GetLabel() < 0) continue;//this is fake - we are not able to match any track | |
346 | TParticle *p = stack->Particle(esdtrack->GetLabel()); | |
347 | if (p==0x0) | |
348 | { | |
349 | Error("ReadNext","Can not find track with such label."); | |
350 | continue; | |
351 | } | |
7814b246 | 352 | |
353 | if (fCheckParticlePID) | |
a5556ea5 | 354 | { |
cea0a066 | 355 | if(Rejected(p->GetPdgCode())) |
7814b246 | 356 | { |
357 | if ( AliVAODParticle::GetDebug() > 5 ) | |
358 | Info("ReadNext","Simulated Particle PID (%d) did not pass the cut.",p->GetPdgCode()); | |
359 | continue; //check if we are intersted with particles of this type | |
360 | } | |
361 | } | |
a5556ea5 | 362 | // if(p->GetPdgCode()<0) charge = -1; |
afa8b37b | 363 | particle = new AliAODParticle(*p,i); |
a5556ea5 | 364 | |
365 | } | |
366 | ||
367 | if(CheckTrack(esdtrack)) continue; | |
368 | ||
369 | //Here we apply Bayes' formula | |
370 | Double_t rc=0.; | |
371 | for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s]; | |
372 | if (rc==0.0) | |
373 | { | |
afa8b37b | 374 | if (AliVAODParticle::GetDebug() > 2) |
a5556ea5 | 375 | Info("ReadNext","Particle rejected since total bayessian PID probab. is zero."); |
376 | continue; | |
377 | } | |
378 | ||
379 | for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc; | |
380 | ||
afa8b37b | 381 | if (AliVAODParticle::GetDebug() > 4) |
a5556ea5 | 382 | { |
383 | Info("ReadNext","###########################################################################"); | |
384 | Info("ReadNext","Momentum: %f %f %f",mom[0],mom[1],mom[2]); | |
385 | Info("ReadNext","Position: %f %f %f",pos[0],pos[1],pos[2]); | |
386 | TString msg("Pid list got from track:"); | |
387 | for (Int_t s = 0;s<kNSpecies;s++) | |
388 | { | |
389 | msg+="\n "; | |
390 | msg+=s; | |
391 | msg+="("; | |
392 | msg+=charge*GetSpeciesPdgCode((ESpecies)s); | |
393 | msg+="): "; | |
394 | msg+=w[s]; | |
395 | msg+=" ("; | |
396 | msg+=pidtable[s]; | |
397 | msg+=")"; | |
398 | } | |
399 | Info("ReadNext","%s",msg.Data()); | |
afa8b37b | 400 | }//if (AliVAODParticle::GetDebug()>4) |
a5556ea5 | 401 | |
402 | AliTrackPoints* tpts = 0x0; | |
403 | if (fNTrackPoints > 0) | |
404 | { | |
7257b3b6 | 405 | tpts = new AliTrackPoints(fNTrackPoints,esdtrack,mf*10.0,fdR); |
0d8a4589 | 406 | // tpts->Move(-vertexpos[0],-vertexpos[1],-vertexpos[2]); |
7814b246 | 407 | } |
408 | ||
409 | AliTrackPoints* itstpts = 0x0; | |
410 | if (fITSTrackPoints) | |
411 | { | |
7257b3b6 | 412 | itstpts = new AliTrackPoints(AliTrackPoints::kITS,esdtrack,mf*10.0); |
0d8a4589 | 413 | // itstpts->Move(-vertexpos[0],-vertexpos[1],-vertexpos[2]); |
a5556ea5 | 414 | } |
415 | ||
416 | AliClusterMap* cmap = 0x0; | |
417 | if ( fClusterMap ) | |
418 | { | |
419 | cmap = new AliClusterMap(esdtrack); | |
420 | } | |
421 | ||
422 | //If this flag fReadMostProbableOnly is false the | |
423 | //loop over species (see "LOOPPIDS") is over all possible PIDs | |
424 | //in other case the most probablle PID is searched | |
425 | //and the loop is limited to that PID only | |
426 | ||
427 | Int_t firstspecie = 0; | |
428 | Int_t lastspecie = kNSpecies; | |
429 | ||
430 | if (fReadMostProbableOnly) | |
431 | { | |
432 | //find the most probable PID | |
433 | Int_t spec = 0; | |
434 | Float_t maxprob = w[0]; | |
435 | for (Int_t s=1; s<AliESDtrack::kSPECIES; s++) | |
436 | { | |
437 | if (w[s]>maxprob) | |
438 | { | |
439 | maxprob = w[s]; | |
440 | spec = s; | |
441 | } | |
442 | } | |
443 | firstspecie = spec; | |
444 | lastspecie = spec + 1; | |
445 | } | |
446 | ||
447 | for (Int_t s = firstspecie; s<lastspecie; s++)//LOOPPIDS | |
448 | { | |
449 | Int_t pdgcode = charge*GetSpeciesPdgCode((ESpecies)s); | |
450 | Float_t pp = w[s]; | |
451 | if (pp == 0.0) | |
452 | { | |
afa8b37b | 453 | if ( AliVAODParticle::GetDebug() > 5 ) |
a5556ea5 | 454 | Info("ReadNext","Probability of being PID %d is zero. Continuing.",pdgcode); |
455 | continue; | |
456 | } | |
457 | ||
cea0a066 | 458 | if(Rejected(pdgcode)) |
a5556ea5 | 459 | { |
afa8b37b | 460 | if ( AliVAODParticle::GetDebug() > 5 ) |
a5556ea5 | 461 | Info("ReadNext","PID (%d) did not pass the cut.",pdgcode); |
462 | continue; //check if we are intersted with particles of this type | |
463 | } | |
464 | ||
465 | Double_t mass = pdgdb->GetParticle(pdgcode)->Mass(); | |
466 | Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track | |
467 | ||
afa8b37b | 468 | AliAODParticle* track = new AliAODParticle(pdgcode, w[s],i, |
a5556ea5 | 469 | mom[0], mom[1], mom[2], tEtot, |
470 | pos[0], pos[1], pos[2], 0.); | |
471 | //copy probabilitis of other species (if not zero) | |
472 | for (Int_t k = 0; k<kNSpecies; k++) | |
473 | { | |
474 | if (k == s) continue; | |
475 | if (w[k] == 0.0) continue; | |
476 | track->SetPIDprobability(charge*GetSpeciesPdgCode( (ESpecies)k ),w[k]); | |
477 | } | |
478 | ||
cea0a066 | 479 | if(Rejected(track))//check if meets all criteria of any of our cuts |
a5556ea5 | 480 | //if it does not delete it and take next good track |
481 | { | |
afa8b37b | 482 | if ( AliVAODParticle::GetDebug() > 4 ) |
a5556ea5 | 483 | Info("ReadNext","Track did not pass the cut"); |
484 | delete track; | |
485 | continue; | |
486 | } | |
487 | ||
488 | //Single Particle cuts on cluster map and track points rather do not have sense | |
489 | if (tpts) | |
490 | { | |
78d7c6d3 | 491 | track->SetTPCTrackPoints(tpts); |
a5556ea5 | 492 | } |
7814b246 | 493 | |
494 | if (itstpts) | |
495 | { | |
496 | track->SetITSTrackPoints(itstpts); | |
497 | } | |
a5556ea5 | 498 | |
499 | if (cmap) | |
500 | { | |
501 | track->SetClusterMap(cmap); | |
502 | } | |
503 | ||
504 | fEventRec->AddParticle(track); | |
505 | if (particle) fEventSim->AddParticle(particle); | |
506 | keeppart = kTRUE; | |
507 | ||
afa8b37b | 508 | if (AliVAODParticle::GetDebug() > 4 ) |
a5556ea5 | 509 | { |
510 | Info("ReadNext","\n\nAdding Particle with incarnation %d",pdgcode); | |
511 | track->Print(); | |
512 | if (particle) particle->Print(); | |
513 | Info("ReadNext","\n----------------------------------------------\n"); | |
514 | } | |
515 | }//for (Int_t s = 0; s<kNSpecies; s++) | |
516 | ||
517 | if (keeppart == kFALSE) | |
518 | { | |
519 | delete particle;//particle was not stored in event | |
520 | delete tpts; | |
7814b246 | 521 | delete itstpts; |
a5556ea5 | 522 | delete cmap; |
523 | } | |
7814b246 | 524 | else |
525 | { | |
f4b459f7 | 526 | if ( fReadSim && stack ) |
7814b246 | 527 | { |
f4b459f7 | 528 | if (particle->P() < 0.00001) |
529 | { | |
530 | Info("ReadNext","###################################"); | |
531 | Info("ReadNext","###################################"); | |
532 | Info("ReadNext","Track Label %d",esdtrack->GetLabel()); | |
533 | TParticle *p = stack->Particle(esdtrack->GetLabel()); | |
534 | Info("ReadNext",""); | |
535 | p->Print(); | |
536 | Info("ReadNext",""); | |
537 | particle->Print(); | |
538 | } | |
539 | } | |
7814b246 | 540 | } |
a5556ea5 | 541 | |
542 | }//for (Int_t i = 0;i<ntr; i++) -- loop over tracks | |
543 | ||
544 | Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).", | |
545 | fEventRec->GetNumberOfParticles(), fEventSim->GetNumberOfParticles(), | |
546 | fNEventsRead,fCurrentEvent,fCurrentDir); | |
547 | fTrackCounter->Fill(fEventRec->GetNumberOfParticles()); | |
0d8a4589 | 548 | |
549 | /******************************************************/ | |
550 | /****** Setting glevet properties *************/ | |
551 | /******************************************************/ | |
552 | ||
553 | if (fEventRec->GetNumberOfParticles() > 0) | |
554 | { | |
555 | fEventRec->SetPrimaryVertex(vertexpos[0],vertexpos[1],vertexpos[2]); | |
556 | } | |
a5556ea5 | 557 | return 0; |
558 | } | |
559 | ||
beb1c41d | 560 | /**********************************************************/ |
561 | Int_t AliReaderESD::ReadESDMuon(AliESD* esd) | |
562 | { | |
563 | ||
564 | Double_t vertexpos[3];//vertex position, assuming no secondary decay | |
565 | ||
566 | const AliESDVertex* vertex = esd->GetVertex(); | |
567 | ||
568 | if (vertex == 0x0) { | |
569 | Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)"); | |
570 | vertexpos[0] = 0.0; | |
571 | vertexpos[1] = 0.0; | |
572 | vertexpos[2] = 0.0; | |
573 | } else { | |
574 | vertex->GetXYZ(vertexpos); | |
575 | } | |
576 | ||
577 | Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ; | |
578 | ||
579 | if (AliVAODParticle::GetDebug() > 0) { | |
580 | Info("ReadESD","Reading Event %d",fCurrentEvent); | |
581 | Info("ReadESD","Found %d tracks.",nTracks); | |
582 | } | |
583 | // settings | |
584 | Float_t Chi2Cut = 100.; | |
585 | Float_t PtCutMin = 1.; | |
586 | Float_t PtCutMax = 10000.; | |
587 | Float_t muonMass = 0.105658389; | |
588 | Int_t pdgcode = -13; | |
589 | Double_t thetaX, thetaY, pYZ; | |
590 | Double_t pxRec1, pyRec1, pzRec1, E1; | |
591 | Int_t charge; | |
592 | ||
593 | Int_t ntrackhits; | |
594 | Double_t fitfmin; | |
595 | ||
596 | TLorentzVector fV1; | |
597 | fEventRec->Reset(); | |
598 | for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) { | |
599 | ||
600 | AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack); | |
601 | ||
602 | thetaX = muonTrack->GetThetaX(); | |
603 | thetaY = muonTrack->GetThetaY(); | |
604 | ||
605 | pYZ = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum()); | |
606 | pzRec1 = - pYZ / TMath::Sqrt(1.0 + TMath::Tan(thetaY)*TMath::Tan(thetaX)); | |
607 | pxRec1 = pzRec1 * TMath::Tan(thetaX); | |
608 | pyRec1 = pzRec1 * TMath::Tan(thetaY); | |
609 | charge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum())); | |
610 | E1 = TMath::Sqrt(muonMass * muonMass + pxRec1 * pxRec1 + pyRec1 * pyRec1 + pzRec1 * pzRec1); | |
611 | fV1.SetPxPyPzE(pxRec1, pyRec1, pzRec1, E1); | |
612 | ||
613 | ntrackhits = muonTrack->GetNHit(); | |
614 | fitfmin = muonTrack->GetChi2(); | |
615 | ||
616 | // transverse momentum | |
617 | Float_t pt1 = fV1.Pt(); | |
618 | ||
619 | // chi2 per d.o.f. | |
620 | Float_t ch1 = fitfmin / (2.0 * ntrackhits - 5); | |
621 | ||
622 | if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) { | |
623 | AliAODParticle* track = new AliAODParticle(pdgcode*charge,1,iTrack, | |
624 | pxRec1, pyRec1,pzRec1, E1, | |
625 | vertexpos[0], vertexpos[1], vertexpos[2], 0.); | |
626 | fEventRec->AddParticle(track); | |
627 | } | |
628 | ||
629 | } | |
630 | fTrackCounter->Fill(fEventRec->GetNumberOfParticles()); | |
631 | return 0; | |
632 | } | |
633 | ||
a5556ea5 | 634 | /**********************************************************/ |
635 | ||
636 | void AliReaderESD::Rewind() | |
637 | { | |
638 | //rewinds reading | |
beb1c41d | 639 | // delete fKeyIterator; |
a5556ea5 | 640 | delete fFile; |
641 | fFile = 0x0; | |
beb1c41d | 642 | // fKeyIterator = 0x0; |
a5556ea5 | 643 | delete fRunLoader; |
644 | fRunLoader = 0x0; | |
645 | fCurrentDir = 0; | |
646 | fNEventsRead = 0; | |
647 | if (fTrackCounter) fTrackCounter->Reset(); | |
648 | } | |
649 | /**********************************************************/ | |
650 | ||
651 | TFile* AliReaderESD::OpenFile(Int_t n) | |
652 | { | |
beb1c41d | 653 | //opens fFile with tree |
a5556ea5 | 654 | |
655 | const TString& dirname = GetDirName(n); | |
656 | if (dirname == "") | |
657 | { | |
658 | Error("OpenFiles","Can not get directory name"); | |
659 | return 0x0; | |
660 | } | |
661 | TString filename = dirname +"/"+ fESDFileName; | |
662 | TFile *ret = TFile::Open(filename.Data()); | |
663 | ||
664 | if ( ret == 0x0) | |
665 | { | |
666 | Error("OpenFiles","Can't open fFile %s",filename.Data()); | |
667 | return 0x0; | |
668 | } | |
669 | if (!ret->IsOpen()) | |
670 | { | |
671 | Error("OpenFiles","Can't open fFile %s",filename.Data()); | |
672 | return 0x0; | |
673 | } | |
674 | ||
675 | if (fReadSim ) | |
676 | { | |
677 | fRunLoader = AliRunLoader::Open(dirname +"/"+ fGAlFileName); | |
678 | if (fRunLoader == 0x0) | |
679 | { | |
680 | Error("OpenFiles","Can't get RunLoader for directory %s",dirname.Data()); | |
681 | delete ret; | |
682 | return 0x0; | |
683 | } | |
684 | ||
685 | fRunLoader->LoadHeader(); | |
686 | if (fRunLoader->LoadKinematics()) | |
687 | { | |
688 | Error("Next","Error occured while loading kinematics."); | |
689 | delete fRunLoader; | |
690 | delete ret; | |
691 | return 0x0; | |
692 | } | |
693 | } | |
694 | ||
695 | return ret; | |
696 | } | |
697 | /**********************************************************/ | |
698 | ||
699 | Int_t AliReaderESD::GetSpeciesPdgCode(ESpecies spec)//skowron | |
700 | { | |
701 | //returns pdg code from the PID index | |
702 | //ask jura about charge | |
703 | switch (spec) | |
704 | { | |
705 | case kESDElectron: | |
706 | return kPositron; | |
707 | break; | |
708 | case kESDMuon: | |
709 | return kMuonPlus; | |
710 | break; | |
711 | case kESDPion: | |
712 | return kPiPlus; | |
713 | break; | |
714 | case kESDKaon: | |
715 | return kKPlus; | |
716 | break; | |
717 | case kESDProton: | |
718 | return kProton; | |
719 | break; | |
720 | default: | |
721 | ::Warning("GetSpeciesPdgCode","Specie with number %d is not defined.",(Int_t)spec); | |
722 | break; | |
723 | } | |
724 | return 0; | |
725 | } | |
726 | /********************************************************************/ | |
727 | Bool_t AliReaderESD::CheckTrack(AliESDtrack* t) const | |
728 | { | |
729 | //Performs check of the track | |
730 | ||
731 | if ( (t->GetConstrainedChi2() < fChi2Min) || (t->GetConstrainedChi2() > fChi2Max) ) return kTRUE; | |
732 | ||
733 | if ( (t->GetTPCclusters(0x0) < fNTPCClustMin) || (t->GetTPCclusters(0x0) > fNTPCClustMax) ) return kTRUE; | |
734 | ||
735 | if (t->GetTPCclusters(0x0) > 0) | |
736 | { | |
737 | Float_t chisqpercl = t->GetTPCchi2()/((Double_t)t->GetTPCclusters(0x0)); | |
738 | if ( (chisqpercl < fTPCChi2PerClustMin) || (chisqpercl > fTPCChi2PerClustMax) ) return kTRUE; | |
739 | } | |
740 | ||
741 | Double_t cc[15]; | |
742 | t->GetConstrainedExternalCovariance(cc); | |
743 | ||
744 | if ( (cc[0] < fC00Min) || (cc[0] > fC00Max) ) return kTRUE; | |
745 | if ( (cc[2] < fC11Min) || (cc[2] > fC11Max) ) return kTRUE; | |
746 | if ( (cc[5] < fC22Min) || (cc[5] > fC22Max) ) return kTRUE; | |
747 | if ( (cc[9] < fC33Min) || (cc[9] > fC33Max) ) return kTRUE; | |
748 | if ( (cc[14] < fC44Min) || (cc[14] > fC44Max) ) return kTRUE; | |
749 | ||
750 | ||
751 | t->GetInnerExternalCovariance(cc); | |
752 | ||
753 | if ( (cc[0] < fTPCC00Min) || (cc[0] > fTPCC00Max) ) return kTRUE; | |
754 | if ( (cc[2] < fTPCC11Min) || (cc[2] > fTPCC11Max) ) return kTRUE; | |
755 | if ( (cc[5] < fTPCC22Min) || (cc[5] > fTPCC22Max) ) return kTRUE; | |
756 | if ( (cc[9] < fTPCC33Min) || (cc[9] > fTPCC33Max) ) return kTRUE; | |
757 | if ( (cc[14] < fTPCC44Min) || (cc[14] > fTPCC44Max) ) return kTRUE; | |
758 | ||
759 | return kFALSE; | |
760 | ||
761 | } | |
762 | /********************************************************************/ | |
763 | ||
764 | void AliReaderESD::SetChi2Range(Float_t min, Float_t max) | |
765 | { | |
766 | //sets range of Chi2 per Cluster | |
767 | fChi2Min = min; | |
768 | fChi2Max = max; | |
769 | } | |
770 | /********************************************************************/ | |
771 | ||
772 | void AliReaderESD::SetTPCNClustersRange(Int_t min,Int_t max) | |
773 | { | |
774 | //sets range of Number Of Clusters that tracks have to have | |
775 | fNTPCClustMin = min; | |
776 | fNTPCClustMax = max; | |
777 | } | |
778 | /********************************************************************/ | |
779 | ||
780 | void AliReaderESD::SetTPCChi2PerCluserRange(Float_t min, Float_t max) | |
781 | { | |
782 | //sets range of Chi2 per Cluster | |
783 | fTPCChi2PerClustMin = min; | |
784 | fTPCChi2PerClustMax = max; | |
785 | } | |
786 | /********************************************************************/ | |
787 | ||
788 | void AliReaderESD::SetC00Range(Float_t min, Float_t max) | |
789 | { | |
790 | //Sets range of C00 parameter of covariance matrix of the track | |
791 | //it defines uncertainty of the momentum | |
792 | fC00Min = min; | |
793 | fC00Max = max; | |
794 | } | |
795 | /********************************************************************/ | |
796 | ||
797 | void AliReaderESD::SetC11Range(Float_t min, Float_t max) | |
798 | { | |
799 | //Sets range of C11 parameter of covariance matrix of the track | |
800 | //it defines uncertainty of the momentum | |
801 | fC11Min = min; | |
802 | fC11Max = max; | |
803 | } | |
804 | /********************************************************************/ | |
805 | ||
806 | void AliReaderESD::SetC22Range(Float_t min, Float_t max) | |
807 | { | |
808 | //Sets range of C22 parameter of covariance matrix of the track | |
809 | //it defines uncertainty of the momentum | |
810 | fC22Min = min; | |
811 | fC22Max = max; | |
812 | } | |
813 | /********************************************************************/ | |
814 | ||
815 | void AliReaderESD::SetC33Range(Float_t min, Float_t max) | |
816 | { | |
817 | //Sets range of C33 parameter of covariance matrix of the track | |
818 | //it defines uncertainty of the momentum | |
819 | fC33Min = min; | |
820 | fC33Max = max; | |
821 | } | |
822 | /********************************************************************/ | |
823 | ||
824 | void AliReaderESD::SetC44Range(Float_t min, Float_t max) | |
825 | { | |
826 | //Sets range of C44 parameter of covariance matrix of the track | |
827 | //it defines uncertainty of the momentum | |
828 | fC44Min = min; | |
829 | fC44Max = max; | |
830 | } |