ANALYSIS/AliAnalysisTaskPIDCombined.cxx

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