2 /**************************************************************************
3 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 //-----------------------------------------------------------------
19 //-----------------------------------------------------------------
21 #include "AliHelperPID.h"
22 #include "AliVEvent.h"
23 #include "AliAODEvent.h"
24 #include "AliMCEvent.h"
25 #include "AliMCEventHandler.h"
30 #include "AliVTrack.h"
31 #include "TParticle.h"
32 #include "AliAODMCParticle.h"
33 #include "AliPIDResponse.h"
34 #include "AliPIDCombined.h"
35 #include "AliAnalysisManager.h"
36 #include "AliInputEventHandler.h"
38 using namespace AliHelperPIDNameSpace;
41 const char * kPIDTypeName[]={"TPC","TOF","TPC-TOF"} ;
42 const char * kDetectorName[]={"ITS","TPC","TOF"} ;
43 const char * kParticleSpeciesName[]={"Pions","Kaons","Protons","Undefined"} ;
45 ClassImp(AliHelperPID)
47 AliHelperPID::AliHelperPID() : TNamed("HelperPID", "PID object"),fisMC(0), fPIDType(kNSigmaTPCTOF), fNSigmaPID(3), fBayesCut(0.8), fPIDResponse(0), fPIDCombined(0),fOutputList(0),fRequestTOFPID(1),fRemoveTracksT0Fill(0),fUseExclusiveNSigma(0),fPtTOFPID(.6),fHasTOFPID(0){
49 for(Int_t ipart=0;ipart<kNSpecies;ipart++)
50 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++)
51 fnsigmas[ipart][ipid]=999.;
53 for(Int_t ipart=0;ipart<kNSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;
55 fOutputList = new TList;
56 fOutputList->SetOwner();
57 fOutputList->SetName("OutputList");
60 for(Int_t ipart=0;ipart<kNSpecies;ipart++){
61 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
64 if(ipid==kNSigmaTPCTOF){miny=0;maxy=50;}
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);
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);
73 for(Int_t ipart=0;ipart<kNSpecies;ipart++){
74 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
77 if(ipid==kNSigmaTPCTOF){miny=0;maxy=20;}
78 TH2F *fHistoNSigma=new TH2F(Form("NSigmaRec_%d_%d",ipart,ipid),
79 Form("n#sigma for reconstructed %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
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);
87 for(Int_t ipart=0;ipart<kNSpecies;ipart++){
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);
98 for(Int_t ipart=0;ipart<kNSpecies;ipart++){
99 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
102 if(ipid==kNSigmaTPCTOF){miny=0;maxy=20;}
103 TH2F *fHistoNSigma=new TH2F(Form("NSigmaDC_%d_%d",ipart,ipid),
104 Form("n#sigma for double counting %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
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);
112 for(Int_t ipart=0;ipart<kNSpecies;ipart++){
113 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
116 if(ipid==kNSigmaTPCTOF){miny=0;maxy=50;}
117 TH2F *fHistoNSigma=new TH2F(Form("NSigmaMC_%d_%d",ipart,ipid),
118 Form("n#sigma for MC %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
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);
126 for(Int_t idet=0;idet<kNDetectors;idet++){
127 for(Int_t ipart=0;ipart<kNSpecies;ipart++){
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);
136 //PID signal plot, before PID cut
137 for(Int_t idet=0;idet<kNDetectors;idet++){
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);
148 //////////////////////////////////////////////////////////////////////////////////////////////////
150 TH2F* AliHelperPID::GetHistogram2D(const char * name){
151 // returns histo named name
152 return (TH2F*) fOutputList->FindObject(name);
155 //////////////////////////////////////////////////////////////////////////////////////////////////
157 Int_t AliHelperPID::GetParticleSpecies(AliVTrack * trk, Bool_t FIllQAHistos){
158 //return the specie according to the minimum nsigma value
159 //no double counting, this has to be evaluated using CheckDoubleCounting()
161 Int_t ID=kSpUndefined;
163 //get the PID response
165 AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
166 AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
167 fPIDResponse = inputHandler->GetPIDResponse();
170 AliFatal("Cannot get pid response");
173 //calculate nsigmas (used also by the bayesian)
174 CalculateNSigmas(trk,FIllQAHistos);//fill the data member fnsigmas with the nsigmas value [ipart][iPID]
177 if(fPIDType==kBayes){//use bayesianPID
180 AliFatal("PIDCombined object has to be set in the steering macro");
183 ID = GetIDBayes(trk,FIllQAHistos);
185 }else{ //use nsigma PID
187 ID=FindMinNSigma(trk,FIllQAHistos);
189 if(fUseExclusiveNSigma){ //if one particle has double counting and exclusive nsigma is requested ID = kSpUndefined
191 HasDC=GetDoubleCounting(trk,FIllQAHistos);
192 for(Int_t ipart=0;ipart<kNSpecies;ipart++){
193 if(HasDC[ipart]==kTRUE) ID = kSpUndefined;
198 //Fill PID signal plot
199 if(ID != kSpUndefined){
200 for(Int_t idet=0;idet<kNDetectors;idet++){
201 TH2F *h=GetHistogram2D(Form("PID_%d_%d",idet,ID));
202 if(idet==kITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
203 if(idet==kTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
204 if(idet==kTOF && fHasTOFPID)h->Fill(trk->P(),TOFBetaCalc(trk)*trk->Charge());
207 //Fill PID signal plot without cuts
208 for(Int_t idet=0;idet<kNDetectors;idet++){
209 TH2F *h=GetHistogram2D(Form("PIDAll_%d",idet));
210 if(idet==kITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
211 if(idet==kTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
212 if(idet==kTOF && fHasTOFPID)h->Fill(trk->P(),TOFBetaCalc(trk)*trk->Charge());
218 //////////////////////////////////////////////////////////////////////////////////////////////////
220 Int_t AliHelperPID::GetParticleSpecies(AliVParticle * part) {
221 // return PID according to MC truth
222 switch(TMath::Abs(part->PdgCode())){
237 //////////////////////////////////////////////////////////////////////////////////////////////////
239 Int_t AliHelperPID::GetIDBayes(AliVTrack * trk, Bool_t FIllQAHistos){
241 Bool_t *IDs=GetAllCompatibleIdentitiesNSigma(trk,FIllQAHistos);
243 Double_t probBayes[AliPID::kSPECIES];
247 if(fHasTOFPID && trk->Pt()>fPtTOFPID){//use TOF information
248 detUsed = CalcPIDCombined(trk, fPIDResponse, AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC, probBayes);
249 if(detUsed != (AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC))return kSpUndefined;//check that TPC and TOF are used
251 detUsed = CalcPIDCombined(trk, fPIDResponse,AliPIDResponse::kDetTPC, probBayes);
252 if(detUsed != AliPIDResponse::kDetTPC)return kSpUndefined;//check that TPC is used
255 //the probability has to be normalized to one, we check it
257 for(Int_t ipart=0;ipart<AliPID::kSPECIES;ipart++)sump+=probBayes[ipart];
258 if(sump<.99 && sump>1.01){//FIXME precision problem in the sum, workaround
259 AliFatal("Bayesian probability not normalized to one");
262 //probabilities are normalized to one, if the cut is above .5 there is no problem
263 if(probBayes[AliPID::kPion]>fBayesCut && IDs[kSpPion]==1){
264 TH2F *h=GetHistogram2D(Form("BayesRec_%d",kSpPion));
265 h->Fill(trk->Pt(),probBayes[AliPID::kPion]);
268 else if(probBayes[AliPID::kKaon]>fBayesCut && IDs[kSpKaon]==1){
269 TH2F *h=GetHistogram2D(Form("BayesRec_%d",kSpKaon));
270 h->Fill(trk->Pt(),probBayes[AliPID::kKaon]);
273 else if(probBayes[AliPID::kProton]>fBayesCut && IDs[kSpProton]==1){
274 TH2F *h=GetHistogram2D(Form("BayesRec_%d",kSpProton));
275 h->Fill(trk->Pt(),probBayes[AliPID::kProton]);
283 //////////////////////////////////////////////////////////////////////////////////////////////////
285 UInt_t AliHelperPID::CalcPIDCombined(const AliVTrack *track,const AliPIDResponse *PIDResponse, Int_t detMask, Double_t* prob) const{
287 // Bayesian PID calculation
289 for(Int_t i=0;i<AliPID::kSPECIES;i++)
293 fPIDCombined->SetDetectorMask(detMask);
295 return fPIDCombined->ComputeProbabilities((AliVTrack*)track, PIDResponse, prob);
298 //////////////////////////////////////////////////////////////////////////////////////////////////
300 void AliHelperPID::CalculateNSigmas(AliVTrack * trk, Bool_t FIllQAHistos){
301 //defines data member fnsigmas
303 // Compute nsigma for each hypthesis
304 AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(trk);
306 Double_t nsigmaTPCkProton = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton);
307 Double_t nsigmaTPCkKaon = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon);
308 Double_t nsigmaTPCkPion = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion);
310 Double_t nsigmaTOFkProton=999.,nsigmaTOFkKaon=999.,nsigmaTOFkPion=999.;
311 Double_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;
315 if(fHasTOFPID && trk->Pt()>fPtTOFPID){//use TOF information
316 nsigmaTOFkProton = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton);
317 nsigmaTOFkKaon = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon);
318 nsigmaTOFkPion = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion);
319 Double_t d2Proton=nsigmaTPCkProton * nsigmaTPCkProton + nsigmaTOFkProton * nsigmaTOFkProton;
320 Double_t d2Kaon=nsigmaTPCkKaon * nsigmaTPCkKaon + nsigmaTOFkKaon * nsigmaTOFkKaon;
321 Double_t d2Pion=nsigmaTPCkPion * nsigmaTPCkPion + nsigmaTOFkPion * nsigmaTOFkPion;
322 //commented, this is done in analogy with AliESDTrackCuts, nsigma combind according to the probability
324 // -----------------------------------
325 // How to get to a n-sigma cut?
327 // The accumulated statistics from 0 to d is
329 // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
330 // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
331 // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
333 // work around precision problem
334 // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
335 //if(TMath::Exp(- d2Proton / 2) > 1e-15)nsigmaTPCTOFkProton = TMath::Sqrt(2)*TMath::ErfInverse(1 - TMath::Exp(- d2Proton / 2));
336 //if(TMath::Exp(- d2Kaon / 2) > 1e-15)nsigmaTPCTOFkKaon = TMath::Sqrt(2)*TMath::ErfInverse(1 - TMath::Exp(- d2Kaon / 2));
337 //if(TMath::Exp(- d2Pion / 2) > 1e-15)nsigmaTPCTOFkPion = TMath::Sqrt(2)*TMath::ErfInverse(1 - TMath::Exp(- d2Pion / 2));
339 //used for the 2PC PID paper (circular cut)
340 nsigmaTPCTOFkProton = TMath::Sqrt(d2Proton);
341 nsigmaTPCTOFkKaon = TMath::Sqrt(d2Kaon);
342 nsigmaTPCTOFkPion = TMath::Sqrt(d2Pion);
345 // if TOF is missing and below fPtTOFPID only the TPC information is used
346 nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
347 nsigmaTPCTOFkKaon = TMath::Abs(nsigmaTPCkKaon);
348 nsigmaTPCTOFkPion = TMath::Abs(nsigmaTPCkPion);
351 //set data member fnsigmas
352 fnsigmas[kSpPion][kNSigmaTPC]=nsigmaTPCkPion;
353 fnsigmas[kSpKaon][kNSigmaTPC]=nsigmaTPCkKaon;
354 fnsigmas[kSpProton][kNSigmaTPC]=nsigmaTPCkProton;
355 fnsigmas[kSpPion][kNSigmaTOF]=nsigmaTOFkPion;
356 fnsigmas[kSpKaon][kNSigmaTOF]=nsigmaTOFkKaon;
357 fnsigmas[kSpProton][kNSigmaTOF]=nsigmaTOFkProton;
358 fnsigmas[kSpPion][kNSigmaTPCTOF]=nsigmaTPCTOFkPion;
359 fnsigmas[kSpKaon][kNSigmaTPCTOF]=nsigmaTPCTOFkKaon;
360 fnsigmas[kSpProton][kNSigmaTPCTOF]=nsigmaTPCTOFkProton;
363 //Fill NSigma SeparationPlot
364 for(Int_t ipart=0;ipart<kNSpecies;ipart++){
365 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
366 if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && trk->Pt()<fPtTOFPID)continue;//not filling TOF and combined if no TOF PID
367 TH2F *h=GetHistogram2D(Form("NSigma_%d_%d",ipart,ipid));
368 h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
374 //////////////////////////////////////////////////////////////////////////////////////////////////
376 Int_t AliHelperPID::FindMinNSigma(AliVTrack * trk,Bool_t FillQAHistos){
379 if(fRequestTOFPID && (!fHasTOFPID) && trk->Pt()>fPtTOFPID)return kSpUndefined;
381 //get the identity of the particle with the minimum Nsigma
382 Double_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
385 nsigmaProton = TMath::Abs(fnsigmas[kSpProton][kNSigmaTPC]);
386 nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPC]) ;
387 nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTPC]) ;
390 nsigmaProton = TMath::Abs(fnsigmas[kSpProton][kNSigmaTOF]);
391 nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTOF]) ;
392 nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTOF]) ;
394 case kNSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
395 nsigmaProton = TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
396 nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF]) ;
397 nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF]) ;
399 case kBayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
400 nsigmaProton = TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
401 nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF]) ;
402 nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF]) ;
406 // guess the particle based on the smaller nsigma (within fNSigmaPID)
407 if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return kSpUndefined;//if is the default value for the three
409 if( ( nsigmaKaon < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon < fNSigmaPID)){
411 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
412 if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
413 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",kSpKaon,ipid));
414 h->Fill(trk->Pt(),fnsigmas[kSpKaon][ipid]);
419 if( ( nsigmaPion < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion < fNSigmaPID)){
421 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
422 if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
423 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",kSpPion,ipid));
424 h->Fill(trk->Pt(),fnsigmas[kSpPion][ipid]);
429 if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)){
431 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
432 if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
433 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",kSpProton,ipid));
434 h->Fill(trk->Pt(),fnsigmas[kSpProton][ipid]);
439 // else, return undefined
443 //////////////////////////////////////////////////////////////////////////////////////////////////
445 Bool_t* AliHelperPID::GetDoubleCounting(AliVTrack * trk,Bool_t FIllQAHistos){
446 //if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE
448 for(Int_t ipart=0;ipart<kNSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track
450 Int_t MinNSigma=FindMinNSigma(trk,kFALSE);//not filling the NSigmaRec histos
454 if(MinNSigma==kSpUndefined)return fHasDoubleCounting;//in case of undefined no Double counting
456 Double_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
459 nsigmaProton = TMath::Abs(fnsigmas[kSpProton][kNSigmaTPC]);
460 nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPC]) ;
461 nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTPC]) ;
464 nsigmaProton = TMath::Abs(fnsigmas[kSpProton][kNSigmaTOF]);
465 nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTOF]) ;
466 nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTOF]) ;
468 case kNSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
469 nsigmaProton = TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
470 nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF]) ;
471 nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF]) ;
473 case kBayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
474 nsigmaProton = TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
475 nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF]) ;
476 nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF]) ;
479 if(nsigmaPion<fNSigmaPID && MinNSigma!=kSpPion)fHasDoubleCounting[kSpPion]=kTRUE;
480 if(nsigmaKaon<fNSigmaPID && MinNSigma!=kSpKaon)fHasDoubleCounting[kSpKaon]=kTRUE;
481 if(nsigmaProton<fNSigmaPID && MinNSigma!=kSpProton)fHasDoubleCounting[kSpProton]=kTRUE;
484 //fill NSigma distr for double counting
485 for(Int_t ipart=0;ipart<kNSpecies;ipart++){
486 if(fHasDoubleCounting[ipart]){
487 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
488 if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
489 TH2F *h=GetHistogram2D(Form("NSigmaDC_%d_%d",ipart,ipid));
490 h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
495 return fHasDoubleCounting;
498 //////////////////////////////////////////////////////////////////////////////////////////////////
500 Bool_t* AliHelperPID::GetAllCompatibleIdentitiesNSigma(AliVTrack * trk,Bool_t FIllQAHistos){
502 Bool_t *IDs=GetDoubleCounting(trk,FIllQAHistos);
503 IDs[FindMinNSigma(trk,FIllQAHistos)]=kTRUE;
508 //////////////////////////////////////////////////////////////////////////////////////////////////
510 Int_t AliHelperPID::GetMCParticleSpecie(AliVEvent* event, AliVTrack * trk, Bool_t FillQAHistos){
511 //return the specie according to the MC truth
514 if(!fisMC)AliFatal("Error: AliHelperPID::GetMCParticleSpecie called on data\n");
519 TString nameoftrack(event->ClassName());
520 if(!nameoftrack.CompareTo("AliESDEvent"))isESD=kTRUE;
521 else if(!nameoftrack.CompareTo("AliAODEvent"))isAOD=kTRUE;
522 else AliFatal("Not processing AODs nor ESDs") ;
525 TClonesArray *arrayMC = 0;
526 arrayMC = (TClonesArray*) event->GetList()->FindObject(AliAODMCParticle::StdBranchName());
527 if (!arrayMC)AliFatal("Error: MC particles branch not found!\n");
528 AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(trk->GetLabel()));
530 AliError("Cannot get MC particle");
533 code=partMC->GetPdgCode();
537 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
538 if (!eventHandler) AliFatal("ERROR: Could not retrieve MC event handler");
539 AliMCEvent* mcEvent = eventHandler->MCEvent();
540 if (!mcEvent)AliFatal("ERROR: Could not retrieve MC event");
541 stack = mcEvent->Stack();
542 if (!stack) AliFatal("ERROR: stack not available\n");
544 part = (TParticle*)stack->Particle(TMath::Abs(trk->GetLabel()));
546 AliError("Cannot get MC particle");
549 code = part->GetPdgCode();
551 switch(TMath::Abs(code)){
554 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
555 if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
556 TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",kSpProton,ipid));
557 h->Fill(trk->Pt(),fnsigmas[kSpProton][ipid]);
564 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
565 if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
566 TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",kSpKaon,ipid));
567 h->Fill(trk->Pt(),fnsigmas[kSpKaon][ipid]);
574 for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
575 if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
576 TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",kSpPion,ipid));
577 h->Fill(trk->Pt(),fnsigmas[kSpPion][ipid]);
588 //////////////////////////////////////////////////////////////////////////////////////////////////
590 void AliHelperPID::CheckTOF(AliVTrack * trk)
592 //check if the particle has TOF Matching
594 //get the PIDResponse
595 if(fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF,trk)==0)fHasTOFPID=kFALSE;
596 else fHasTOFPID=kTRUE;
598 //in addition to TOF status we look at the pt
599 if(trk->Pt()<fPtTOFPID)fHasTOFPID=kFALSE;
601 if(fRemoveTracksT0Fill)
603 Int_t startTimeMask = fPIDResponse->GetTOFResponse().GetStartTimeMask(trk->P());
604 if (startTimeMask < 0)fHasTOFPID=kFALSE;
608 //////////////////////////////////////////////////////////////////////////////////////////////////
610 Double_t AliHelperPID::TOFBetaCalc(AliVTrack *track) const{
611 //TOF beta calculation
612 Double_t tofTime=track->GetTOFsignal();
614 Double_t c=TMath::C()*1.E-9;// m/ns
615 Float_t startTime = fPIDResponse->GetTOFResponse().GetStartTime(((AliVTrack*)track)->P());//in ps
616 Double_t length= fPIDResponse->GetTOFResponse().GetExpectedSignal(track,AliPID::kElectron)*1E-3*c;
617 tofTime -= startTime; // subtract startTime to the signal
618 Double_t tof= tofTime*1E-3; // ns, average T0 fill subtracted, no info from T0detector
623 //////////////////////////////////////////////////////////////////////////////////////////////////
625 Double_t AliHelperPID::GetMass(AliHelperParticleSpecies_t id) const{
626 //return Mass according to AliHelperParticleSpecies_t. If undefined return -999.
628 if (id == kSpProton) { mass=9.38271999999999995e-01; }
629 if (id == kSpKaon) { mass=4.93676999999999977e-01; }
630 if (id == kSpPion) { mass=1.39570000000000000e-01; }
634 //////////////////////////////////////////////////////////////////////////////////////////////////
636 Long64_t AliHelperPID::Merge(TCollection* list)
638 // Merging interface.
639 // Returns the number of merged objects (including this).
645 TIterator* iter = list->MakeIterator();
647 // collections of TList with output histos
650 while ((obj = iter->Next())) {
651 AliHelperPID* entry = dynamic_cast<AliHelperPID*> (obj);
654 TList * outputlist = entry->GetOutputList();
655 collections.Add(outputlist);
658 fOutputList->Merge(&collections);