Bugfix: removing ","
[u/mrichter/AliRoot.git] / 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