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