PWGJE/UserTasks/AliAnalysisTaskPIDV0base.cxx

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