AliPhysics  58f3d52 (58f3d52)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskPIDV0base.cxx
Go to the documentation of this file.
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 
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 
38 
39 Double_t AliAnalysisTaskPIDV0base::fgCutGeo = 1.;
40 Double_t AliAnalysisTaskPIDV0base::fgCutNcr = 0.85;
41 Double_t AliAnalysisTaskPIDV0base::fgCutNcl = 0.7;
42 
43 UShort_t AliAnalysisTaskPIDV0base::fgCutPureNcl = 60;
44 
45 //________________________________________________________________________
48  , fEvent(0x0)
49  , fESD(0x0)
50  , fMC(0x0)
51  , fPIDResponse(0x0)
52  , fV0KineCuts(0x0)
53  , fAnaUtils(0x0)
54  , fRunMode(kJetPIDMode)
55  , fPileUpRejectionType(kPileUpRejectionOff)
56  , fMinPlpContribSPD(3)
57  , fIsPbpOrpPb(kFALSE)
58  , fUsePhiCut(kFALSE)
59  , fTPCcutType(kNoCut)
60  , fZvtxCutEvent(10.0)
61  , fEtaCut(0.9)
62  , fPhiCutLow(0x0)
63  , fPhiCutHigh(0x0)
64  , fRandom(0x0)
65  , fTrackFilter(0x0)
66  , fNumTagsStored(0)
67  , fV0tags(0x0)
68  , fStoreMotherIndex(kFALSE)
69  , fV0motherIndex(0x0)
70 {
71  // default Constructor
72 
73  fRandom = new TRandom3(0); // 0 means random seed
74 
75 
76  // Question: Is this the right place to initialize these functions?
77  // Will it work on proof? i.e. will they be streamed to the workers?
78  // Also one should add getters and setters
79  fPhiCutLow = new TF1("StdPhiCutLow", "0.1/x/x+pi/18.0-0.025", 0, 100);
80  fPhiCutHigh = new TF1("StandardPhiCutHigh", "0.12/x+pi/18.0+0.035", 0, 100);
81 }
82 
83 //________________________________________________________________________
85  : AliAnalysisTaskSE(name)
86  , fEvent(0x0)
87  , fESD(0x0)
88  , fMC(0x0)
89  , fPIDResponse(0x0)
90  , fV0KineCuts(0x0)
91  , fAnaUtils(0x0)
92  , fRunMode(kJetPIDMode)
93  , fPileUpRejectionType(kPileUpRejectionOff)
94  , fMinPlpContribSPD(3)
95  , fIsPbpOrpPb(kFALSE)
96  , fUsePhiCut(kFALSE)
97  , fTPCcutType(kNoCut)
98  , fZvtxCutEvent(10.0)
99  , fEtaCut(0.9)
100  , fPhiCutLow(0x0)
101  , fPhiCutHigh(0x0)
102  , fRandom(0x0)
103  , fTrackFilter(0x0)
104  , fNumTagsStored(0)
105  , fV0tags(0x0)
106  , fStoreMotherIndex(kFALSE)
107  , fV0motherIndex(0x0)
108 {
109  // Constructor
110 
111  fRandom = new TRandom3(0); // 0 means random seed
112 
113 
114  // Question: Is this the right place to initialize these functions?
115  // Will it work on proof? i.e. will they be streamed to the workers?
116  // Also one should add getters and setters
117  fPhiCutLow = new TF1("StdPhiCutLow", "0.1/x/x+pi/18.0-0.025", 0, 100);
118  fPhiCutHigh = new TF1("StandardPhiCutHigh", "0.12/x+pi/18.0+0.035", 0, 100);
119 
120  DefineInput (0, TChain::Class());
121  DefineOutput(0, TTree::Class());
122 }
123 
124 
125 //________________________________________________________________________
127 {
128  // dtor
129 
130  delete fPhiCutLow;
131  fPhiCutLow = 0;
132 
133  delete fPhiCutHigh;
134  fPhiCutHigh = 0;
135 
136  delete fTrackFilter;
137  fTrackFilter = 0;
138 
139  delete fRandom;
140  fRandom = 0;
141 
142  delete fV0KineCuts;
143  fV0KineCuts = 0;
144 
145  delete [] fV0tags;
146  fV0tags = 0;
147  fNumTagsStored = 0;
148 
149  delete [] fV0motherIndex;
150  fV0motherIndex = 0;
151 
152  delete fAnaUtils;
153  fAnaUtils = 0;
154 }
155 
156 
157 //________________________________________________________________________
159 {
160  // Create histograms
161  // Called once
162 
163  // Input hander
164  AliAnalysisManager* man = AliAnalysisManager::GetAnalysisManager();
165  AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
166 
167  if (!inputHandler) {
168  AliFatal("Input handler needed");
169  fPIDResponse = 0x0;
170 
171  return;
172  }
173 
174  // PID response object
175  fPIDResponse = inputHandler->GetPIDResponse();
176  if (!fPIDResponse)
177  AliError("PIDResponse object was not created");
178 
179  // V0 Kine cuts
180  fV0KineCuts = new AliESDv0KineCuts;
181  fV0KineCuts->SetGammaCutChi2NDF(5.);
182  // Only accept V0el with prod. radius within 45 cm -> PID will by systematically biased for larger values!
183  Float_t gammaProdVertexRadiusCuts[2] = { 3.0, 45. };
184  fV0KineCuts->SetGammaCutVertexR(&gammaProdVertexRadiusCuts[0]);
185 
186  // Default analysis utils
187  fAnaUtils = new AliAnalysisUtils();
188 
189  //Sets the minimum number of contributors to the SPD vertext (Pile-Up)
190  fAnaUtils->SetMinPlpContribSPD(GetMinPlpContribSPD());
191 
192  // Not used yet, but to be save, forward vertex z cut to analysis utils object
193  fAnaUtils->SetMaxVtxZ(fZvtxCutEvent);
194 }
195 
196 
197 //________________________________________________________________________
199 {
200  // Main loop
201  // Called for each event
202 }
203 
204 //________________________________________________________________________
206 {
207  // Called once at the end of the query
208 }
209 
210 
211 //_____________________________________________________________________________
213 {
214  // Get phiPrime which is the cut variable to reject high pT tracks near edges
215 
216  if(magField < 0) // for negatve polarity field
217  phi = TMath::TwoPi() - phi;
218 
219  if(charge < 0) // for negatve charge
220  phi = TMath::TwoPi() - phi;
221 
222  phi += TMath::Pi() / 18.0; // to center gap in the middle
223  phi = fmod(phi, TMath::Pi() / 9.0);
224  return phi;
225 }
226 
227 
228 //______________________________________________________________________________
229 Bool_t AliAnalysisTaskPIDV0base::PhiPrimeCut(Double_t trackPt, Double_t trackPhi, Short_t trackCharge, Double_t magField) const
230 {
231  // Apply phi' cut to given track parameters
232 
233  if (trackPt < 2.0)
234  return kTRUE;
235 
236  Double_t phiPrime = GetPhiPrime(trackPhi, magField, trackCharge);
237 
238  if (phiPrime < fPhiCutHigh->Eval(trackPt) && phiPrime > fPhiCutLow->Eval(trackPt))
239  return kFALSE; // reject track
240 
241  return kTRUE;
242 }
243 
244 
245 //______________________________________________________________________________
246 Bool_t AliAnalysisTaskPIDV0base::PhiPrimeCut(const AliVTrack* track, Double_t magField) const
247 {
248  return PhiPrimeCut(track->Pt(), track->Phi(), track->Charge(), magField);
249 }
250 
251 
252 //______________________________________________________________________________
253 Bool_t AliAnalysisTaskPIDV0base::GetVertexIsOk(AliVEvent* event, Bool_t doVtxZcut) const
254 {
255  // Check whether vertex fulfills quality requirements.
256  // Apply cut on z-position of vertex if doVtxZcut = kTRUE.
257 
258  AliAODEvent* aod = 0x0;
259  AliESDEvent* esd = 0x0;
260 
261  aod = dynamic_cast<AliAODEvent*>(event);
262  if (!aod) {
263  esd = dynamic_cast<AliESDEvent*>(event);
264 
265  if (!esd) {
266  AliError("Event seems to be neither AOD nor ESD!");
267  return kFALSE;
268  }
269  }
270 
271  if (GetRunMode() == kJetPIDMode) {
272 
273  if (fIsPbpOrpPb) {
274  const AliVVertex* trkVtx = event->GetPrimaryVertex();
275 
276  if (!trkVtx || !(trkVtx->GetStatus()))
277  return kFALSE;
278 
279  TString vtxTtl = trkVtx->GetTitle();
280  if (!vtxTtl.Contains("VertexerTracks"))
281  return kFALSE;
282 
283  Float_t zvtx = trkVtx->GetZ();
284  const AliVVertex* spdVtx = event->GetPrimaryVertexSPD();
285 
286  if (!spdVtx || !(spdVtx->GetStatus()))
287  return kFALSE;
288 
289  TString vtxTyp = spdVtx->GetTitle();
290  Double_t cov[6] = {0};
291  spdVtx->GetCovarianceMatrix(cov);
292  Double_t zRes = TMath::Sqrt(cov[5]);
293  if (vtxTyp.Contains("vertexer:Z") && (zRes > 0.25))
294  return kFALSE;
295 
296  if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ()) > 0.5)
297  return kFALSE;
298 
299  if (doVtxZcut) {
300  if (TMath::Abs(zvtx) > fZvtxCutEvent) //Default: 10 cm
301  return kFALSE;
302  }
303 
304  return kTRUE;
305  }
306  else {
307  // pp and PbPb
308  const AliVVertex* primaryVertex = 0x0;
309  if (aod) {
310  primaryVertex = dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertex());
311  if (!primaryVertex || primaryVertex->GetNContributors() <= 0)
312  return kFALSE;
313 
314  // Reject TPC vertices
315  TString primVtxTitle(primaryVertex->GetTitle());
316  if (primVtxTitle.Contains("TPCVertex",TString::kIgnoreCase))
317  return kFALSE;
318  }
319  else {
320  primaryVertex = dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertexTracks());
321  if (!primaryVertex)
322  return kFALSE;
323 
324  if (primaryVertex->GetNContributors() <= 0) {
325  // Try SPD vertex
326  primaryVertex = dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertexSPD());
327  if (!primaryVertex || primaryVertex->GetNContributors() <= 0)
328  return kFALSE;
329  }
330  }
331 
332  if (doVtxZcut) {
333  if (TMath::Abs(primaryVertex->GetZ()) > fZvtxCutEvent) //Default: 10 cm
334  return kFALSE;
335  }
336  }
337  }
338 
339  if (GetRunMode() == kLightFlavorMode) {
340  const AliESDVertex *trkVertex = esd->GetPrimaryVertexTracks();
341  const AliESDVertex *spdVertex = esd->GetPrimaryVertexSPD();
342 
343  Bool_t hasSPD = spdVertex->GetStatus();
344  Bool_t hasTrk = trkVertex->GetStatus();
345 
346  if (!hasSPD)
347  return kFALSE;
348 
349  if (spdVertex->IsFromVertexerZ() && !(spdVertex->GetDispersion()<0.04 && spdVertex->GetZRes()<0.25))
350  return kFALSE;
351 
352  if (TMath::Abs(spdVertex->GetZ() - trkVertex->GetZ())>0.5)
353  return kFALSE;
354 
355  if (doVtxZcut) {
356  const AliVVertex *vertex = event->GetPrimaryVertex();
357  if (TMath::Abs(vertex->GetZ())>GetZvtxCutEvent())
358  return kFALSE;
359  }
360  }
361 
362  return kTRUE;
363 }
364 
365 
366 //______________________________________________________________________________
367 Bool_t AliAnalysisTaskPIDV0base::GetIsPileUp(AliVEvent* event, PileUpRejectionType pileUpRejection) const
368 {
369  // Check whether event is a pile-up event according to current AnalysisUtils object.
370  // If rejection type is kPileUpRejectionOff, kFALSE is returned.
371  // In case of errors, the error is displayed and kTRUE is returned.
372 
373  PileUpRejectionType functionPURejectionType = GetPileUpRejectionType();
374 
375  if (!(pileUpRejection == kPileUpRejectionClass))
376  functionPURejectionType == pileUpRejection;
377 
378  if (functionPURejectionType == kPileUpRejectionOff)
379  return kFALSE;
380 
381  if (!event) {
382  AliError("No event!");
383  return kTRUE;
384  }
385 
386  if (!fAnaUtils) {
387  AliError("AnalysisUtils object not available!");
388  return kTRUE;
389  }
390 
391  if (functionPURejectionType == kPileUpRejectionSPD)
392  return fAnaUtils->IsPileUpSPD(event);
393  else if (functionPURejectionType == kPileUpRejectionMV)
394  return fAnaUtils->IsPileUpMV(event);
395 
396  return kTRUE;
397 }
398 
399 
400 //______________________________________________________________________________
402 {
403  //
404  // Fill the PID tag list
405  //
406 
407  // If no event forwarded as parameter (default), cast current input event.
408  // Dynamic cast to ESD events (DO NOTHING for AOD events)
409  if (!event)
410  event = dynamic_cast<AliESDEvent *>(InputEvent());
411 
412  // If this fails, do nothing
413  if (!event) {
414  AliError("Failed to retrieve ESD event. No V0's processed (only works for ESDs by now).");
415  return;
416  }
417 
418  if (!fV0KineCuts) {
419  AliError("V0KineCuts not available!");
420  return;
421  }
422 
423  TString beamType(event->GetBeamType());
424 
425  if (beamType.CompareTo("Pb-Pb") == 0 || beamType.CompareTo("A-A") == 0) {
426  fV0KineCuts->SetMode(AliESDv0KineCuts::kPurity, AliESDv0KineCuts::kPbPb);
427  }
428  else {
429  fV0KineCuts->SetMode(AliESDv0KineCuts::kPurity, AliESDv0KineCuts::kPP);
430  }
431 
432  // V0 selection
433  // set event
434  fV0KineCuts->SetEvent(event);
435 
436  const Int_t numTracks = event->GetNumberOfTracks();
437  fV0tags = new Char_t[numTracks];
438  for (Int_t i = 0; i < numTracks; i++)
439  fV0tags[i] = 0;
440 
441  if (fStoreMotherIndex) {
442  fV0motherIndex = new Int_t[numTracks];
443  for (Int_t i = 0; i < numTracks; i++)
444  fV0motherIndex[i] = -1;
445  }
446 
447  fNumTagsStored = numTracks;
448 
449  // loop over V0 particles
450  for (Int_t iv0 = 0; iv0 < event->GetNumberOfV0s(); iv0++) {
451  AliESDv0* v0 = (AliESDv0*)event->GetV0(iv0);
452 
453  if (!v0)
454  continue;
455 
456  // Reject onFly V0's <-> Only take V0's from offline V0 finder
457  if (v0->GetOnFlyStatus())
458  continue;
459 
460  // Get the particle selection
461  Bool_t foundV0 = kFALSE;
462  Int_t pdgV0 = 0, pdgP = 0, pdgN = 0;
463 
464  foundV0 = fV0KineCuts->ProcessV0(v0, pdgV0, pdgP, pdgN);
465  if (!foundV0)
466  continue;
467 
468  Int_t iTrackP = v0->GetPindex(); // positive track
469  Int_t iTrackN = v0->GetNindex(); // negative track
470 
471 
472  // Fill the Object arrays
473  // positive particles
474  if (pdgP == -11) {
475  fV0tags[iTrackP] = 14;
476  }
477  else if (pdgP == 211) {
478  fV0tags[iTrackP] = 15;
479  }
480  else if(pdgP == 2212) {
481  fV0tags[iTrackP] = 16;
482  }
483 
484  if (fStoreMotherIndex)
485  fV0motherIndex[iTrackP] = iv0;
486 
487  // negative particles
488  if( pdgN == 11){
489  fV0tags[iTrackN] = -14;
490  }
491  else if( pdgN == -211){
492  fV0tags[iTrackN] = -15;
493  }
494  else if( pdgN == -2212){
495  fV0tags[iTrackN] = -16;
496  }
497 
498  if (fStoreMotherIndex)
499  fV0motherIndex[iTrackN] = iv0;
500  }
501 }
502 
503 
504 //______________________________________________________________________________
506 {
507  //
508  // Clear the PID tag list
509  //
510 
511  delete [] fV0tags;
512  fV0tags = 0;
513 
514  delete [] fV0motherIndex;
515  fV0motherIndex = 0;
516 
517  fNumTagsStored = 0;
518 }
519 
520 
521 //______________________________________________________________________________
523 {
524  //
525  // Get the tag for the corresponding trackIndex. Returns -99 in case of invalid index/tag list.
526  //
527 
528  if (trackIndex < 0 || trackIndex >= fNumTagsStored || !fV0tags)
529  return -99;
530 
531  return fV0tags[trackIndex];
532 }
533 
534 
535 //______________________________________________________________________________
537 {
538  //
539  // Get the index of the V0 mother for the corresponding trackIndex. Returns -99 in case of invalid index/mother index list.
540  //
541 
542  if (!fStoreMotherIndex || trackIndex < 0 || trackIndex >= fNumTagsStored || !fV0motherIndex)
543  return -99;
544 
545  return fV0motherIndex[trackIndex];
546 }
547 
548 
549 //________________________________________________________________________
550 Bool_t AliAnalysisTaskPIDV0base::TPCCutMIGeo(const AliVTrack* track, const AliVEvent* evt, TTreeStream* streamer)
551 {
552  //
553  // TPC Cut MIGeo
554  //
555 
556  if (!track || !evt)
557  return kFALSE;
558 
559  const Short_t sign = track->Charge();
560  Double_t xyz[50];
561  Double_t pxpypz[50];
562  Double_t cv[100];
563 
564  track->GetXYZ(xyz);
565  track->GetPxPyPz(pxpypz);
566 
567  AliExternalTrackParam* par = new AliExternalTrackParam(xyz, pxpypz, cv, sign);
568  const AliESDtrack dummy;
569 
570  const Double_t magField = evt->GetMagneticField();
571  Double_t varGeom = dummy.GetLengthInActiveZone(par, 3, 236, magField, 0, 0);
572  Double_t varNcr = track->GetTPCClusterInfo(3, 1);
573  Double_t varNcls = track->GetTPCsignalN();
574 
575  const Double_t varEval = 130. - 5. * TMath::Abs(1. / track->Pt());
576  Bool_t cutGeom = varGeom > fgCutGeo * varEval;
577  Bool_t cutNcr = varNcr > fgCutNcr * varEval;
578  Bool_t cutNcls = varNcls > fgCutNcl * varEval;
579 
580  Bool_t kout = cutGeom && cutNcr && cutNcls;
581 
582  if (streamer) {
583  Double_t dedx = track->GetTPCsignal();
584 
585  (*streamer)<<"tree"<<
586  "param.="<< par<<
587  "varGeom="<<varGeom<<
588  "varNcr="<<varNcr<<
589  "varNcls="<<varNcls<<
590 
591  "cutGeom="<<cutGeom<<
592  "cutNcr="<<cutNcr<<
593  "cutNcls="<<cutNcls<<
594 
595  "kout="<<kout<<
596  "dedx="<<dedx<<
597 
598  "\n";
599  }
600 
601  delete par;
602 
603  return kout;
604 }
605 
606 
607 //________________________________________________________________________
609 {
610  //
611  // TPC Cut on TPCsignalN() only
612  //
613 
614  if (!track)
615  return kFALSE;
616 
617  return (track->GetTPCsignalN() >= fgCutPureNcl);
618 }
Int_t charge
virtual void Terminate(const Option_t *)
virtual Double_t GetZvtxCutEvent() const
virtual Int_t GetMinPlpContribSPD() const
double Double_t
Definition: External.C:58
AliAnalysisFilter * fTrackFilter
Can be used to statistically determine the shape in the pt bins e.g.
PileUpRejectionType GetPileUpRejectionType() const
static Bool_t TPCCutMIGeo(const AliVTrack *track, const AliVEvent *evt, TTreeStream *streamer=0x0)
char Char_t
Definition: External.C:18
Bool_t fStoreMotherIndex
Pointer to array with tags for identified particles from V0 decays.
virtual Double_t GetPhiPrime(Double_t phi, Double_t magField, Int_t charge) const
virtual Bool_t GetIsPileUp(AliVEvent *event, PileUpRejectionType pileUpRejection=kPileUpRejectionClass) const
virtual Bool_t PhiPrimeCut(const AliVTrack *track, Double_t magField) const
AliPIDResponse * fPIDResponse
MC object.
virtual Int_t GetV0motherIndex(Int_t trackIndex) const
AliESDv0KineCuts * fV0KineCuts
PID response Handler.
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
AliAnalysisUtils * fAnaUtils
ESD V0 kine cuts.
virtual void UserExec(Option_t *option)
static Bool_t TPCnclCut(const AliVTrack *track)
short Short_t
Definition: External.C:23
virtual Char_t GetV0tag(Int_t trackIndex) const
virtual Bool_t GetVertexIsOk(AliVEvent *event, Bool_t doVtxZcut=kTRUE) const
unsigned short UShort_t
Definition: External.C:28
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
void FillV0PIDlist(AliESDEvent *esdEvent=0x0)