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   fPIDResponse(0x0),\r
  54   fPIDCombined(0x0),\r
  55   fTrackCuts(0x0),\r
  56   fTrackFilter(0x0)\r
  57 {\r
  58   //\r
  59   // Constructor\r
  60   //\r
  61 }\r
  62 \r
  63 //_________________________________________________________________________________\r
  64 AliAnalysisTaskPIDCombined::AliAnalysisTaskPIDCombined(const char *name) :\r
  65   AliAnalysisTaskSE(name),\r
  66   fHistList(),\r
  67   fPIDResponse(0x0),\r
  68   fPIDCombined(0x0),\r
  69   fTrackCuts(0x0),\r
  70   fTrackFilter(0x0)\r
  71 {\r
  72   //\r
  73   // Constructor\r
  74   //\r
  75   DefineInput(0,TChain::Class());\r
  76   DefineOutput(1, TList::Class());\r
  77 }\r
  78 \r
  79 //_________________________________________________________________________________\r
  80 void AliAnalysisTaskPIDCombined::UserCreateOutputObjects()\r
  81 {\r
  82   //\r
  83   // Initialise the framework objects\r
  84   //\r
  85 \r
  86 \r
  87   // ------- track cuts\r
  88   fTrackCuts = new AliESDtrackCuts("fTrackCuts", "Standard");\r
  89   fTrackCuts->SetAcceptKinkDaughters(kFALSE);\r
  90   fTrackCuts->SetMinNClustersTPC(80);\r
  91   fTrackCuts->SetMaxChi2PerClusterTPC(4);\r
  92   fTrackCuts->SetMaxDCAToVertexXY(3);\r
  93   fTrackCuts->SetMaxDCAToVertexZ(3);\r
  94   fTrackCuts->SetRequireTPCRefit(kTRUE);\r
  95   fTrackCuts->SetRequireITSRefit(kTRUE);\r
  96   fTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);\r
  97   fTrackFilter = new AliAnalysisFilter("trackFilter");\r
  98   fTrackFilter->AddCuts(fTrackCuts);\r
  99 \r
 100 \r
 101 \r
 102   // ------- setup PIDCombined\r
 103   fPIDCombined=new AliPIDCombined;\r
 104   fPIDCombined->SetDefaultTPCPriors();\r
 105   fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);\r
 106 \r
 107   // no light nuclei\r
 108   fPIDCombined->SetSelectedSpecies(AliPID::kSPECIES);\r
 109 \r
 110 \r
 111   fHistList.Add(new TH1D("nEvents","Number of Evnets;Selection",2,0,2));\r
 112 \r
 113   for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){\r
 114 \r
 115  \r
 116     fProbTPC[ispec]=new TH2D(Form("prob%s_mom_TPC",AliPID::ParticleName(ispec)),\r
 117                                    Form("%s probability vs. momentum;momentum;probability",AliPID::ParticleName(ispec)),\r
 118                                    100,0.,20.,50,0.,1.);\r
 119     fHistList.Add(fProbTPC[ispec]);\r
 120     fProbTPCnSigma[ispec]=new TH2D(Form("prob%s_nSigma_TPC",AliPID::ParticleName(ispec)),\r
 121                                    Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),\r
 122                                    20,-5.,5.,50,0.,1.);\r
 123     fHistList.Add(fProbTPCnSigma[ispec]);\r
 124 \r
 125     for (Int_t ibin=0;ibin<kPtBins;ibin++) {\r
 126       fProbTPCnSigTPCMom[ibin][ispec]=new TH2D(Form("prob%s_nSigma_TPC (%s)",AliPID::ParticleName(ispec),fgkBinMomDesc[ibin]),\r
 127                                                Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),\r
 128                                                20,-5.,5.,50,0.,1.);\r
 129       fHistList.Add(fProbTPCnSigTPCMom[ibin][ispec]);\r
 130     }\r
 131 \r
 132 \r
 133 \r
 134     fProbTOF[ispec]=new TH2D(Form("prob%s_mom_TOF",AliPID::ParticleName(ispec)),\r
 135                                    Form("%s probability vs. momentum;momentum;probability",AliPID::ParticleName(ispec)),\r
 136                                    100,0.,20.,50,0.,1.);\r
 137     fHistList.Add(fProbTOF[ispec]);\r
 138     fProbTOFnSigma[ispec]=new TH2D(Form("prob%s_nSigma_TOF",AliPID::ParticleName(ispec)),\r
 139                                    Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),\r
 140                                    20,-5.,5.,50,0.,1.);\r
 141     fHistList.Add(fProbTOFnSigma[ispec]);\r
 142     for (Int_t ibin=0;ibin<kPtBins;ibin++) {\r
 143       fProbTOFnSigTOFMom[ibin][ispec]=new TH2D(Form("prob%s_nSigma_TOF (%s)",AliPID::ParticleName(ispec),fgkBinMomDesc[ibin]),\r
 144                                                Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),\r
 145                                                20,-5.,5.,50,0.,1.);\r
 146       fHistList.Add(fProbTOFnSigTOFMom[ibin][ispec]);\r
 147     }\r
 148 \r
 149 \r
 150 \r
 151     fProbTPCTOF[ispec]=new TH2D(Form("prob%s_mom_TPCTOF",AliPID::ParticleName(ispec)),\r
 152                                    Form("%s probability vs. momentum;momentum;probability",AliPID::ParticleName(ispec)),\r
 153                                    100,0.,20.,50,0.,1.);\r
 154     fHistList.Add(fProbTPCTOF[ispec]);\r
 155     fProbTPCTOFnSigmaTPC[ispec]=new TH2D(Form("probTPCTOF%s_nSigma_TPC",AliPID::ParticleName(ispec)),\r
 156                                    Form("%s TPCTOF probability vs. n#sigmaTPC;n#sigma;probability",AliPID::ParticleName(ispec)),\r
 157                                    20,-5.,5.,50,0.,1.);\r
 158     fHistList.Add(fProbTPCTOFnSigmaTPC[ispec]);\r
 159     for (Int_t ibin=0;ibin<kPtBins;ibin++) {\r
 160       fProbTPCTOFnSigTPCMom[ibin][ispec]=new TH2D(Form("probTPCTOF%s_nSigma_TPC (%s)",AliPID::ParticleName(ispec),fgkBinMomDesc[ibin]),\r
 161                                                Form("%s probability vs. n#sigma;n#sigma;probability",AliPID::ParticleName(ispec)),\r
 162                                                20,-5.,5.,50,0.,1.);\r
 163       fHistList.Add(fProbTPCTOFnSigTPCMom[ibin][ispec]);\r
 164     }\r
 165 \r
 166 \r
 167 \r
 168     // basic priors\r
 169     fPriors[ispec]=new TH1F(Form("%s_priors",AliPID::ParticleName(ispec)),\r
 170                             Form("%s priors vs momentum",AliPID::ParticleName(ispec)),\r
 171                             100,0.,20.);\r
 172     fHistList.Add(fPriors[ispec]);\r
 173     switch (ispec) {\r
 174     case AliPID::kElectron:\r
 175       for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.02);\r
 176       break;\r
 177     case AliPID::kMuon:\r
 178       for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.02);\r
 179       break;\r
 180     case AliPID::kPion:\r
 181       for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.56);\r
 182       break;\r
 183     case AliPID::kKaon:\r
 184       for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.20);\r
 185       break;\r
 186     case AliPID::kProton:\r
 187       for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.20);\r
 188       break;\r
 189     default:\r
 190       break;\r
 191     }\r
 192     fPIDCombined->SetPriorDistribution((AliPID::EParticleType)ispec,fPriors[ispec]);\r
 193 \r
 194     // priors used\r
 195     fPriorsUsed[ispec] = new TH2D(Form("%s_priorsUsed",AliPID::ParticleName(ispec)),\r
 196                             Form("%s priors vs transverse momentum;p_{t} (GeV/c);priors",AliPID::ParticleName(ispec)),\r
 197                                   100,0.,20.,101,0,1.01);      \r
 198     fHistList.Add(fPriorsUsed[ispec]);\r
 199   }\r
 200 \r
 201 \r
 202 \r
 203   fHistList.SetOwner();\r
 204   PostData(1,&fHistList);\r
 205 \r
 206 \r
 207 }\r
 208 \r
 209 //_________________________________________________________________________________\r
 210 void AliAnalysisTaskPIDCombined::UserExec(Option_t *)\r
 211 {\r
 212   //\r
 213   // Main loop. Called for every event\r
 214   //\r
 215   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();\r
 216   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());\r
 217   fPIDResponse=inputHandler->GetPIDResponse();\r
 218   if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");\r
 219 \r
 220   //  Printf(" ---------------------- UserExec PID task ---------------------");\r
 221   \r
 222   FillHistogram("nEvents",0.);\r
 223 \r
 224   AliVEvent *event=InputEvent();\r
 225   AliVTrack *track=0x0;\r
 226   Int_t ntracks=event->GetNumberOfTracks();\r
 227 \r
 228   Double_t probTPC[AliPID::kSPECIES]={0.};\r
 229   Double_t probTOF[AliPID::kSPECIES]={0.};\r
 230   Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r
 231   \r
 232   //loop over all tracks\r
 233   for (Int_t itrack=0; itrack<ntracks; ++itrack){\r
 234 \r
 235     track=(AliVTrack*)event->GetTrack(itrack);\r
 236 \r
 237     if ( fTrackFilter->IsSelected(track) ) {\r
 238 \r
 239       Double_t mom=track->GetTPCmomentum();\r
 240       Int_t ibin=GetMomBin(mom);\r
 241 \r
 242       fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r
 243       UInt_t detUsed = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);\r
 244      \r
 245       if (detUsed != 0) {  // TPC is available\r
 246         for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec) {\r
 247           Double_t nSigmaTPC = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)ispec);\r
 248           fProbTPC[ispec]->Fill(mom,probTPC[ispec]);\r
 249           fProbTPCnSigma[ispec]->Fill(nSigmaTPC,probTPC[ispec]);\r
 250           fProbTPCnSigTPCMom[ibin][ispec]->Fill(nSigmaTPC,probTPC[ispec]);\r
 251         }\r
 252 \r
 253         // compute priors for TPC+TOF, even if we ask just TOF for PID\r
 254         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r
 255         detUsed = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);\r
 256         Double_t priors[5];     // check priors used for TOF\r
 257         fPIDCombined->GetPriors(track,priors,fPIDResponse,detUsed);\r
 258         for(Int_t ispec=0;ispec<5;ispec++) fPriorsUsed[ispec]->Fill(TMath::Abs(track->Pt()),priors[ispec]);\r
 259 \r
 260         if (detUsed != 0) {  // TOF is available\r
 261           for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec) {\r
 262             Double_t nSigmaTOF = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)ispec);\r
 263             fProbTOF[ispec]->Fill(mom,probTOF[ispec]);\r
 264             fProbTOFnSigma[ispec]->Fill(nSigmaTOF,probTOF[ispec]);\r
 265             fProbTOFnSigTOFMom[ibin][ispec]->Fill(nSigmaTOF,probTOF[ispec]);\r
 266           }\r
 267         }\r
 268 \r
 269         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r
 270         detUsed = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);\r
 271         if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask() ) {\r
 272           for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec) {\r
 273             Double_t nSigmaTPC = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)ispec);\r
 274             fProbTPCTOF[ispec]->Fill(mom,probTPCTOF[ispec]);\r
 275             fProbTPCTOFnSigmaTPC[ispec]->Fill(nSigmaTPC,probTPCTOF[ispec]);\r
 276             fProbTPCTOFnSigTPCMom[ibin][ispec]->Fill(nSigmaTPC,probTPCTOF[ispec]);\r
 277           }\r
 278         }\r
 279 \r
 280       }\r
 281 \r
 282 \r
 283     }\r
 284   }\r
 285 \r
 286   PostData(1, &fHistList);\r
 287 }\r
 288 \r
 289 //_________________________________________________________________________________\r
 290 void AliAnalysisTaskPIDCombined::FillHistogram(const char* name, Double_t x, Double_t weight)\r
 291 {\r
 292   //\r
 293   // Fill 1D histogram by name\r
 294   //\r
 295   ((TH1*)fHistList.FindObject(name))->Fill(x,weight);\r
 296 }\r
 297 \r
 298 //_________________________________________________________________________________\r
 299 void AliAnalysisTaskPIDCombined::FillHistogram(const char* name, Double_t x, Double_t y, Double_t weight)\r
 300 {\r
 301   //\r
 302   // Fill 2D histogram by name\r
 303   //\r
 304   ((TH2*)fHistList.FindObject(name))->Fill(x,y,weight);\r
 305 }\r
 306 \r
 307 \r
 308 //_________________________________________________________________________________\r
 309 Int_t AliAnalysisTaskPIDCombined::GetMomBin(Float_t mom)\r
 310 {\r
 311   //\r
 312   // Given momentum return histogram to be filled\r
 313   //\r
 314   if (mom>0. && mom < 0.5) return 0;\r
 315   if (mom>=0.5 && mom < 0.7) return 1;\r
 316   if (mom>=0.7 && mom < 1.0) return 2;\r
 317   if (mom>=1.0 && mom < 1.5) return 3;\r
 318   if (mom>=1.5 && mom < 2.0) return 4;\r
 319   return kPtBins-1;\r
 320 }\r
 321 \r