PWGPP/TPC/AliTPCPIDEtaTree.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TF1.h"
   4 #include "TAxis.h"
   5 #include "TH2I.h"
   6
   7 #include "THnSparse.h"
   8
   9 #include "AliMCParticle.h"
  10 //#include "AliStack.h"
  11
  12 #include "AliAnalysisTask.h"
  13 #include "AliAnalysisManager.h"
  14
  15 #include "AliESDEvent.h"
  16 #include "AliMCEvent.h"
  17 #include "AliESDInputHandler.h"
  18 #include "AliInputEventHandler.h"
  19
  20 #include "AliVVertex.h"
  21 #include "AliAnalysisFilter.h"
  22 #include "AliPID.h"
  23 #include "AliPIDResponse.h"
  24 #include "AliTPCPIDResponse.h"
  25 //#include "AliTPCParamSR.h"
  26
  27 #include "AliTPCPIDEtaTree.h"
  28
  29 /*
  30 This task determines the eta dependence of the TPC signal.
  31 For this purpose, only tracks fulfilling some standard quality cuts are taken into account.
  32 The obtained data can be used to derive the functional behaviour of the eta dependence.
  33 Such a function can be plugged into this task to correct for the eta dependence and to see
  34 if there is then really no eta dependence left.
  35
  36 Class written by Benjamin Hess.
  37 Contact: bhess@cern.ch
  38 */
  39
  40 ClassImp(AliTPCPIDEtaTree)
  41
  42 //________________________________________________________________________
  43 AliTPCPIDEtaTree::AliTPCPIDEtaTree()
  44   : AliTPCPIDBase()
  45   , fStoreMultiplicity(kFALSE)
  46   , fStoreNumOfSubthresholdclusters(kFALSE)
  47   , fStoreNumClustersInActiveVolume(kFALSE)
  48   , fDoAdditionalQA(kFALSE)
  49   , fPtpcPionCut(1.0)
  50   , fPtpc(0)
  51   , fPt(0)
  52   , fDeDx(0)
  53   , fDeDxExpected(0)
  54   , fTanTheta(0)
  55   //, fSinAlpha(0)
  56   //, fY(0)
  57   , fPhiPrime(0)
  58   , fTPCsignalN(0)
  59   , fTPCsignalNsubthreshold(0)
  60   , fNumTPCClustersInActiveVolume(0)
  61   , fPIDtype(0)
  62   , fMultiplicity(0)
  63   , fCorrectdEdxEtaDependence(kFALSE)
  64   , fCorrectdEdxMultiplicityDependence(kFALSE)
  65   , fTree(0x0)
  66   , fTreePions(0x0)
  67   , fOutputContainer(0x0)
  68   , fhTOFqa(0x0)
  69   , fhMultiplicityQA(0x0)
  70 {
  71   // default Constructor
  72 }
  73
  74 //________________________________________________________________________
  75 AliTPCPIDEtaTree::AliTPCPIDEtaTree(const char *name)
  76   : AliTPCPIDBase(name)
  77   , fStoreMultiplicity(kFALSE)
  78   , fStoreNumOfSubthresholdclusters(kFALSE)
  79   , fStoreNumClustersInActiveVolume(kFALSE)
  80   , fDoAdditionalQA(kFALSE)
  81   , fPtpcPionCut(1.0)
  82   , fPtpc(0)
  83   , fPt(0)
  84   , fDeDx(0)
  85   , fDeDxExpected(0)
  86   , fTanTheta(0)
  87   //, fSinAlpha(0)
  88   //, fY(0)
  89   , fPhiPrime(0)
  90   , fTPCsignalN(0)
  91   , fTPCsignalNsubthreshold(0)
  92   , fNumTPCClustersInActiveVolume(0)
  93   , fPIDtype(0)
  94   , fMultiplicity(0)
  95   , fCorrectdEdxEtaDependence(kFALSE)
  96   , fCorrectdEdxMultiplicityDependence(kFALSE)
  97   , fTree(0x0)
  98   , fTreePions(0x0)
  99   , fOutputContainer(0x0)
 100   , fhTOFqa(0x0)
 101   , fhMultiplicityQA(0x0)
 102 {
 103   // Constructor
 104
 105   // Define input and output slots here
 106   DefineInput(0, TChain::Class());
 107   DefineOutput(1, TTree::Class());
 108   DefineOutput(2, TTree::Class());
 109   DefineOutput(3, TObjArray::Class());
 110 }
 111
 112
 113 //________________________________________________________________________
 114 AliTPCPIDEtaTree::~AliTPCPIDEtaTree()
 115 {
 116   // dtor
 117
 118   //delete fTree;
 119   //fTree = 0x0;
 120
 121   delete fOutputContainer;
 122   fOutputContainer = 0x0;
 123 }
 124
 125
 126 //________________________________________________________________________
 127 void AliTPCPIDEtaTree::UserCreateOutputObjects()
 128 {
 129   // Create histograms
 130   // Called once
 131
 132   AliTPCPIDBase::UserCreateOutputObjects();
 133
 134   if (fDoAdditionalQA) {
 135     OpenFile(3);
 136
 137     fOutputContainer = new TObjArray(2);
 138     fOutputContainer->SetName(GetName());
 139     fOutputContainer->SetOwner(kTRUE);
 140
 141     const Int_t nBins = 6;
 142     // p_vert, p_TPC, eta, nSigmaTOF, beta, multiplicity
 143     Int_t bins[nBins] =    {   48,  48,    9,   30, 110,    4  };
 144     Double_t xmin[nBins] = {  0.3, 0.3, -0.9, -5.0, 0.0,   0.  };
 145     Double_t xmax[nBins] = {  2.7, 2.7,  0.9,  5.0, 1.1, 3200  };
 146
 147
 148     fhTOFqa = new THnSparseI("hTOFqa","", nBins, bins, xmin, xmax);
 149     fhTOFqa->GetAxis(0)->SetTitle("p_{Vertex} (GeV/c)");
 150     fhTOFqa->GetAxis(1)->SetTitle("p_{TPC_inner} (GeV/c)");
 151     fhTOFqa->GetAxis(2)->SetTitle("#eta");
 152     fhTOFqa->GetAxis(3)->SetTitle("n #sigma_{p} TOF");
 153     fhTOFqa->GetAxis(4)->SetTitle("#beta");
 154     fhTOFqa->GetAxis(5)->SetTitle("Multiplicity");
 155
 156     fOutputContainer->Add(fhTOFqa);
 157
 158
 159
 160     fhMultiplicityQA = new TH2I("hMultiplicityQA", "Multiplicity check; Contributors to primary vertex per event; Total number of tracks per Event",
 161                                 120, 0, 6000, 400, 0, 20000);
 162     fOutputContainer->Add(fhMultiplicityQA);
 163   }
 164   else {
 165     fOutputContainer = new TObjArray(1);
 166     fOutputContainer->SetName(GetName());
 167     fOutputContainer->SetOwner(kTRUE);
 168   }
 169
 170   OpenFile(1);
 171
 172   fTree = new TTree("fTree", "Tree for analysis of #eta dependence of TPC signal");
 173   fTree->Branch("pTPC", &fPtpc);
 174   //fTree->Branch("pT", &fPt);
 175   fTree->Branch("dEdx", &fDeDx);
 176   fTree->Branch("dEdxExpected", &fDeDxExpected);
 177   fTree->Branch("tanTheta", &fTanTheta);
 178   //fTree->Branch("sinAlpha", &fSinAlpha);
 179   //fTree->Branch("y", &fY);
 180   //TODO fTree->Branch("phiPrime", &fPhiPrime);
 181   fTree->Branch("tpcSignalN", &fTPCsignalN);
 182
 183   if (fStoreNumOfSubthresholdclusters)
 184     fTree->Branch("tpcSignalNsubthreshold", &fTPCsignalNsubthreshold);
 185
 186   if (fStoreNumClustersInActiveVolume)
 187     fTree->Branch("numTPCClustersInActiveVolume", &fNumTPCClustersInActiveVolume);
 188
 189   fTree->Branch("pidType", &fPIDtype);
 190
 191   if (fStoreMultiplicity)  {
 192     fTree->Branch("fMultiplicity", &fMultiplicity);
 193   }
 194
 195   OpenFile(2);
 196
 197   fTreePions = new TTree("fTreePions", "Tree for analysis of #eta dependence of TPC signal - Pions");
 198   fTreePions->Branch("pTPC", &fPtpc);
 199   fTreePions->Branch("pT", &fPt);
 200   fTreePions->Branch("dEdx", &fDeDx);
 201   fTreePions->Branch("dEdxExpected", &fDeDxExpected);
 202   fTreePions->Branch("tanTheta", &fTanTheta);
 203   fTreePions->Branch("tpcSignalN", &fTPCsignalN);
 204
 205   if (fStoreNumOfSubthresholdclusters)
 206     fTreePions->Branch("tpcSignalNsubthreshold", &fTPCsignalNsubthreshold);
 207
 208   if (fStoreMultiplicity)  {
 209     fTreePions->Branch("fMultiplicity", &fMultiplicity);
 210   }
 211
 212   PostData(1, fTree);
 213   PostData(2, fTreePions);
 214   PostData(3, fOutputContainer);
 215 }
 216
 217
 218 //________________________________________________________________________
 219 void AliTPCPIDEtaTree::UserExec(Option_t *)
 220 {
 221   // Main loop
 222   // Called for each event
 223
 224   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 225   if (!fESD) {
 226     Printf("ERROR: fESD not available");
 227     return;
 228   }
 229
 230   fMC = dynamic_cast<AliMCEvent*>(MCEvent());
 231
 232   if (!fPIDResponse || !fV0KineCuts)
 233     return;
 234
 235   if (!GetVertexIsOk(fESD))
 236     return;
 237
 238   const AliVVertex* primaryVertex = fESD->GetPrimaryVertexTracks();
 239   if (!primaryVertex)
 240     return;
 241   //fMultiplicity = primaryVertex->GetNContributors();
 242
 243
 244   if(primaryVertex->GetNContributors() <= 0)
 245     return;
 246
 247   fMultiplicity = fESD->GetNumberOfTracks();
 248
 249   if (fDoAdditionalQA) {
 250     Int_t nTotTracks = fESD->GetNumberOfTracks();
 251     fhMultiplicityQA->Fill(primaryVertex->GetNContributors(), nTotTracks);
 252   }
 253
 254   Double_t magField = fESD->GetMagneticField();
 255
 256   // Fill V0 arrays with V0s
 257   FillV0PIDlist();
 258
 259   //AliTPCParamSR par;
 260   //par.Update();
 261
 262   // Track loop to fill a Train spectrum
 263   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 264     Bool_t isPr = kFALSE;
 265     Bool_t isPi = kFALSE;
 266
 267     AliESDtrack* track = fESD->GetTrack(iTracks);
 268     if (!track) {
 269       Printf("ERROR: Could not receive track %d", iTracks);
 270       continue;
 271     }
 272
 273     if (TMath::Abs(track->Eta()) > fEtaCut)
 274       continue;
 275
 276     // Do not distinguish positively and negatively charged V0's
 277     Char_t v0tagAllCharges = TMath::Abs(GetV0tag(iTracks));
 278     if (v0tagAllCharges == -99) {
 279       AliError(Form("Problem with V0 tag list (requested V0 track for track %d from %d (list status %d))!", iTracks, fESD->GetNumberOfTracks(),
 280                     fV0tags != 0x0));
 281       continue;
 282     }
 283
 284     Bool_t isV0prNotMC = (v0tagAllCharges == 16) && !fMC; // Only accept V0 protons for data, not for MC
 285     Bool_t isV0piNotMC = (v0tagAllCharges == 15) && !fMC; // Only accept V0 pions for data, not for MC
 286
 287     Bool_t isV0NotMC = isV0prNotMC || isV0piNotMC;
 288
 289     // Apply track cut
 290     if(!isV0NotMC && fTrackFilter && !fTrackFilter->IsSelected(track))
 291       continue;
 292
 293     // Note: For V0's, the cut on ncl can be done via the tree (value is stored). One should not cut on geometry
 294     // for V0's since they have different topology
 295     if (!isV0NotMC && GetUseTPCCutMIGeo()) {
 296       if (!TPCCutMIGeo(track, InputEvent()))
 297         continue;
 298     }
 299
 300
 301     fPtpc = track->GetTPCmomentum();
 302
 303     if (fPtpc > 5.)
 304       continue;
 305
 306     fPt = track->Pt();
 307     fPhiPrime = GetPhiPrime(track->Phi(), magField, track->Charge());
 308
 309     Int_t label = track->GetLabel();
 310
 311     AliMCParticle* mcTrack = 0x0;
 312
 313     if (fMC) {
 314       if (label < 0)
 315         continue; // If MC is available, reject tracks with negative ESD label
 316       mcTrack = dynamic_cast<AliMCParticle*>(fMC->GetTrack(TMath::Abs(label)));
 317       if (!mcTrack) {
 318         Printf("ERROR: Could not receive mcTrack with label %d for track %d", label, iTracks);
 319         continue;
 320       }
 321
 322       /*// Only accept MC primaries
 323       if (!fMC->Stack()->IsPhysicalPrimary(TMath::Abs(label))) {
 324         continue;
 325       }*/
 326     }
 327
 328     if (fMC) {
 329       Int_t pdg = mcTrack->PdgCode();
 330
 331       if (TMath::Abs(pdg) == 2212) // Proton
 332         isPr = kTRUE;
 333       else if ((pdg == 111 || TMath::Abs(pdg) == 211) && fPtpc <= fPtpcPionCut) // Pion below desired momentum threshold
 334         isPi = kTRUE;
 335       else
 336         continue; // Only take protons and pions
 337
 338       fPIDtype = kMCid;
 339       /*
 340       if (pdg == 111 || TMath::Abs(pdg) == 211) {//Pion
 341         binMC = 3;
 342       }
 343       else if (TMath::Abs(pdg) == 311 || TMath::Abs(pdg) == 321) {//Kaon
 344         binMC = 1;
 345       }
 346       else if (TMath::Abs(pdg) == 2212) {//Proton
 347         binMC = 4;
 348       }
 349       else if (TMath::Abs(pdg) == 11) {//Electron
 350         binMC = 0;
 351       }
 352       else if (TMath::Abs(pdg) == 13) {//Muon
 353         binMC = 2;
 354       }*/
 355     }
 356
 357     fDeDx = track->GetTPCsignal();
 358
 359     UInt_t status = track->GetStatus();
 360     Bool_t hasTOFout  = status&AliESDtrack::kTOFout;
 361     Bool_t hasTOFtime = status&AliESDtrack::kTIME;
 362     Bool_t hasTOF     = kFALSE;
 363     if (hasTOFout && hasTOFtime)
 364       hasTOF = kTRUE;
 365     Float_t length = track->GetIntegratedLength();
 366     // Length check only for primaries, not for V0's!
 367     if (length < 350 && !isV0NotMC)
 368       hasTOF = kFALSE;
 369
 370     if (!fMC) {
 371       // Note: Normally, the track cuts include a cut on primaries. Therefore, if the particle is a V0, it would usually
 372       // NOT be found via the primary cuts. This means that the PIDtype describes more or less disjoint sets
 373       if (isV0prNotMC) {
 374         // V0
 375         isPr = kTRUE;
 376
 377         if (!hasTOF) {
 378           fPIDtype = kV0idNoTOF;
 379         }
 380         else {
 381           if (TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton)) < 3.0)
 382             fPIDtype = kV0idPlusTOFaccepted;
 383           else
 384             fPIDtype = kV0idPlusTOFrejected;
 385         }
 386       }
 387       else if (isV0piNotMC) {
 388         // V0 pion
 389         if (fPtpc > fPtpcPionCut || fDeDx > 140.) // Reject pions above desired momentum threshold and also reject too high dEdx
 390           continue;
 391
 392         isPi = kTRUE;
 393
 394         if (!hasTOF) {
 395           fPIDtype = kV0idNoTOF;
 396         }
 397         else {
 398           if (TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion)) < 3.0)
 399             fPIDtype = kV0idPlusTOFaccepted;
 400           else
 401             fPIDtype = kV0idPlusTOFrejected;
 402         }
 403       }
 404       else {
 405         // Non-V0
 406         isPr = kTRUE; // If particle is accepted, it is a proton (for pions, only V0s are used)
 407
 408         if (fPtpc < 4.0 && //TODO was 2.7 // Do not accept non-V0s above this threshold -> High contamination!!!
 409             (fDeDx >= 50. / TMath::Power(fPtpc, 1.3))) {// Pattern recognition instead of TPC cut to be ~independent of old TPC expected dEdx
 410           if (fPtpc < 0.6) {
 411             fPIDtype = kTPCid;
 412           }
 413           // fPtpc >= 0.6
 414           else if (hasTOF && TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton)) < 3.0) {
 415               fPIDtype = kTPCandTOFid;
 416           }
 417           else
 418             continue; // Reject particle
 419         }
 420         else
 421           continue; // Reject particle
 422       }
 423     }
 424
 425     if (fDoAdditionalQA) {
 426       Double_t tofTime = track->GetTOFsignal();//in ps
 427       Double_t tof = tofTime * 1E-3; // ns, average T0 fill subtracted, no info from T0detector
 428       Double_t length2 = track->GetIntegratedLength();
 429       Double_t c = TMath::C() * 1.E-9;// m/ns
 430       length2 = length2 * 0.01; // in meters
 431       tof = tof * c;
 432       Double_t beta = length2 / tof;
 433
 434       Double_t nSigmaTOF = hasTOF ? fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton) : 999;
 435
 436       // p_vert, p_TPC, eta, nSigmaTOF, beta, multiplicity
 437       Double_t entry[6] = { track->P(), fPtpc, track->Eta(), nSigmaTOF, beta, (Double_t)fMultiplicity };
 438       fhTOFqa->Fill(entry);
 439     }
 440
 441     // Prepare entry for tree (some quantities have already been set)
 442     // Turn eta correction off -> Important here is the pure spline value and the selection via cuts. The correction
 443     // can be re-done manually, if needed. But it cannot be undone!
 444
 445     if (isPr)
 446       fDeDxExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
 447                                                                        fCorrectdEdxEtaDependence, fCorrectdEdxMultiplicityDependence);
 448     else if (isPi)
 449       fDeDxExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
 450                                                                        fCorrectdEdxEtaDependence, fCorrectdEdxMultiplicityDependence);
 451     else
 452       fDeDxExpected = -1.;
 453     fTanTheta = track->GetInnerParam()->GetTgl();
 454     //fSinAlpha = track->GetInnerParam()->GetSnp();
 455     //fY = track->GetInnerParam()->GetY();
 456     fTPCsignalN = track->GetTPCsignalN();
 457
 458     if (fStoreNumOfSubthresholdclusters) {
 459       const AliTPCdEdxInfo *clsInfo = track->GetTPCdEdxInfo();
 460
 461       if (clsInfo) {
 462         Double32_t signal[4] = {0}; // 0: IROC only; 1: OROC short padas; 2: OROC long pads; 3:OROC combined
 463         Char_t ncl[3] = {0};        // clusters above threshold; 0: IROC; 1: OROC short; 2: OROC long
 464         Char_t nrows[3] = {0};      // clusters above and below threshold; 0: IROC; 1: OROC short; 2: OROC long
 465
 466         clsInfo->GetTPCSignalRegionInfo(signal, ncl, nrows);
 467
 468         fTPCsignalNsubthreshold = nrows[0] + nrows[1] + nrows[2] - ncl[0] - ncl[1] - ncl[2];
 469       }
 470       else
 471         fTPCsignalNsubthreshold = 200;// Set to invalid value
 472     }
 473
 474
 475     if (fStoreNumClustersInActiveVolume)
 476       fNumTPCClustersInActiveVolume = track->GetLengthInActiveZone(1, 1.8, 220, magField);
 477     else
 478       fNumTPCClustersInActiveVolume = 200.;//Set to invalid value
 479
 480
 481
 482     /* START
 483     const Double_t tanTheta = track->GetInnerParam()->GetTgl();
 484
 485     // Constant in formula for B in kGauss (factor 0.1 to convert B from Tesla to kGauss),
 486     // pT in GeV/c (factor c*1e-9 to arrive at GeV/c) and curvature in 1/cm (factor 0.01 to get from m to cm)
 487     const Double_t constant = TMath::C()* 1e-9 * 0.1 * 0.01;
 488     const Double_t curvature = magField * constant / track->Pt(); // in 1./cm
 489
 490     //const Double_t angleIROC = TMath::ASin(TMath::Min(TMath::Abs(par.GetPadRowRadii(0,  par.GetNRow(0)  / 2.) * curvature) * 0.5, 1.));
 491     //const Double_t angleOROC = TMath::ASin(TMath::Min(TMath::Abs(par.GetPadRowRadii(36, par.GetNRow(36) / 2.) * curvature) * 0.5, 1.));
 492
 493     Double_t averageddzdr = 0.;
 494     Int_t nParts = 0;
 495
 496     for (Double_t r = 85; r < 245; r++) {
 497       Double_t sinPhiLocal = TMath::Abs(r*curvature*0.5);
 498
 499       // Cut on |sin(phi)| as during reco
 500       if (TMath::Abs(sinPhiLocal) <= 0.95) {
 501         const Double_t phiLocal = TMath::ASin(sinPhiLocal);
 502         const Double_t tanPhiLocal = TMath::Tan(phiLocal);
 503
 504         averageddzdr += tanTheta * TMath::Sqrt(1. + tanPhiLocal * tanPhiLocal);
 505         nParts++;
 506       }
 507     }
 508
 509     if (nParts > 0)
 510       averageddzdr /= nParts;
 511     else {
 512       AliError("Problems during determination of dz/dr. Skipping track!");
 513       continue;
 514     }
 515
 516
 517     //printf("padRow 0 / 36: %f / %f\n", par.GetPadRowRadii(0,  par.GetNRow(0) / 2.), par.GetPadRowRadii(36, par.GetNRow(36) / 2.));
 518     //printf("pT: %f\nFactor/magField(kGs)/curvature^-1: %f / %f /%f\nIROC/OROC/averageOld/average: %f / %f / %f / //%f\ntanThetaGlobalFromTheta/tanTheta/dzdr: %f / %f / %f\n\n",
 519     //        track->Pt(), constant, magField, 1./curvature, angleIROC, angleOROC, angleAverage, averageAngle, TMath::Tan(-track->Theta() + TMath::Pi() / 2.0), tanTheta, dzdr);
 520     //printf("pT: %f\nFactor/magField(kGs)/curvature^-1: %f / %f /%f\ntanThetaGlobalFromTheta/tanTheta/Averageddzdr: %f / %f / %f\n\n",
 521     //        track->Pt(), constant, magField, 1./curvature, TMath::Tan(-track->Theta() + TMath::Pi() / 2.0), tanTheta, averageddzdr);
 522
 523
 524     fTanTheta = averageddzdr;
 525     */
 526
 527     if (isPr)
 528       fTree->Fill();
 529     else if (isPi)
 530       fTreePions->Fill();
 531   } //track loop
 532
 533   // Post output data.
 534   PostData(1, fTree);
 535   PostData(2, fTreePions);
 536
 537   if (fDoAdditionalQA)
 538     PostData(3, fOutputContainer);
 539
 540   // Clear the V0 PID arrays
 541   ClearV0PIDlist();
 542 }
 543
 544 //________________________________________________________________________
 545 void AliTPCPIDEtaTree::Terminate(const Option_t *)
 546 {
 547   // Called once at the end of the query
 548 }