PWGPP/TPC/AliTPCPIDBase.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 "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 }