]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/Tools/AliHelperPID.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
class namechange
[u/mrichter/AliRoot.git] / PWG / Tools / AliHelperPID.cxx
CommitLineData
32e49607 1
2/**************************************************************************
3 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4 * *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
7 * *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
16
17//-----------------------------------------------------------------
18// AliHelperPID class
19//-----------------------------------------------------------------
20
21#include "AliHelperPID.h"
22#include "AliVEvent.h"
23#include "AliAODEvent.h"
24#include "AliMCEvent.h"
25#include "AliMCEventHandler.h"
26#include "TH1F.h"
27#include "TH2F.h"
28#include "TList.h"
29#include "AliStack.h"
30#include "AliVTrack.h"
31#include "TParticle.h"
32#include "AliAODMCParticle.h"
33#include "AliPIDResponse.h"
d2062f6b 34#include "AliPIDCombined.h"
32e49607 35#include "AliAnalysisManager.h"
36#include "AliInputEventHandler.h"
37
38using namespace AliHelperPIDNameSpace;
39using namespace std;
40
41const char * kPIDTypeName[]={"TPC","TOF","TPC-TOF"} ;
d2062f6b 42const char * kDetectorName[]={"ITS","TPC","TOF"} ;
32e49607 43const char * kParticleSpeciesName[]={"Pions","Kaons","Protons","Undefined"} ;
44
45ClassImp(AliHelperPID)
46
1fdecdaf 47AliHelperPID::AliHelperPID() : TNamed("HelperPID", "PID object"),fisMC(0), fPIDType(kNSigmaTPCTOF), fNSigmaPID(3), fBayesCut(0.8), fPIDResponse(0x0), fPIDCombined(0x0),fOutputList(0x0),fRequestTOFPID(1),fRemoveTracksT0Fill(0),fUseExclusiveNSigma(0),fPtTOFPID(.6),fHasTOFPID(0){
32e49607 48
49 for(Int_t ipart=0;ipart<kNSpecies;ipart++)
50 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++)
51 fnsigmas[ipart][ipid]=999.;
52
53 for(Int_t ipart=0;ipart<kNSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;
54
55 fOutputList = new TList;
56 fOutputList->SetOwner();
57 fOutputList->SetName("OutputList");
58
59 //nsigma plot
60 for(Int_t ipart=0;ipart<kNSpecies;ipart++){
61 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
62 Double_t miny=-30;
63 Double_t maxy=30;
64 if(ipid==kNSigmaTPCTOF){miny=0;maxy=50;}
d2062f6b 65 TH2F *fHistoNSigma=new TH2F(Form("NSigma_%d_%d",ipart,ipid),Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
32e49607 66 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
67 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
68 fOutputList->Add(fHistoNSigma);
69 }
70 }
71
72 //nsigmaRec plot
73 for(Int_t ipart=0;ipart<kNSpecies;ipart++){
74 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
75 Double_t miny=-10;
76 Double_t maxy=10;
77 if(ipid==kNSigmaTPCTOF){miny=0;maxy=20;}
78 TH2F *fHistoNSigma=new TH2F(Form("NSigmaRec_%d_%d",ipart,ipid),
d2062f6b 79 Form("n#sigma for reconstructed %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
32e49607 80 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
81 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
82 fOutputList->Add(fHistoNSigma);
83 }
84 }
85
d2062f6b 86 //BayesRec plot
87 for(Int_t ipart=0;ipart<kNSpecies;ipart++){
88 Double_t miny=0.;
89 Double_t maxy=1;
90 TH2F *fHistoBayes=new TH2F(Form("BayesRec_%d",ipart),
91 Form("probability for reconstructed %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
92 fHistoBayes->GetXaxis()->SetTitle("P_{T} (GeV / c)");
93 fHistoBayes->GetYaxis()->SetTitle(Form("Bayes prob %s",kParticleSpeciesName[ipart]));
94 fOutputList->Add(fHistoBayes);
95 }
96
32e49607 97 //nsigmaDC plot
98 for(Int_t ipart=0;ipart<kNSpecies;ipart++){
99 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
100 Double_t miny=-10;
101 Double_t maxy=10;
102 if(ipid==kNSigmaTPCTOF){miny=0;maxy=20;}
103 TH2F *fHistoNSigma=new TH2F(Form("NSigmaDC_%d_%d",ipart,ipid),
d2062f6b 104 Form("n#sigma for double counting %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
32e49607 105 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
106 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
107 fOutputList->Add(fHistoNSigma);
108 }
109 }
110
111 //nsigmaMC plot
112 for(Int_t ipart=0;ipart<kNSpecies;ipart++){
113 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
114 Double_t miny=-30;
115 Double_t maxy=30;
116 if(ipid==kNSigmaTPCTOF){miny=0;maxy=50;}
117 TH2F *fHistoNSigma=new TH2F(Form("NSigmaMC_%d_%d",ipart,ipid),
d2062f6b 118 Form("n#sigma for MC %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
32e49607 119 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
120 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
121 fOutputList->Add(fHistoNSigma);
122 }
123 }
124
125 //PID signal plot
126 for(Int_t idet=0;idet<kNDetectors;idet++){
d2062f6b 127 for(Int_t ipart=0;ipart<kNSpecies;ipart++){
128 Double_t maxy=500;
129 if(idet==kTOF)maxy=1.1;
130 TH2F *fHistoPID=new TH2F(Form("PID_%d_%d",idet,ipart),Form("%s signal - %s",kDetectorName[idet],kParticleSpeciesName[ipart]),200,0,10,500,-maxy,maxy);
131 fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
132 fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
133 fOutputList->Add(fHistoPID);
134 }
32e49607 135 }
214084ec 136 //PID signal plot, before PID cut
137 for(Int_t idet=0;idet<kNDetectors;idet++){
138 Double_t maxy=500;
139 if(idet==kTOF)maxy=1.1;
140 TH2F *fHistoPID=new TH2F(Form("PIDAll_%d",idet),Form("%s signal",kDetectorName[idet]),200,0,10,500,-maxy,maxy);
141 fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
142 fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
143 fOutputList->Add(fHistoPID);
144 }
32e49607 145}
146
147//////////////////////////////////////////////////////////////////////////////////////////////////
148
149TH2F* AliHelperPID::GetHistogram2D(const char * name){
150 // returns histo named name
151 return (TH2F*) fOutputList->FindObject(name);
152}
153
154//////////////////////////////////////////////////////////////////////////////////////////////////
155
156Int_t AliHelperPID::GetParticleSpecies(AliVTrack * trk, Bool_t FIllQAHistos){
157 //return the specie according to the minimum nsigma value
158 //no double counting, this has to be evaluated using CheckDoubleCounting()
32e49607 159
d2062f6b 160 Int_t ID=kSpUndefined;
161
162 //get the PID response
163 if(!fPIDResponse) {
164 AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
165 AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
166 fPIDResponse = inputHandler->GetPIDResponse();
167 }
168 if(!fPIDResponse) {
169 AliFatal("Cannot get pid response");
170 }
171
172 //calculate nsigmas (used also by the bayesian)
32e49607 173 CalculateNSigmas(trk,FIllQAHistos);//fill the data member fnsigmas with the nsigmas value [ipart][iPID]
174
d2062f6b 175 //Do PID
176 if(fPIDType==kBayes){//use bayesianPID
177
178 if(!fPIDCombined) {
52fffa3c 179 // ------- setup PIDCombined
180 fPIDCombined=new AliPIDCombined;
181 fPIDCombined->SetDefaultTPCPriors();
182 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
183 }
184 if(!fPIDCombined) {
1fdecdaf 185 AliFatal("PIDCombined object not found");
d2062f6b 186 }
187
188 ID = GetIDBayes(trk,FIllQAHistos);
189
190 }else{ //use nsigma PID
191
192 ID=FindMinNSigma(trk,FIllQAHistos);
193
194 if(fUseExclusiveNSigma){ //if one particle has double counting and exclusive nsigma is requested ID = kSpUndefined
195 Bool_t *HasDC;
196 HasDC=GetDoubleCounting(trk,FIllQAHistos);
197 for(Int_t ipart=0;ipart<kNSpecies;ipart++){
198 if(HasDC[ipart]==kTRUE) ID = kSpUndefined;
199 }
5838bc91 200 }
201 }
32e49607 202
d2062f6b 203 //Fill PID signal plot
204 if(ID != kSpUndefined){
205 for(Int_t idet=0;idet<kNDetectors;idet++){
206 TH2F *h=GetHistogram2D(Form("PID_%d_%d",idet,ID));
207 if(idet==kITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
208 if(idet==kTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
209 if(idet==kTOF && fHasTOFPID)h->Fill(trk->P(),TOFBetaCalc(trk)*trk->Charge());
210 }
211 }
214084ec 212 //Fill PID signal plot without cuts
213 for(Int_t idet=0;idet<kNDetectors;idet++){
214 TH2F *h=GetHistogram2D(Form("PIDAll_%d",idet));
215 if(idet==kITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
216 if(idet==kTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
217 if(idet==kTOF && fHasTOFPID)h->Fill(trk->P(),TOFBetaCalc(trk)*trk->Charge());
218 }
219
d2062f6b 220 return ID;
32e49607 221}
222
223//////////////////////////////////////////////////////////////////////////////////////////////////
224
7d3a3d9b 225Int_t AliHelperPID::GetParticleSpecies(AliVParticle * part) {
226 // return PID according to MC truth
227 switch(TMath::Abs(part->PdgCode())){
228 case 2212:
229 return kSpProton;
230 break;
231 case 321:
232 return kSpKaon;
233 break;
234 case 211:
235 return kSpPion;
236 break;
237 default:
238 return kSpUndefined;
239 }
240}
241
242//////////////////////////////////////////////////////////////////////////////////////////////////
243
d2062f6b 244Int_t AliHelperPID::GetIDBayes(AliVTrack * trk, Bool_t FIllQAHistos){
245
246 Bool_t *IDs=GetAllCompatibleIdentitiesNSigma(trk,FIllQAHistos);
247
248 Double_t probBayes[AliPID::kSPECIES];
249
250 UInt_t detUsed= 0;
251 CheckTOF(trk);
252 if(fHasTOFPID && trk->Pt()>fPtTOFPID){//use TOF information
253 detUsed = CalcPIDCombined(trk, fPIDResponse, AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC, probBayes);
254 if(detUsed != (AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC))return kSpUndefined;//check that TPC and TOF are used
255 }else{
256 detUsed = CalcPIDCombined(trk, fPIDResponse,AliPIDResponse::kDetTPC, probBayes);
257 if(detUsed != AliPIDResponse::kDetTPC)return kSpUndefined;//check that TPC is used
32e49607 258 }
d2062f6b 259
260 //the probability has to be normalized to one, we check it
261 Double_t sump=0.;
262 for(Int_t ipart=0;ipart<AliPID::kSPECIES;ipart++)sump+=probBayes[ipart];
263 if(sump<.99 && sump>1.01){//FIXME precision problem in the sum, workaround
264 AliFatal("Bayesian probability not normalized to one");
265 }
266
267 //probabilities are normalized to one, if the cut is above .5 there is no problem
268 if(probBayes[AliPID::kPion]>fBayesCut && IDs[kSpPion]==1){
269 TH2F *h=GetHistogram2D(Form("BayesRec_%d",kSpPion));
270 h->Fill(trk->Pt(),probBayes[AliPID::kPion]);
271 return kSpPion;
32e49607 272 }
d2062f6b 273 else if(probBayes[AliPID::kKaon]>fBayesCut && IDs[kSpKaon]==1){
274 TH2F *h=GetHistogram2D(Form("BayesRec_%d",kSpKaon));
275 h->Fill(trk->Pt(),probBayes[AliPID::kKaon]);
276 return kSpKaon;
277 }
278 else if(probBayes[AliPID::kProton]>fBayesCut && IDs[kSpProton]==1){
279 TH2F *h=GetHistogram2D(Form("BayesRec_%d",kSpProton));
280 h->Fill(trk->Pt(),probBayes[AliPID::kProton]);
281 return kSpProton;
282 }
283 else{
284 return kSpUndefined;
285 }
286}
287
288//////////////////////////////////////////////////////////////////////////////////////////////////
289
290UInt_t AliHelperPID::CalcPIDCombined(const AliVTrack *track,const AliPIDResponse *PIDResponse, Int_t detMask, Double_t* prob) const{
291 //
292 // Bayesian PID calculation
293 //
294 for(Int_t i=0;i<AliPID::kSPECIES;i++)
295 {
296 prob[i]=0.;
297 }
298 fPIDCombined->SetDetectorMask(detMask);
299
300 return fPIDCombined->ComputeProbabilities((AliVTrack*)track, PIDResponse, prob);
301}
302
303//////////////////////////////////////////////////////////////////////////////////////////////////
304
305void AliHelperPID::CalculateNSigmas(AliVTrack * trk, Bool_t FIllQAHistos){
306 //defines data member fnsigmas
32e49607 307
308 // Compute nsigma for each hypthesis
309 AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(trk);
310 // --- TPC
311 Double_t nsigmaTPCkProton = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton);
312 Double_t nsigmaTPCkKaon = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon);
313 Double_t nsigmaTPCkPion = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion);
314 // --- TOF
315 Double_t nsigmaTOFkProton=999.,nsigmaTOFkKaon=999.,nsigmaTOFkPion=999.;
316 Double_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;
317
318 CheckTOF(trk);
319
320 if(fHasTOFPID && trk->Pt()>fPtTOFPID){//use TOF information
321 nsigmaTOFkProton = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton);
322 nsigmaTOFkKaon = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon);
323 nsigmaTOFkPion = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion);
aed2199f 324 Double_t d2Proton=nsigmaTPCkProton * nsigmaTPCkProton + nsigmaTOFkProton * nsigmaTOFkProton;
325 Double_t d2Kaon=nsigmaTPCkKaon * nsigmaTPCkKaon + nsigmaTOFkKaon * nsigmaTOFkKaon;
326 Double_t d2Pion=nsigmaTPCkPion * nsigmaTPCkPion + nsigmaTOFkPion * nsigmaTOFkPion;
327 //commented, this is done in analogy with AliESDTrackCuts, nsigma combind according to the probability
32e49607 328 // --- combined
aed2199f 329 // -----------------------------------
330 // How to get to a n-sigma cut?
331 //
332 // The accumulated statistics from 0 to d is
333 //
334 // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
335 // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
336 // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
337 //
338 // work around precision problem
339 // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
340 //if(TMath::Exp(- d2Proton / 2) > 1e-15)nsigmaTPCTOFkProton = TMath::Sqrt(2)*TMath::ErfInverse(1 - TMath::Exp(- d2Proton / 2));
341 //if(TMath::Exp(- d2Kaon / 2) > 1e-15)nsigmaTPCTOFkKaon = TMath::Sqrt(2)*TMath::ErfInverse(1 - TMath::Exp(- d2Kaon / 2));
342 //if(TMath::Exp(- d2Pion / 2) > 1e-15)nsigmaTPCTOFkPion = TMath::Sqrt(2)*TMath::ErfInverse(1 - TMath::Exp(- d2Pion / 2));
343
344 //used for the 2PC PID paper (circular cut)
345 nsigmaTPCTOFkProton = TMath::Sqrt(d2Proton);
346 nsigmaTPCTOFkKaon = TMath::Sqrt(d2Kaon);
347 nsigmaTPCTOFkPion = TMath::Sqrt(d2Pion);
32e49607 348 }else{
349 // --- combined
350 // if TOF is missing and below fPtTOFPID only the TPC information is used
351 nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
352 nsigmaTPCTOFkKaon = TMath::Abs(nsigmaTPCkKaon);
353 nsigmaTPCTOFkPion = TMath::Abs(nsigmaTPCkPion);
354 }
355
356 //set data member fnsigmas
357 fnsigmas[kSpPion][kNSigmaTPC]=nsigmaTPCkPion;
358 fnsigmas[kSpKaon][kNSigmaTPC]=nsigmaTPCkKaon;
359 fnsigmas[kSpProton][kNSigmaTPC]=nsigmaTPCkProton;
360 fnsigmas[kSpPion][kNSigmaTOF]=nsigmaTOFkPion;
361 fnsigmas[kSpKaon][kNSigmaTOF]=nsigmaTOFkKaon;
362 fnsigmas[kSpProton][kNSigmaTOF]=nsigmaTOFkProton;
363 fnsigmas[kSpPion][kNSigmaTPCTOF]=nsigmaTPCTOFkPion;
364 fnsigmas[kSpKaon][kNSigmaTPCTOF]=nsigmaTPCTOFkKaon;
365 fnsigmas[kSpProton][kNSigmaTPCTOF]=nsigmaTPCTOFkProton;
366
367 if(FIllQAHistos){
368 //Fill NSigma SeparationPlot
369 for(Int_t ipart=0;ipart<kNSpecies;ipart++){
370 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
371 if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && trk->Pt()<fPtTOFPID)continue;//not filling TOF and combined if no TOF PID
372 TH2F *h=GetHistogram2D(Form("NSigma_%d_%d",ipart,ipid));
373 h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
374 }
375 }
32e49607 376 }
377}
378
379//////////////////////////////////////////////////////////////////////////////////////////////////
380
381Int_t AliHelperPID::FindMinNSigma(AliVTrack * trk,Bool_t FillQAHistos){
382
383 CheckTOF(trk);
384 if(fRequestTOFPID && (!fHasTOFPID) && trk->Pt()>fPtTOFPID)return kSpUndefined;
385
386 //get the identity of the particle with the minimum Nsigma
387 Double_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
388 switch (fPIDType){
389 case kNSigmaTPC:
390 nsigmaProton = TMath::Abs(fnsigmas[kSpProton][kNSigmaTPC]);
391 nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPC]) ;
392 nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTPC]) ;
393 break;
394 case kNSigmaTOF:
395 nsigmaProton = TMath::Abs(fnsigmas[kSpProton][kNSigmaTOF]);
396 nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTOF]) ;
397 nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTOF]) ;
398 break;
399 case kNSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
400 nsigmaProton = TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
401 nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF]) ;
402 nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF]) ;
403 break;
d2062f6b 404 case kBayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
405 nsigmaProton = TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
406 nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF]) ;
407 nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF]) ;
408 break;
32e49607 409 }
410
411 // guess the particle based on the smaller nsigma (within fNSigmaPID)
412 if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return kSpUndefined;//if is the default value for the three
413
414 if( ( nsigmaKaon < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon < fNSigmaPID)){
415 if(FillQAHistos){
416 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
417 if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
418 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",kSpKaon,ipid));
419 h->Fill(trk->Pt(),fnsigmas[kSpKaon][ipid]);
420 }
421 }
422 return kSpKaon;
423 }
424 if( ( nsigmaPion < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion < fNSigmaPID)){
425 if(FillQAHistos){
426 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
427 if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
428 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",kSpPion,ipid));
429 h->Fill(trk->Pt(),fnsigmas[kSpPion][ipid]);
430 }
431 }
432 return kSpPion;
433 }
434 if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)){
435 if(FillQAHistos){
436 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
437 if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
438 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",kSpProton,ipid));
439 h->Fill(trk->Pt(),fnsigmas[kSpProton][ipid]);
440 }
441 }
442 return kSpProton;
443 }
444 // else, return undefined
445 return kSpUndefined;
446}
447
448//////////////////////////////////////////////////////////////////////////////////////////////////
449
450Bool_t* AliHelperPID::GetDoubleCounting(AliVTrack * trk,Bool_t FIllQAHistos){
451 //if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE
452 //fill DC histos
8ccf185c 453 for(Int_t ipart=0;ipart<kNSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track
32e49607 454
455 Int_t MinNSigma=FindMinNSigma(trk,kFALSE);//not filling the NSigmaRec histos
456
457 CheckTOF(trk);
458
459 if(MinNSigma==kSpUndefined)return fHasDoubleCounting;//in case of undefined no Double counting
460
461 Double_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
462 switch (fPIDType) {
463 case kNSigmaTPC:
464 nsigmaProton = TMath::Abs(fnsigmas[kSpProton][kNSigmaTPC]);
465 nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPC]) ;
466 nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTPC]) ;
467 break;
468 case kNSigmaTOF:
469 nsigmaProton = TMath::Abs(fnsigmas[kSpProton][kNSigmaTOF]);
470 nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTOF]) ;
471 nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTOF]) ;
472 break;
473 case kNSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
474 nsigmaProton = TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
475 nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF]) ;
476 nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF]) ;
477 break;
d2062f6b 478 case kBayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
479 nsigmaProton = TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
480 nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF]) ;
481 nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF]) ;
482 break;
32e49607 483 }
484 if(nsigmaPion<fNSigmaPID && MinNSigma!=kSpPion)fHasDoubleCounting[kSpPion]=kTRUE;
485 if(nsigmaKaon<fNSigmaPID && MinNSigma!=kSpKaon)fHasDoubleCounting[kSpKaon]=kTRUE;
486 if(nsigmaProton<fNSigmaPID && MinNSigma!=kSpProton)fHasDoubleCounting[kSpProton]=kTRUE;
487
488 if(FIllQAHistos){
489 //fill NSigma distr for double counting
490 for(Int_t ipart=0;ipart<kNSpecies;ipart++){
491 if(fHasDoubleCounting[ipart]){
492 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
493 if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
494 TH2F *h=GetHistogram2D(Form("NSigmaDC_%d_%d",ipart,ipid));
495 h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
496 }
497 }
498 }
499 }
500 return fHasDoubleCounting;
501}
502
503//////////////////////////////////////////////////////////////////////////////////////////////////
504
d2062f6b 505Bool_t* AliHelperPID::GetAllCompatibleIdentitiesNSigma(AliVTrack * trk,Bool_t FIllQAHistos){
506
507 Bool_t *IDs=GetDoubleCounting(trk,FIllQAHistos);
508 IDs[FindMinNSigma(trk,FIllQAHistos)]=kTRUE;
509 return IDs;
510
511}
512
513//////////////////////////////////////////////////////////////////////////////////////////////////
514
32e49607 515Int_t AliHelperPID::GetMCParticleSpecie(AliVEvent* event, AliVTrack * trk, Bool_t FillQAHistos){
516 //return the specie according to the MC truth
517 CheckTOF(trk);
518
519 if(!fisMC)AliFatal("Error: AliHelperPID::GetMCParticleSpecie called on data\n");
520
521 Int_t code=999;
522 Bool_t isAOD=kFALSE;
523 Bool_t isESD=kFALSE;
524 TString nameoftrack(event->ClassName());
525 if(!nameoftrack.CompareTo("AliESDEvent"))isESD=kTRUE;
526 else if(!nameoftrack.CompareTo("AliAODEvent"))isAOD=kTRUE;
527 else AliFatal("Not processing AODs nor ESDs") ;
528
529 if(isAOD){
530 TClonesArray *arrayMC = 0;
531 arrayMC = (TClonesArray*) event->GetList()->FindObject(AliAODMCParticle::StdBranchName());
532 if (!arrayMC)AliFatal("Error: MC particles branch not found!\n");
533 AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(trk->GetLabel()));
534 if (!partMC){
535 AliError("Cannot get MC particle");
536 return kSpUndefined;
537 }
538 code=partMC->GetPdgCode();
539 }
540 if(isESD){
541 AliStack* stack=0x0;
542 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
543 if (!eventHandler) AliFatal("ERROR: Could not retrieve MC event handler");
544 AliMCEvent* mcEvent = eventHandler->MCEvent();
545 if (!mcEvent)AliFatal("ERROR: Could not retrieve MC event");
546 stack = mcEvent->Stack();
547 if (!stack) AliFatal("ERROR: stack not available\n");
548 TParticle *part=0;
549 part = (TParticle*)stack->Particle(TMath::Abs(trk->GetLabel()));
550 if (!part){
551 AliError("Cannot get MC particle");
552 return kSpUndefined;
553 }
554 code = part->GetPdgCode();
555 }
556 switch(TMath::Abs(code)){
557 case 2212:
558 if(FillQAHistos){
559 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
560 if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
561 TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",kSpProton,ipid));
562 h->Fill(trk->Pt(),fnsigmas[kSpProton][ipid]);
563 }
564 }
565 return kSpProton;
566 break;
567 case 321:
568 if(FillQAHistos){
569 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
570 if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
571 TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",kSpKaon,ipid));
572 h->Fill(trk->Pt(),fnsigmas[kSpKaon][ipid]);
573 }
574 }
575 return kSpKaon;
576 break;
577 case 211:
578 if(FillQAHistos){
579 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
580 if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
581 TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",kSpPion,ipid));
582 h->Fill(trk->Pt(),fnsigmas[kSpPion][ipid]);
583 }
584 }
585 return kSpPion;
586 break;
587 default:
588 return kSpUndefined;
589 }
590
591}
592
593//////////////////////////////////////////////////////////////////////////////////////////////////
594
595void AliHelperPID::CheckTOF(AliVTrack * trk)
596{
597 //check if the particle has TOF Matching
d2062f6b 598
599 //get the PIDResponse
600 if(fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF,trk)==0)fHasTOFPID=kFALSE;
32e49607 601 else fHasTOFPID=kTRUE;
d2062f6b 602
603 //in addition to TOF status we look at the pt
32e49607 604 if(trk->Pt()<fPtTOFPID)fHasTOFPID=kFALSE;
b46d749a 605
606 if(fRemoveTracksT0Fill)
607 {
b46d749a 608 Int_t startTimeMask = fPIDResponse->GetTOFResponse().GetStartTimeMask(trk->P());
609 if (startTimeMask < 0)fHasTOFPID=kFALSE;
610 }
32e49607 611}
612
613//////////////////////////////////////////////////////////////////////////////////////////////////
614
d2062f6b 615Double_t AliHelperPID::TOFBetaCalc(AliVTrack *track) const{
616 //TOF beta calculation
617 Double_t tofTime=track->GetTOFsignal();
618
619 Double_t c=TMath::C()*1.E-9;// m/ns
620 Float_t startTime = fPIDResponse->GetTOFResponse().GetStartTime(((AliVTrack*)track)->P());//in ps
621 Double_t length= fPIDResponse->GetTOFResponse().GetExpectedSignal(track,AliPID::kElectron)*1E-3*c;
622 tofTime -= startTime; // subtract startTime to the signal
623 Double_t tof= tofTime*1E-3; // ns, average T0 fill subtracted, no info from T0detector
624 tof=tof*c;
625 return length/tof;
626}
627
628//////////////////////////////////////////////////////////////////////////////////////////////////
629
7d3a3d9b 630Double_t AliHelperPID::GetMass(AliHelperParticleSpecies_t id) const{
631 //return Mass according to AliHelperParticleSpecies_t. If undefined return -999.
632 Double_t mass=-999.;
633 if (id == kSpProton) { mass=9.38271999999999995e-01; }
634 if (id == kSpKaon) { mass=4.93676999999999977e-01; }
635 if (id == kSpPion) { mass=1.39570000000000000e-01; }
636 return mass;
637}
638
639//////////////////////////////////////////////////////////////////////////////////////////////////
640
32e49607 641Long64_t AliHelperPID::Merge(TCollection* list)
642{
643 // Merging interface.
644 // Returns the number of merged objects (including this).
645 Printf("Merging");
646 if (!list)
647 return 0;
648 if (list->IsEmpty())
649 return 1;
650 TIterator* iter = list->MakeIterator();
651 TObject* obj;
652 // collections of TList with output histos
653 TList collections;
654 Int_t count = 0;
655 while ((obj = iter->Next())) {
656 AliHelperPID* entry = dynamic_cast<AliHelperPID*> (obj);
657 if (entry == 0)
658 continue;
659 TList * outputlist = entry->GetOutputList();
660 collections.Add(outputlist);
661 count++;
662 }
663 fOutputList->Merge(&collections);
664 delete iter;
665 Printf("OK");
666 return count+1;
667}