]>
Commit | Line | Data |
---|---|---|
88b71b9f | 1 | #include "TChain.h" |
2 | #include "TF1.h" | |
3 | #include "TRandom3.h" | |
4 | ||
5 | #include "AliMCParticle.h" | |
6 | ||
7 | #include "AliAnalysisTask.h" | |
8 | #include "AliAnalysisManager.h" | |
9 | ||
10 | #include "AliVEvent.h" | |
11 | #include "AliESDEvent.h" | |
12 | #include "AliAODEvent.h" | |
13 | #include "AliMCEvent.h" | |
14 | #include "AliESDInputHandler.h" | |
15 | #include "AliInputEventHandler.h" | |
16 | #include "AliVTrack.h" | |
17 | #include "AliExternalTrackParam.h" | |
18 | #include "AliVVertex.h" | |
19 | #include "AliAnalysisFilter.h" | |
20 | #include "AliPID.h" | |
21 | #include "AliPIDResponse.h" | |
22 | #include "AliESDv0KineCuts.h" | |
23 | #include "AliESDv0.h" | |
24 | ||
25 | #include "AliTPCPIDBase.h" | |
26 | ||
27 | /* | |
28 | This class is a base class for all other | |
29 | analysis tasks that use V0's. | |
30 | It provides basics for V0 identification. | |
31 | In addition, some further basic functions are provided. | |
32 | ||
33 | Class written by Benjamin Hess. | |
34 | Contact: bhess@cern.ch | |
35 | */ | |
36 | ||
37 | ClassImp(AliTPCPIDBase) | |
38 | ||
39 | Double_t AliTPCPIDBase::fgCutGeo = 1.; | |
40 | Double_t AliTPCPIDBase::fgCutNcr = 0.85; | |
41 | Double_t AliTPCPIDBase::fgCutNcl = 0.7; | |
42 | ||
43 | UShort_t AliTPCPIDBase::fgCutPureNcl = 60; | |
44 | ||
45 | //________________________________________________________________________ | |
46 | AliTPCPIDBase::AliTPCPIDBase() | |
47 | : AliAnalysisTaskSE() | |
48 | , fEvent(0x0) | |
49 | , fESD(0x0) | |
50 | , fMC(0x0) | |
51 | , fPIDResponse(0x0) | |
52 | , fV0KineCuts(0x0) | |
53 | , fIsPbpOrpPb(kFALSE) | |
54 | , fUsePhiCut(kFALSE) | |
55 | , fTPCcutType(kNoCut) | |
56 | , fZvtxCutEvent(10.0) | |
57 | , fEtaCut(0.9) | |
58 | , fPhiCutLow(0x0) | |
59 | , fPhiCutHigh(0x0) | |
60 | , fRandom(0x0) | |
61 | , fTrackFilter(0x0) | |
62 | , fNumTagsStored(0) | |
63 | , fV0tags(0x0) | |
64 | , fStoreMotherIndex(kFALSE) | |
65 | , fV0motherIndex(0x0) | |
66 | { | |
67 | // default Constructor | |
68 | ||
69 | fRandom = new TRandom3(0); // 0 means random seed | |
70 | ||
71 | ||
72 | // Question: Is this the right place to initialize these functions? | |
73 | // Will it work on proof? i.e. will they be streamed to the workers? | |
74 | // Also one should add getters and setters | |
75 | fPhiCutLow = new TF1("StdPhiCutLow", "0.1/x/x+pi/18.0-0.025", 0, 100); | |
76 | fPhiCutHigh = new TF1("StandardPhiCutHigh", "0.12/x+pi/18.0+0.035", 0, 100); | |
77 | } | |
78 | ||
79 | //________________________________________________________________________ | |
80 | AliTPCPIDBase::AliTPCPIDBase(const char *name) | |
81 | : AliAnalysisTaskSE(name) | |
82 | , fEvent(0x0) | |
83 | , fESD(0x0) | |
84 | , fMC(0x0) | |
85 | , fPIDResponse(0x0) | |
86 | , fV0KineCuts(0x0) | |
87 | , fIsPbpOrpPb(kFALSE) | |
88 | , fUsePhiCut(kFALSE) | |
89 | , fTPCcutType(kNoCut) | |
90 | , fZvtxCutEvent(10.0) | |
91 | , fEtaCut(0.9) | |
92 | , fPhiCutLow(0x0) | |
93 | , fPhiCutHigh(0x0) | |
94 | , fRandom(0x0) | |
95 | , fTrackFilter(0x0) | |
96 | , fNumTagsStored(0) | |
97 | , fV0tags(0x0) | |
98 | , fStoreMotherIndex(kFALSE) | |
99 | , fV0motherIndex(0x0) | |
100 | { | |
101 | // Constructor | |
102 | ||
103 | fRandom = new TRandom3(0); // 0 means random seed | |
104 | ||
105 | ||
106 | // Question: Is this the right place to initialize these functions? | |
107 | // Will it work on proof? i.e. will they be streamed to the workers? | |
108 | // Also one should add getters and setters | |
109 | fPhiCutLow = new TF1("StdPhiCutLow", "0.1/x/x+pi/18.0-0.025", 0, 100); | |
110 | fPhiCutHigh = new TF1("StandardPhiCutHigh", "0.12/x+pi/18.0+0.035", 0, 100); | |
111 | ||
112 | DefineInput (0, TChain::Class()); | |
113 | DefineOutput(0, TTree::Class()); | |
114 | } | |
115 | ||
116 | ||
117 | //________________________________________________________________________ | |
118 | AliTPCPIDBase::~AliTPCPIDBase() | |
119 | { | |
120 | // dtor | |
121 | ||
122 | delete fPhiCutLow; | |
123 | fPhiCutLow = 0; | |
124 | ||
125 | delete fPhiCutHigh; | |
126 | fPhiCutHigh = 0; | |
127 | ||
128 | delete fTrackFilter; | |
129 | fTrackFilter = 0; | |
130 | ||
131 | delete fRandom; | |
132 | fRandom = 0; | |
133 | ||
134 | delete fV0KineCuts; | |
135 | fV0KineCuts = 0; | |
136 | ||
137 | delete fV0tags; | |
138 | fV0tags = 0; | |
139 | fNumTagsStored = 0; | |
140 | ||
141 | delete fV0motherIndex; | |
142 | fV0motherIndex = 0; | |
143 | } | |
144 | ||
145 | ||
146 | //________________________________________________________________________ | |
147 | void AliTPCPIDBase::UserCreateOutputObjects() | |
148 | { | |
149 | // Create histograms | |
150 | // Called once | |
151 | ||
152 | // Input hander | |
153 | AliAnalysisManager* man = AliAnalysisManager::GetAnalysisManager(); | |
154 | AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler()); | |
155 | ||
156 | if (!inputHandler) { | |
157 | AliFatal("Input handler needed"); | |
158 | fPIDResponse = 0x0; | |
159 | ||
160 | return; | |
161 | } | |
162 | ||
163 | // PID response object | |
164 | fPIDResponse = inputHandler->GetPIDResponse(); | |
165 | if (!fPIDResponse) | |
166 | AliError("PIDResponse object was not created"); | |
167 | ||
168 | // V0 Kine cuts | |
169 | fV0KineCuts = new AliESDv0KineCuts; | |
170 | fV0KineCuts->SetGammaCutChi2NDF(5.); | |
171 | // Only accept V0el with prod. radius within 45 cm -> PID will by systematically biased for larger values! | |
172 | Float_t gammaProdVertexRadiusCuts[2] = { 3.0, 45. }; | |
173 | fV0KineCuts->SetGammaCutVertexR(&gammaProdVertexRadiusCuts[0]); | |
174 | } | |
175 | ||
176 | ||
177 | //________________________________________________________________________ | |
178 | void AliTPCPIDBase::UserExec(Option_t *) | |
179 | { | |
180 | // Main loop | |
181 | // Called for each event | |
182 | } | |
183 | ||
184 | //________________________________________________________________________ | |
185 | void AliTPCPIDBase::Terminate(const Option_t *) | |
186 | { | |
187 | // Called once at the end of the query | |
188 | } | |
189 | ||
190 | ||
191 | //_____________________________________________________________________________ | |
192 | Double_t AliTPCPIDBase::GetPhiPrime(Double_t phi, Double_t magField, Int_t charge) const | |
193 | { | |
194 | // Get phiPrime which is the cut variable to reject high pT tracks near edges | |
195 | ||
196 | if(magField < 0) // for negatve polarity field | |
197 | phi = TMath::TwoPi() - phi; | |
198 | ||
199 | if(charge < 0) // for negatve charge | |
200 | phi = TMath::TwoPi() - phi; | |
201 | ||
202 | phi += TMath::Pi() / 18.0; // to center gap in the middle | |
203 | phi = fmod(phi, TMath::Pi() / 9.0); | |
204 | return phi; | |
205 | } | |
206 | ||
207 | ||
208 | //______________________________________________________________________________ | |
209 | Bool_t AliTPCPIDBase::PhiPrimeCut(Double_t trackPt, Double_t trackPhi, Short_t trackCharge, Double_t magField) const | |
210 | { | |
211 | // Apply phi' cut to given track parameters | |
212 | ||
213 | if (trackPt < 2.0) | |
214 | return kTRUE; | |
215 | ||
216 | Double_t phiPrime = GetPhiPrime(trackPhi, magField, trackCharge); | |
217 | ||
218 | if (phiPrime < fPhiCutHigh->Eval(trackPt) && phiPrime > fPhiCutLow->Eval(trackPt)) | |
219 | return kFALSE; // reject track | |
220 | ||
221 | return kTRUE; | |
222 | } | |
223 | ||
224 | ||
225 | //______________________________________________________________________________ | |
226 | Bool_t AliTPCPIDBase::PhiPrimeCut(const AliVTrack* track, Double_t magField) const | |
227 | { | |
228 | return PhiPrimeCut(track->Pt(), track->Phi(), track->Charge(), magField); | |
229 | } | |
230 | ||
231 | ||
232 | //______________________________________________________________________________ | |
233 | Bool_t AliTPCPIDBase::GetVertexIsOk(AliVEvent* event, Bool_t doVtxZcut) const | |
234 | { | |
235 | // Check whether vertex ful-fills quality requirements. | |
236 | // Apply cut on z-position of vertex if doVtxZcut = kTRUE. | |
237 | ||
238 | AliAODEvent* aod = 0x0; | |
239 | AliESDEvent* esd = 0x0; | |
240 | ||
241 | aod = dynamic_cast<AliAODEvent*>(event); | |
242 | if (!aod) { | |
243 | esd = dynamic_cast<AliESDEvent*>(event); | |
244 | ||
245 | if (!esd) { | |
246 | AliError("Event seems to be neither AOD nor ESD!"); | |
247 | return kFALSE; | |
248 | } | |
249 | } | |
250 | ||
251 | if (fIsPbpOrpPb) { | |
252 | const AliVVertex* trkVtx = (aod ? dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertex()) : | |
253 | dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertex())); | |
254 | ||
255 | if (!trkVtx || trkVtx->GetNContributors() <= 0) | |
256 | return kFALSE; | |
257 | ||
258 | TString vtxTtl = trkVtx->GetTitle(); | |
259 | if (!vtxTtl.Contains("VertexerTracks")) | |
260 | return kFALSE; | |
261 | ||
262 | Float_t zvtx = trkVtx->GetZ(); | |
263 | const AliVVertex* spdVtx = (aod ? dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertexSPD()) : | |
264 | dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertexSPD())); | |
265 | if (spdVtx->GetNContributors() <= 0) | |
266 | return kFALSE; | |
267 | ||
268 | TString vtxTyp = spdVtx->GetTitle(); | |
269 | Double_t cov[6] = {0}; | |
270 | spdVtx->GetCovarianceMatrix(cov); | |
271 | Double_t zRes = TMath::Sqrt(cov[5]); | |
272 | if (vtxTyp.Contains("vertexer:Z") && (zRes > 0.25)) | |
273 | return kFALSE; | |
274 | ||
275 | if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ()) > 0.5) | |
276 | return kFALSE; | |
277 | ||
278 | if (doVtxZcut) { | |
279 | if (TMath::Abs(zvtx) > fZvtxCutEvent) //Default: 10 cm | |
280 | return kFALSE; | |
281 | } | |
282 | ||
283 | return kTRUE; | |
284 | } | |
285 | ||
286 | ||
287 | // pp and PbPb | |
288 | const AliVVertex* primaryVertex = (aod ? dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertex()) : | |
289 | dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertexTracks())); | |
290 | ||
291 | if (!primaryVertex || primaryVertex->GetNContributors() <= 0) | |
292 | return kFALSE; | |
293 | ||
294 | if (doVtxZcut) { | |
295 | if (TMath::Abs(primaryVertex->GetZ()) > fZvtxCutEvent) //Default: 10 cm | |
296 | return kFALSE; | |
297 | } | |
298 | ||
299 | return kTRUE; | |
300 | } | |
301 | ||
302 | ||
303 | //______________________________________________________________________________ | |
304 | void AliTPCPIDBase::FillV0PIDlist(AliESDEvent* event) | |
305 | { | |
306 | // | |
307 | // Fill the PID tag list | |
308 | // | |
309 | ||
310 | // If no event forwarded as parameter (default), cast current input event. | |
311 | // Dynamic cast to ESD events (DO NOTHING for AOD events) | |
312 | if (!event) | |
313 | event = dynamic_cast<AliESDEvent *>(InputEvent()); | |
314 | ||
315 | // If this fails, do nothing | |
316 | if (!event) { | |
317 | AliError("Failed to retrieve ESD event. No V0's processed (only works for ESDs by now)."); | |
318 | return; | |
319 | } | |
320 | ||
321 | if (!fV0KineCuts) { | |
322 | AliError("V0KineCuts not available!"); | |
323 | return; | |
324 | } | |
325 | ||
326 | TString beamType(event->GetBeamType()); | |
327 | ||
328 | if (beamType.CompareTo("Pb-Pb") == 0 || beamType.CompareTo("A-A") == 0) { | |
329 | fV0KineCuts->SetMode(AliESDv0KineCuts::kPurity, AliESDv0KineCuts::kPbPb); | |
330 | } | |
331 | else { | |
332 | fV0KineCuts->SetMode(AliESDv0KineCuts::kPurity, AliESDv0KineCuts::kPP); | |
333 | } | |
334 | ||
335 | // V0 selection | |
336 | // set event | |
337 | fV0KineCuts->SetEvent(event); | |
338 | ||
339 | const Int_t numTracks = event->GetNumberOfTracks(); | |
340 | fV0tags = new Char_t[numTracks]; | |
341 | for (Int_t i = 0; i < numTracks; i++) | |
342 | fV0tags[i] = 0; | |
343 | ||
344 | if (fStoreMotherIndex) { | |
345 | fV0motherIndex = new Int_t[numTracks]; | |
346 | for (Int_t i = 0; i < numTracks; i++) | |
347 | fV0motherIndex[i] = -1; | |
348 | } | |
349 | ||
350 | fNumTagsStored = numTracks; | |
351 | ||
352 | // loop over V0 particles | |
353 | for (Int_t iv0 = 0; iv0 < event->GetNumberOfV0s(); iv0++) { | |
354 | AliESDv0* v0 = (AliESDv0*)event->GetV0(iv0); | |
355 | ||
356 | if (!v0) | |
357 | continue; | |
358 | ||
359 | // Reject onFly V0's <-> Only take V0's from offline V0 finder | |
360 | if (v0->GetOnFlyStatus()) | |
361 | continue; | |
362 | ||
363 | // Get the particle selection | |
364 | Bool_t foundV0 = kFALSE; | |
365 | Int_t pdgV0 = 0, pdgP = 0, pdgN = 0; | |
366 | ||
367 | foundV0 = fV0KineCuts->ProcessV0(v0, pdgV0, pdgP, pdgN); | |
368 | if (!foundV0) | |
369 | continue; | |
370 | ||
371 | Int_t iTrackP = v0->GetPindex(); // positive track | |
372 | Int_t iTrackN = v0->GetNindex(); // negative track | |
373 | ||
374 | ||
375 | // Fill the Object arrays | |
376 | // positive particles | |
377 | if (pdgP == -11) { | |
378 | fV0tags[iTrackP] = 14; | |
379 | } | |
380 | else if (pdgP == 211) { | |
381 | fV0tags[iTrackP] = 15; | |
382 | } | |
383 | else if(pdgP == 2212) { | |
384 | fV0tags[iTrackP] = 16; | |
385 | } | |
386 | ||
387 | if (fStoreMotherIndex) | |
388 | fV0motherIndex[iTrackP] = iv0; | |
389 | ||
390 | // negative particles | |
391 | if( pdgN == 11){ | |
392 | fV0tags[iTrackN] = -14; | |
393 | } | |
394 | else if( pdgN == -211){ | |
395 | fV0tags[iTrackN] = -15; | |
396 | } | |
397 | else if( pdgN == -2212){ | |
398 | fV0tags[iTrackN] = -16; | |
399 | } | |
400 | ||
401 | if (fStoreMotherIndex) | |
402 | fV0motherIndex[iTrackN] = iv0; | |
403 | } | |
404 | } | |
405 | ||
406 | ||
407 | //______________________________________________________________________________ | |
408 | void AliTPCPIDBase::ClearV0PIDlist() | |
409 | { | |
410 | // | |
411 | // Clear the PID tag list | |
412 | // | |
413 | ||
414 | delete fV0tags; | |
415 | fV0tags = 0; | |
416 | ||
417 | delete fV0motherIndex; | |
418 | fV0motherIndex = 0; | |
419 | ||
420 | fNumTagsStored = 0; | |
421 | } | |
422 | ||
423 | ||
424 | //______________________________________________________________________________ | |
425 | Char_t AliTPCPIDBase::GetV0tag(Int_t trackIndex) const | |
426 | { | |
427 | // | |
428 | // Get the tag for the corresponding trackIndex. Returns -99 in case of invalid index/tag list. | |
429 | // | |
430 | ||
431 | if (trackIndex < 0 || trackIndex >= fNumTagsStored || !fV0tags) | |
432 | return -99; | |
433 | ||
434 | return fV0tags[trackIndex]; | |
435 | } | |
436 | ||
437 | ||
438 | //______________________________________________________________________________ | |
439 | Int_t AliTPCPIDBase::GetV0motherIndex(Int_t trackIndex) const | |
440 | { | |
441 | // | |
442 | // Get the index of the V0 mother for the corresponding trackIndex. Returns -99 in case of invalid index/mother index list. | |
443 | // | |
444 | ||
445 | if (!fStoreMotherIndex || trackIndex < 0 || trackIndex >= fNumTagsStored || !fV0motherIndex) | |
446 | return -99; | |
447 | ||
448 | return fV0motherIndex[trackIndex]; | |
449 | } | |
450 | ||
451 | ||
452 | //________________________________________________________________________ | |
453 | Bool_t AliTPCPIDBase::TPCCutMIGeo(const AliVTrack* track, const AliVEvent* evt, TTreeStream* streamer) | |
454 | { | |
455 | // | |
456 | // TPC Cut MIGeo | |
457 | // | |
458 | ||
459 | if (!track || !evt) | |
460 | return kFALSE; | |
461 | ||
462 | const Short_t sign = track->Charge(); | |
463 | Double_t xyz[50]; | |
464 | Double_t pxpypz[50]; | |
465 | Double_t cv[100]; | |
466 | ||
467 | track->GetXYZ(xyz); | |
468 | track->GetPxPyPz(pxpypz); | |
469 | ||
470 | AliExternalTrackParam* par = new AliExternalTrackParam(xyz, pxpypz, cv, sign); | |
471 | const AliESDtrack dummy; | |
472 | ||
473 | const Double_t magField = evt->GetMagneticField(); | |
474 | Double_t varGeom = dummy.GetLengthInActiveZone(par, 3, 236, magField, 0, 0); | |
475 | Double_t varNcr = track->GetTPCClusterInfo(3, 1); | |
476 | Double_t varNcls = track->GetTPCsignalN(); | |
477 | ||
478 | const Double_t varEval = 130. - 5. * TMath::Abs(1. / track->Pt()); | |
479 | Bool_t cutGeom = varGeom > fgCutGeo * varEval; | |
480 | Bool_t cutNcr = varNcr > fgCutNcr * varEval; | |
481 | Bool_t cutNcls = varNcls > fgCutNcl * varEval; | |
482 | ||
483 | Bool_t kout = cutGeom && cutNcr && cutNcls; | |
484 | ||
485 | if (streamer) { | |
486 | Double_t dedx = track->GetTPCsignal(); | |
487 | ||
488 | (*streamer)<<"tree"<< | |
489 | "param.="<< par<< | |
490 | "varGeom="<<varGeom<< | |
491 | "varNcr="<<varNcr<< | |
492 | "varNcls="<<varNcls<< | |
493 | ||
494 | "cutGeom="<<cutGeom<< | |
495 | "cutNcr="<<cutNcr<< | |
496 | "cutNcls="<<cutNcls<< | |
497 | ||
498 | "kout="<<kout<< | |
499 | "dedx="<<dedx<< | |
500 | ||
501 | "\n"; | |
502 | } | |
503 | ||
504 | delete par; | |
505 | ||
506 | return kout; | |
507 | } | |
508 | ||
509 | ||
510 | //________________________________________________________________________ | |
511 | Bool_t AliTPCPIDBase::TPCnclCut(const AliVTrack* track) | |
512 | { | |
513 | // | |
514 | // TPC Cut on TPCsignalN() only | |
515 | // | |
516 | ||
517 | if (!track) | |
518 | return kFALSE; | |
519 | ||
520 | return (track->GetTPCsignalN() >= fgCutPureNcl); | |
521 | } |