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