ANALYSIS/AliAnalysisTaskPIDCombined.cxx

   1 /*************************************************************************\r
   2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
   3 *                                                                        *\r
   4 * Author: The ALICE Off-line Project.                                    *\r
   5 * Contributors are mentioned in the code where appropriate.              *\r
   6 *                                                                        *\r
   7 * Permission to use, copy, modify and distribute this software and its   *\r
   8 * documentation strictly for non-commercial purposes is hereby granted   *\r
   9 * without fee, provided that the above copyright notice appears in all   *\r
  10 * copies and that both the copyright notice and this permission notice   *\r
  11 * appear in the supporting documentation. The authors make no claims     *\r
  12 * about the suitability of this software for any purpose. It is          *\r
  13 * provided "as is" without express or implied warranty.                  *\r
  14 **************************************************************************/\r
  15 \r
  16 ///////////////////////////////////////////////////////////////////////////\r
  17 //                                                                       //\r
  18 //                        Basic Analysis Task                            //\r
  19 //                      for PID        Analysis                          //\r
  20 //                                                                       //\r
  21 ///////////////////////////////////////////////////////////////////////////\r
  22 \r
  23 #include <TChain.h>\r
  24 #include <TH1D.h>\r
  25 #include <TH2D.h>\r
  26 \r
  27 #include <AliVEvent.h>\r
  28 #include <AliInputEventHandler.h>\r
  29 #include <AliAnalysisManager.h>\r
  30 \r
  31 #include <AliLog.h>\r
  32 #include <AliPID.h>\r
  33 #include <AliPIDCombined.h>\r
  34 #include <AliPIDResponse.h>\r
  35 \r
  36 #include "AliAnalysisTaskPIDCombined.h"\r
  37 \r
  38 const char *AliAnalysisTaskPIDCombined::fgkBinMomDesc[AliAnalysisTaskPIDCombined::kPtBins] = {\r
  39   " 0 <= p < 0.5 GeV/c",\r
  40   " 0.5 <= p < 0.7 GeV/c",\r
  41   " 0.7 <= p < 1.0 GeV/c",\r
  42   " 1.0 <= p < 1.5 GeV/c",\r
  43   " 1.5 <= p < 2.0 GeV/c",\r
  44   " p >= 2.0 GeV/c"\r
  45 };\r
  46 \r
  47 ClassImp(AliAnalysisTaskPIDCombined)\r
  48 \r
  49 //_________________________________________________________________________________\r
  50 AliAnalysisTaskPIDCombined::AliAnalysisTaskPIDCombined() :\r
  51   AliAnalysisTaskSE(),\r
  52   fHistList(),\r
  53   fProbTPCnSigma(),\r
  54   fProbTOFnSigma(),\r
  55   fProbTPCTOFnSigmaTPC(),\r
  56   fProbTPC(),\r
  57   fProbTOF(),\r
  58   fProbTPCTOF(),\r
  59   fPriors(),\r
  60   fProbTPCTOFnSigTPCMom(),\r
  61   fProbTPCnSigTPCMom(),\r
  62   fProbTOFnSigTOFMom(),\r
  63   fPriorsUsed(),\r
  64   fPIDResponse(0x0),\r
  65   fPIDCombined(0x0),\r
  66   fTrackCuts(0x0),\r
  67   fTrackFilter(0x0),\r
  68   fDeDx(NULL),
  69   fDeDxTuned(NULL)
  70 {\r
  71   //\r
  72   // Constructor\r
  73   //\r
  74 }\r
  75 \r
  76 //_________________________________________________________________________________\r
  77 AliAnalysisTaskPIDCombined::AliAnalysisTaskPIDCombined(const char *name) :\r
  78   AliAnalysisTaskSE(name),\r
  79   fHistList(),\r
  80   fProbTPCnSigma(),\r
  81   fProbTOFnSigma(),\r
  82   fProbTPCTOFnSigmaTPC(),\r
  83   fProbTPC(),\r
  84   fProbTOF(),\r
  85   fProbTPCTOF(),\r
  86   fPriors(),\r
  87   fProbTPCTOFnSigTPCMom(),\r
  88   fProbTPCnSigTPCMom(),\r
  89   fProbTOFnSigTOFMom(),\r
  90   fPriorsUsed(),\r
  91   fPIDResponse(0x0),\r
  92   fPIDCombined(0x0),\r
  93   fTrackCuts(0x0),\r
  94   fTrackFilter(0x0),\r
  95   fDeDx(NULL),
  96   fDeDxTuned(NULL)
  97 {\r
  98   //\r
  99   // Constructor\r
 100   //\r
 101   DefineInput(0,TChain::Class());\r
 102   DefineOutput(1, TList::Class());\r
 103 }\r
 104 \r
 105 //_________________________________________________________________________________\r
 106 void AliAnalysisTaskPIDCombined::UserCreateOutputObjects()\r
 107 {\r
 108   //\r
 109   // Initialise the framework objects\r
 110   //\r
 111 \r
 112 \r
 113   // ------- track cuts\r
 114   fTrackCuts = new AliESDtrackCuts("fTrackCuts", "Standard");\r
 115   fTrackCuts->SetAcceptKinkDaughters(kFALSE);\r
 116   fTrackCuts->SetMinNClustersTPC(80);\r
 117   fTrackCuts->SetMaxChi2PerClusterTPC(4);\r
 118   fTrackCuts->SetMaxDCAToVertexXY(3);\r
 119   fTrackCuts->SetMaxDCAToVertexZ(3);\r
 120   fTrackCuts->SetRequireTPCRefit(kTRUE);\r
 121   fTrackCuts->SetRequireITSRefit(kTRUE);\r
 122   fTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);\r
 123   fTrackFilter = new AliAnalysisFilter("trackFilter");\r
 124   fTrackFilter->AddCuts(fTrackCuts);\r
 125 \r
 126 \r
 127 \r
 128   // ------- setup PIDCombined\r
 129   fPIDCombined=new AliPIDCombined;\r
 130   fPIDCombined->SetDefaultTPCPriors();\r
 131   fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);\r
 132 \r
 133   // no light nuclei - no need to call it, this is default\r
 134   //  fPIDCombined->SetSelectedSpecies(AliPID::kSPECIES);\r
 135 \r
 136 \r
 137   fHistList.Add(new TH1D("nEvents","Number of Evnets;Selection",2,0,2));\r
 138 \r
 139   for (Int_t ispec=0; ispec<5; ++ispec){\r
 140 \r
 141  \r
 142     fProbTPC[ispec]=new TH2D(Form("prob%s_mom_TPC",AliPID::ParticleName(ispec)),\r
 143                                    Form("%s probability vs. momentum;momentum;probability",AliPID::ParticleName(ispec)),\r
 144                                    100,0.,20.,50,0.,1.);\r
 145     fHistList.Add(fProbTPC[ispec]);\r
 146     fProbTPCnSigma[ispec]=new TH2D(Form("prob%s_nSigma_TPC",AliPID::ParticleName(ispec)),\r
 147                                    Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),\r
 148                                    20,-5.,5.,50,0.,1.);\r
 149     fHistList.Add(fProbTPCnSigma[ispec]);\r
 150 \r
 151     for (Int_t ibin=0;ibin<kPtBins;ibin++) {\r
 152       fProbTPCnSigTPCMom[ibin][ispec]=new TH2D(Form("prob%s_nSigma_TPC (%s)",AliPID::ParticleName(ispec),fgkBinMomDesc[ibin]),\r
 153                                                Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),\r
 154                                                20,-5.,5.,50,0.,1.);\r
 155       fHistList.Add(fProbTPCnSigTPCMom[ibin][ispec]);\r
 156     }\r
 157 \r
 158 \r
 159 \r
 160     fProbTOF[ispec]=new TH2D(Form("prob%s_mom_TOF",AliPID::ParticleName(ispec)),\r
 161                                    Form("%s probability vs. momentum;momentum;probability",AliPID::ParticleName(ispec)),\r
 162                                    100,0.,20.,50,0.,1.);\r
 163     fHistList.Add(fProbTOF[ispec]);\r
 164     fProbTOFnSigma[ispec]=new TH2D(Form("prob%s_nSigma_TOF",AliPID::ParticleName(ispec)),\r
 165                                    Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),\r
 166                                    20,-5.,5.,50,0.,1.);\r
 167     fHistList.Add(fProbTOFnSigma[ispec]);\r
 168     for (Int_t ibin=0;ibin<kPtBins;ibin++) {\r
 169       fProbTOFnSigTOFMom[ibin][ispec]=new TH2D(Form("prob%s_nSigma_TOF (%s)",AliPID::ParticleName(ispec),fgkBinMomDesc[ibin]),\r
 170                                                Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),\r
 171                                                20,-5.,5.,50,0.,1.);\r
 172       fHistList.Add(fProbTOFnSigTOFMom[ibin][ispec]);\r
 173     }\r
 174 \r
 175 \r
 176 \r
 177     fProbTPCTOF[ispec]=new TH2D(Form("prob%s_mom_TPCTOF",AliPID::ParticleName(ispec)),\r
 178                                    Form("%s probability vs. momentum;momentum;probability",AliPID::ParticleName(ispec)),\r
 179                                    100,0.,20.,50,0.,1.);\r
 180     fHistList.Add(fProbTPCTOF[ispec]);\r
 181     fProbTPCTOFnSigmaTPC[ispec]=new TH2D(Form("probTPCTOF%s_nSigma_TPC",AliPID::ParticleName(ispec)),\r
 182                                    Form("%s TPCTOF probability vs. n#sigmaTPC;n#sigma;probability",AliPID::ParticleName(ispec)),\r
 183                                    20,-5.,5.,50,0.,1.);\r
 184     fHistList.Add(fProbTPCTOFnSigmaTPC[ispec]);\r
 185     for (Int_t ibin=0;ibin<kPtBins;ibin++) {\r
 186       fProbTPCTOFnSigTPCMom[ibin][ispec]=new TH2D(Form("probTPCTOF%s_nSigma_TPC (%s)",AliPID::ParticleName(ispec),fgkBinMomDesc[ibin]),\r
 187                                                Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),\r
 188                                                20,-5.,5.,50,0.,1.);\r
 189       fHistList.Add(fProbTPCTOFnSigTPCMom[ibin][ispec]);\r
 190     }\r
 191 \r
 192 \r
 193 \r
 194     // basic priors\r
 195     fPriors[ispec]=new TH1F(Form("%s_priors",AliPID::ParticleName(ispec)),\r
 196                             Form("%s priors vs momentum",AliPID::ParticleName(ispec)),\r
 197                             100,0.,20.);\r
 198     fHistList.Add(fPriors[ispec]);\r
 199     switch (ispec) {\r
 200     case AliPID::kElectron:\r
 201       for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.02);\r
 202       break;\r
 203     case AliPID::kMuon:\r
 204       for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.02);\r
 205       break;\r
 206     case AliPID::kPion:\r
 207       for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.56);\r
 208       break;\r
 209     case AliPID::kKaon:\r
 210       for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.20);\r
 211       break;\r
 212     case AliPID::kProton:\r
 213       for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.20);\r
 214       break;\r
 215     default:\r
 216       break;\r
 217     }\r
 218     fPIDCombined->SetPriorDistribution((AliPID::EParticleType)ispec,fPriors[ispec]);\r
 219 \r
 220     // priors used\r
 221     fPriorsUsed[ispec] = new TH2D(Form("%s_priorsUsed",AliPID::ParticleName(ispec)),\r
 222                             Form("%s priors vs transverse momentum;p_{t} (GeV/c);priors",AliPID::ParticleName(ispec)),\r
 223                                   100,0.,20.,101,0,1.01);      \r
 224     fHistList.Add(fPriorsUsed[ispec]);\r
 225   }\r
 226 \r
 227 \r
 228   fDeDx = new TH2D("hDeDx",";p_{TPC};dE/dx (a.u.)",500,0,5,500,0,500);
 229   fHistList.Add(fDeDx);
 230   fDeDxTuned = new TH2D("hDeDxTuned",";p_{TPC};dE/dx (a.u.)",500,0,5,500,0,500);
 231   fHistList.Add(fDeDxTuned);
 232
 233   fHistList.SetOwner();\r
 234   PostData(1,&fHistList);\r
 235 \r
 236 \r
 237 }\r
 238 \r
 239 //_________________________________________________________________________________\r
 240 void AliAnalysisTaskPIDCombined::UserExec(Option_t *)\r
 241 {\r
 242   //\r
 243   // Main loop. Called for every event\r
 244   //\r
 245   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();\r
 246   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());\r
 247   fPIDResponse=inputHandler->GetPIDResponse();\r
 248   if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");\r
 249 \r
 250   //  Printf(" ---------------------- UserExec PID task ---------------------");\r
 251   \r
 252   FillHistogram("nEvents",0.);\r
 253 \r
 254   AliVEvent *event=InputEvent();\r
 255   AliVTrack *track=0x0;\r
 256   Int_t ntracks=event->GetNumberOfTracks();\r
 257 \r
 258   Double_t probTPC[AliPID::kSPECIES]={0.};\r
 259   Double_t probTOF[AliPID::kSPECIES]={0.};\r
 260   Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r
 261   \r
 262   //loop over all tracks\r
 263   for (Int_t itrack=0; itrack<ntracks; ++itrack){\r
 264 \r
 265     track=(AliVTrack*)event->GetTrack(itrack);\r
 266 \r
 267     if ( fTrackFilter->IsSelected(track) ) {\r
 268 \r
 269       Double_t mom=track->GetTPCmomentum();\r
 270       Int_t ibin=GetMomBin(mom);\r
 271 \r
 272       fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r
 273       UInt_t detUsed = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);\r
 274      \r
 275       if (detUsed != 0) {  // TPC is available\r
 276         for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec) {\r
 277           Double_t nSigmaTPC = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)ispec);\r
 278           fProbTPC[ispec]->Fill(mom,probTPC[ispec]);\r
 279           fProbTPCnSigma[ispec]->Fill(nSigmaTPC,probTPC[ispec]);\r
 280           fProbTPCnSigTPCMom[ibin][ispec]->Fill(nSigmaTPC,probTPC[ispec]);\r
 281         }\r
 282 \r
 283         // compute priors for TPC+TOF, even if we ask just TOF for PID\r
 284         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r
 285         detUsed = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);\r
 286         Double_t priors[5];     // check priors used for TOF\r
 287         fPIDCombined->GetPriors(track,priors,fPIDResponse,detUsed);\r
 288         for(Int_t ispec=0;ispec<5;ispec++) fPriorsUsed[ispec]->Fill(TMath::Abs(track->Pt()),priors[ispec]);\r
 289 \r
 290         if (detUsed != 0) {  // TOF is available\r
 291           for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec) {\r
 292             Double_t nSigmaTOF = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)ispec);\r
 293             fProbTOF[ispec]->Fill(mom,probTOF[ispec]);\r
 294             fProbTOFnSigma[ispec]->Fill(nSigmaTOF,probTOF[ispec]);\r
 295             fProbTOFnSigTOFMom[ibin][ispec]->Fill(nSigmaTOF,probTOF[ispec]);\r
 296           }\r
 297         }\r
 298 \r
 299         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r
 300         detUsed = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);\r
 301         if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask() ) {\r
 302           for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec) {\r
 303             Double_t nSigmaTPC = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)ispec);\r
 304             fProbTPCTOF[ispec]->Fill(mom,probTPCTOF[ispec]);\r
 305             fProbTPCTOFnSigmaTPC[ispec]->Fill(nSigmaTPC,probTPCTOF[ispec]);\r
 306             fProbTPCTOFnSigTPCMom[ibin][ispec]->Fill(nSigmaTPC,probTPCTOF[ispec]);\r
 307           }\r
 308         }\r
 309 \r
 310       }\r
 311 \r
 312       fPIDResponse->GetTPCsignalTunedOnData(track);
 313
 314       fDeDx->Fill(mom,track->GetTPCsignal());
 315       fDeDxTuned->Fill(mom,track->GetTPCsignalTunedOnData());
 316 \r
 317     }\r
 318   }\r
 319 \r
 320   PostData(1, &fHistList);\r
 321 }\r
 322 \r
 323 //_________________________________________________________________________________\r
 324 void AliAnalysisTaskPIDCombined::FillHistogram(const char* name, Double_t x, Double_t weight)\r
 325 {\r
 326   //\r
 327   // Fill 1D histogram by name\r
 328   //\r
 329   ((TH1*)fHistList.FindObject(name))->Fill(x,weight);\r
 330 }\r
 331 \r
 332 //_________________________________________________________________________________\r
 333 void AliAnalysisTaskPIDCombined::FillHistogram(const char* name, Double_t x, Double_t y, Double_t weight)\r
 334 {\r
 335   //\r
 336   // Fill 2D histogram by name\r
 337   //\r
 338   ((TH2*)fHistList.FindObject(name))->Fill(x,y,weight);\r
 339 }\r
 340 \r
 341 \r
 342 //_________________________________________________________________________________\r
 343 Int_t AliAnalysisTaskPIDCombined::GetMomBin(Float_t mom)\r
 344 {\r
 345   //\r
 346   // Given momentum return histogram to be filled\r
 347   //\r
 348   if (mom>0. && mom < 0.5) return 0;\r
 349   if (mom>=0.5 && mom < 0.7) return 1;\r
 350   if (mom>=0.7 && mom < 1.0) return 2;\r
 351   if (mom>=1.0 && mom < 1.5) return 3;\r
 352   if (mom>=1.5 && mom < 2.0) return 4;\r
 353   return kPtBins-1;\r
 354 }\r
 355 \r