don't lie in the log!
[u/mrichter/AliRoot.git] / PWGPP / TPC / AliTPCPIDBase.cxx
CommitLineData
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/*
28This class is a base class for all other
29analysis tasks that use V0's.
30It provides basics for V0 identification.
31In addition, some further basic functions are provided.
32
33Class written by Benjamin Hess.
34Contact: bhess@cern.ch
35*/
36
37ClassImp(AliTPCPIDBase)
38
39Double_t AliTPCPIDBase::fgCutGeo = 1.;
40Double_t AliTPCPIDBase::fgCutNcr = 0.85;
41Double_t AliTPCPIDBase::fgCutNcl = 0.7;
42
43UShort_t AliTPCPIDBase::fgCutPureNcl = 60;
44
45//________________________________________________________________________
46AliTPCPIDBase::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//________________________________________________________________________
80AliTPCPIDBase::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//________________________________________________________________________
118AliTPCPIDBase::~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//________________________________________________________________________
147void 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//________________________________________________________________________
178void AliTPCPIDBase::UserExec(Option_t *)
179{
180 // Main loop
181 // Called for each event
182}
183
184//________________________________________________________________________
185void AliTPCPIDBase::Terminate(const Option_t *)
186{
187 // Called once at the end of the query
188}
189
190
191//_____________________________________________________________________________
192Double_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//______________________________________________________________________________
209Bool_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//______________________________________________________________________________
226Bool_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//______________________________________________________________________________
233Bool_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//______________________________________________________________________________
304void 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//______________________________________________________________________________
408void 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//______________________________________________________________________________
425Char_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//______________________________________________________________________________
439Int_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//________________________________________________________________________
453Bool_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//________________________________________________________________________
511Bool_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}