ac72fd97a02a28f133d99302b8c1b7f0b90b0b7b
[u/mrichter/AliRoot.git] / PWG / Tools / AliHelperPID.cxx
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"   
34 #include "AliPIDCombined.h"   
35 #include "AliAnalysisManager.h"
36 #include "AliInputEventHandler.h"
37
38 using namespace AliHelperPIDNameSpace;
39 using namespace std;
40
41 const char * kPIDTypeName[]={"TPC","TOF","TPC-TOF"} ;
42 const char * kDetectorName[]={"ITS","TPC","TOF"} ;
43 const char * kParticleSpeciesName[]={"Pions","Kaons","Protons","Undefined"} ;
44
45 ClassImp(AliHelperPID)
46
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){
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;}
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);
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),
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);
83     }
84   }
85   
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   
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),
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);
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),
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);
122     }
123   }
124   
125   //PID signal plot
126   for(Int_t idet=0;idet<kNDetectors;idet++){
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     }
135   }
136 }
137
138 //////////////////////////////////////////////////////////////////////////////////////////////////
139
140 TH2F* AliHelperPID::GetHistogram2D(const char * name){
141   // returns histo named name
142   return (TH2F*) fOutputList->FindObject(name);
143 }
144
145 //////////////////////////////////////////////////////////////////////////////////////////////////
146
147 Int_t AliHelperPID::GetParticleSpecies(AliVTrack * trk, Bool_t FIllQAHistos){ 
148   //return the specie according to the minimum nsigma value
149   //no double counting, this has to be evaluated using CheckDoubleCounting()
150   
151   Int_t ID=kSpUndefined;
152   
153   //get the PID response
154   if(!fPIDResponse) {
155     AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
156     AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
157     fPIDResponse = inputHandler->GetPIDResponse();
158   }
159   if(!fPIDResponse) {
160     AliFatal("Cannot get pid response");
161   }
162   
163   //calculate nsigmas (used also by the bayesian)
164   CalculateNSigmas(trk,FIllQAHistos);//fill the data member fnsigmas with the nsigmas value [ipart][iPID]
165   
166   //Do PID
167   if(fPIDType==kBayes){//use bayesianPID
168     
169     if(!fPIDCombined) {
170       AliFatal("PIDCombined object has to be set in the steering macro");
171     }
172     
173     ID = GetIDBayes(trk,FIllQAHistos);
174     
175   }else{ //use nsigma PID
176     
177     ID=FindMinNSigma(trk,FIllQAHistos);
178     
179     if(fUseExclusiveNSigma){ //if one particle has double counting and exclusive nsigma is requested ID = kSpUndefined
180       Bool_t *HasDC;
181       HasDC=GetDoubleCounting(trk,FIllQAHistos);
182       for(Int_t ipart=0;ipart<kNSpecies;ipart++){
183         if(HasDC[ipart]==kTRUE)  ID = kSpUndefined;
184       }
185     }
186   }
187   
188   //Fill PID signal plot
189   if(ID != kSpUndefined){
190     for(Int_t idet=0;idet<kNDetectors;idet++){
191       TH2F *h=GetHistogram2D(Form("PID_%d_%d",idet,ID));
192       if(idet==kITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
193       if(idet==kTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
194       if(idet==kTOF && fHasTOFPID)h->Fill(trk->P(),TOFBetaCalc(trk)*trk->Charge());
195     }
196   }
197   return ID;
198 }
199
200 //////////////////////////////////////////////////////////////////////////////////////////////////
201
202 Int_t AliHelperPID::GetIDBayes(AliVTrack * trk, Bool_t FIllQAHistos){ 
203   
204   Bool_t *IDs=GetAllCompatibleIdentitiesNSigma(trk,FIllQAHistos);
205   
206   Double_t probBayes[AliPID::kSPECIES];
207   
208   UInt_t detUsed= 0;
209   CheckTOF(trk);
210   if(fHasTOFPID && trk->Pt()>fPtTOFPID){//use TOF information
211     detUsed = CalcPIDCombined(trk, fPIDResponse, AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC, probBayes);
212     if(detUsed != (AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC))return kSpUndefined;//check that TPC and TOF are used
213   }else{
214     detUsed = CalcPIDCombined(trk, fPIDResponse,AliPIDResponse::kDetTPC, probBayes);
215     if(detUsed != AliPIDResponse::kDetTPC)return kSpUndefined;//check that TPC is used
216   }
217   
218   //the probability has to be normalized to one, we check it
219   Double_t sump=0.;
220   for(Int_t ipart=0;ipart<AliPID::kSPECIES;ipart++)sump+=probBayes[ipart];
221   if(sump<.99 && sump>1.01){//FIXME precision problem in the sum, workaround
222     AliFatal("Bayesian probability not normalized to one");
223   }
224   
225   //probabilities are normalized to one, if the cut is above .5 there is no problem
226   if(probBayes[AliPID::kPion]>fBayesCut && IDs[kSpPion]==1){
227     TH2F *h=GetHistogram2D(Form("BayesRec_%d",kSpPion));
228     h->Fill(trk->Pt(),probBayes[AliPID::kPion]);
229     return kSpPion;
230   }
231   else if(probBayes[AliPID::kKaon]>fBayesCut && IDs[kSpKaon]==1){
232     TH2F *h=GetHistogram2D(Form("BayesRec_%d",kSpKaon));
233     h->Fill(trk->Pt(),probBayes[AliPID::kKaon]);
234     return kSpKaon;
235   }
236   else if(probBayes[AliPID::kProton]>fBayesCut && IDs[kSpProton]==1){
237     TH2F *h=GetHistogram2D(Form("BayesRec_%d",kSpProton));
238     h->Fill(trk->Pt(),probBayes[AliPID::kProton]);
239     return kSpProton;
240   }
241   else{
242     return kSpUndefined;
243   }
244 }
245
246 //////////////////////////////////////////////////////////////////////////////////////////////////
247
248 UInt_t AliHelperPID::CalcPIDCombined(const AliVTrack *track,const AliPIDResponse *PIDResponse, Int_t detMask, Double_t* prob) const{
249   //
250   // Bayesian PID calculation
251   //
252   for(Int_t i=0;i<AliPID::kSPECIES;i++)
253     {
254       prob[i]=0.;
255     }
256   fPIDCombined->SetDetectorMask(detMask);
257   
258   return fPIDCombined->ComputeProbabilities((AliVTrack*)track, PIDResponse, prob);
259 }
260
261 //////////////////////////////////////////////////////////////////////////////////////////////////
262
263 void AliHelperPID::CalculateNSigmas(AliVTrack * trk, Bool_t FIllQAHistos){ 
264   //defines data member fnsigmas
265   
266   // Compute nsigma for each hypthesis
267   AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(trk);
268   // --- TPC
269   Double_t nsigmaTPCkProton = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton);
270   Double_t nsigmaTPCkKaon   = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon); 
271   Double_t nsigmaTPCkPion   = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion); 
272   // --- TOF
273   Double_t nsigmaTOFkProton=999.,nsigmaTOFkKaon=999.,nsigmaTOFkPion=999.;
274   Double_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;
275
276   CheckTOF(trk);
277   
278   if(fHasTOFPID && trk->Pt()>fPtTOFPID){//use TOF information
279     nsigmaTOFkProton = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton);
280     nsigmaTOFkKaon   = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon); 
281     nsigmaTOFkPion   = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion); 
282     Double_t d2Proton=nsigmaTPCkProton * nsigmaTPCkProton + nsigmaTOFkProton * nsigmaTOFkProton;
283     Double_t d2Kaon=nsigmaTPCkKaon * nsigmaTPCkKaon + nsigmaTOFkKaon * nsigmaTOFkKaon;
284     Double_t d2Pion=nsigmaTPCkPion * nsigmaTPCkPion + nsigmaTOFkPion * nsigmaTOFkPion;
285     //commented, this is done in analogy with AliESDTrackCuts, nsigma combind according to the probability
286     // --- combined
287     // -----------------------------------
288     // How to get to a n-sigma cut?
289     //
290     // The accumulated statistics from 0 to d is
291     //
292     // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
293     // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
294     // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
295     //
296     // work around precision problem
297     // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
298     //if(TMath::Exp(- d2Proton / 2) > 1e-15)nsigmaTPCTOFkProton  =  TMath::Sqrt(2)*TMath::ErfInverse(1 - TMath::Exp(- d2Proton / 2));
299     //if(TMath::Exp(- d2Kaon / 2) > 1e-15)nsigmaTPCTOFkKaon    = TMath::Sqrt(2)*TMath::ErfInverse(1 - TMath::Exp(- d2Kaon / 2));
300     //if(TMath::Exp(- d2Pion / 2) > 1e-15)nsigmaTPCTOFkPion     = TMath::Sqrt(2)*TMath::ErfInverse(1 - TMath::Exp(- d2Pion / 2));
301     
302     //used for the 2PC PID paper (circular cut)
303     nsigmaTPCTOFkProton  =  TMath::Sqrt(d2Proton);
304     nsigmaTPCTOFkKaon    =  TMath::Sqrt(d2Kaon);
305     nsigmaTPCTOFkPion    =  TMath::Sqrt(d2Pion);
306   }else{
307     // --- combined
308     // if TOF is missing and below fPtTOFPID only the TPC information is used
309     nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
310     nsigmaTPCTOFkKaon   = TMath::Abs(nsigmaTPCkKaon);
311     nsigmaTPCTOFkPion   = TMath::Abs(nsigmaTPCkPion);
312   }
313   
314   //set data member fnsigmas
315   fnsigmas[kSpPion][kNSigmaTPC]=nsigmaTPCkPion;
316   fnsigmas[kSpKaon][kNSigmaTPC]=nsigmaTPCkKaon;
317   fnsigmas[kSpProton][kNSigmaTPC]=nsigmaTPCkProton;
318   fnsigmas[kSpPion][kNSigmaTOF]=nsigmaTOFkPion;
319   fnsigmas[kSpKaon][kNSigmaTOF]=nsigmaTOFkKaon;
320   fnsigmas[kSpProton][kNSigmaTOF]=nsigmaTOFkProton;
321   fnsigmas[kSpPion][kNSigmaTPCTOF]=nsigmaTPCTOFkPion;
322   fnsigmas[kSpKaon][kNSigmaTPCTOF]=nsigmaTPCTOFkKaon;
323   fnsigmas[kSpProton][kNSigmaTPCTOF]=nsigmaTPCTOFkProton;
324   
325   if(FIllQAHistos){
326     //Fill NSigma SeparationPlot
327     for(Int_t ipart=0;ipart<kNSpecies;ipart++){
328       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
329         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && trk->Pt()<fPtTOFPID)continue;//not filling TOF and combined if no TOF PID
330         TH2F *h=GetHistogram2D(Form("NSigma_%d_%d",ipart,ipid));
331         h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
332       }
333     }
334   }
335 }
336
337 //////////////////////////////////////////////////////////////////////////////////////////////////
338
339 Int_t AliHelperPID::FindMinNSigma(AliVTrack * trk,Bool_t FillQAHistos){ 
340   
341   CheckTOF(trk);  
342   if(fRequestTOFPID && (!fHasTOFPID) && trk->Pt()>fPtTOFPID)return kSpUndefined;
343   
344   //get the identity of the particle with the minimum Nsigma
345   Double_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
346   switch (fPIDType){
347   case kNSigmaTPC:
348     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPC]);
349     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPC])  ;
350     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPC])  ;
351     break;
352   case kNSigmaTOF:
353     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTOF]);
354     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTOF])  ;
355     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTOF])  ;
356     break;
357   case kNSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
358     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
359     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF])  ;
360     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF])  ;
361     break;
362   case kBayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
363     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
364     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF])  ;
365     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF])  ;
366     break;
367   }
368   
369   // guess the particle based on the smaller nsigma (within fNSigmaPID)
370   if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return kSpUndefined;//if is the default value for the three
371   
372   if( ( nsigmaKaon   < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon   < fNSigmaPID)){
373     if(FillQAHistos){
374       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
375         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
376         TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",kSpKaon,ipid));
377         h->Fill(trk->Pt(),fnsigmas[kSpKaon][ipid]);
378       }
379     }
380     return kSpKaon;
381   }
382   if( ( nsigmaPion   < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion   < fNSigmaPID)){
383     if(FillQAHistos){
384       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
385         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
386         TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",kSpPion,ipid));
387         h->Fill(trk->Pt(),fnsigmas[kSpPion][ipid]);
388       }
389     }
390     return kSpPion;
391   }
392   if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)){
393     if(FillQAHistos){
394       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
395         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
396         TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",kSpProton,ipid));
397         h->Fill(trk->Pt(),fnsigmas[kSpProton][ipid]);
398       }
399     }
400     return kSpProton;
401   }
402   // else, return undefined
403   return kSpUndefined;
404 }
405
406 //////////////////////////////////////////////////////////////////////////////////////////////////
407
408 Bool_t* AliHelperPID::GetDoubleCounting(AliVTrack * trk,Bool_t FIllQAHistos){ 
409   //if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE
410   //fill DC histos
411   for(Int_t ipart=0;ipart<=kNSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track
412   
413   Int_t MinNSigma=FindMinNSigma(trk,kFALSE);//not filling the NSigmaRec histos
414   
415   CheckTOF(trk);
416   
417   if(MinNSigma==kSpUndefined)return fHasDoubleCounting;//in case of undefined no Double counting
418   
419   Double_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
420   switch (fPIDType) {
421   case kNSigmaTPC:
422     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPC]);
423     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPC])  ;
424     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPC])  ;
425     break;
426   case kNSigmaTOF:
427     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTOF]);
428     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTOF])  ;
429     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTOF])  ;
430     break;
431   case kNSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
432     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
433     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF])  ;
434     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF])  ;
435     break;
436   case kBayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
437     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
438     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF])  ;
439     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF])  ;
440     break;
441   }
442   if(nsigmaPion<fNSigmaPID && MinNSigma!=kSpPion)fHasDoubleCounting[kSpPion]=kTRUE;
443   if(nsigmaKaon<fNSigmaPID && MinNSigma!=kSpKaon)fHasDoubleCounting[kSpKaon]=kTRUE;
444   if(nsigmaProton<fNSigmaPID && MinNSigma!=kSpProton)fHasDoubleCounting[kSpProton]=kTRUE;
445   
446   if(FIllQAHistos){
447     //fill NSigma distr for double counting
448     for(Int_t ipart=0;ipart<kNSpecies;ipart++){
449       if(fHasDoubleCounting[ipart]){
450         for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
451           if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
452           TH2F *h=GetHistogram2D(Form("NSigmaDC_%d_%d",ipart,ipid));
453           h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
454         }
455       }
456     }
457   }
458   return fHasDoubleCounting;
459 }
460
461 //////////////////////////////////////////////////////////////////////////////////////////////////
462
463 Bool_t* AliHelperPID::GetAllCompatibleIdentitiesNSigma(AliVTrack * trk,Bool_t FIllQAHistos){ 
464   
465   Bool_t *IDs=GetDoubleCounting(trk,FIllQAHistos);
466   IDs[FindMinNSigma(trk,FIllQAHistos)]=kTRUE;
467   return IDs;
468   
469 }
470
471 //////////////////////////////////////////////////////////////////////////////////////////////////
472
473 Int_t AliHelperPID::GetMCParticleSpecie(AliVEvent* event, AliVTrack * trk, Bool_t FillQAHistos){ 
474   //return the specie according to the MC truth
475   CheckTOF(trk);
476   
477   if(!fisMC)AliFatal("Error: AliHelperPID::GetMCParticleSpecie called on data\n");
478   
479   Int_t code=999;
480   Bool_t isAOD=kFALSE;
481   Bool_t isESD=kFALSE;
482   TString nameoftrack(event->ClassName());  
483   if(!nameoftrack.CompareTo("AliESDEvent"))isESD=kTRUE;
484   else if(!nameoftrack.CompareTo("AliAODEvent"))isAOD=kTRUE;
485   else AliFatal("Not processing AODs nor ESDs") ;
486   
487   if(isAOD){
488     TClonesArray *arrayMC = 0;
489     arrayMC = (TClonesArray*) event->GetList()->FindObject(AliAODMCParticle::StdBranchName());
490     if (!arrayMC)AliFatal("Error: MC particles branch not found!\n");
491     AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(trk->GetLabel()));
492     if (!partMC){
493       AliError("Cannot get MC particle");
494       return kSpUndefined;
495     }
496     code=partMC->GetPdgCode();
497   }
498   if(isESD){
499     AliStack* stack=0x0;
500     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
501     if (!eventHandler) AliFatal("ERROR: Could not retrieve MC event handler");
502     AliMCEvent* mcEvent = eventHandler->MCEvent();
503     if (!mcEvent)AliFatal("ERROR: Could not retrieve MC event");
504     stack = mcEvent->Stack();
505     if (!stack) AliFatal("ERROR: stack not available\n");
506     TParticle *part=0;
507     part = (TParticle*)stack->Particle(TMath::Abs(trk->GetLabel()));
508     if (!part){
509       AliError("Cannot get MC particle");
510       return kSpUndefined;
511     }
512     code = part->GetPdgCode();
513   }
514   switch(TMath::Abs(code)){
515   case 2212:
516     if(FillQAHistos){
517       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
518         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
519         TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",kSpProton,ipid));
520         h->Fill(trk->Pt(),fnsigmas[kSpProton][ipid]);
521       }
522     }
523     return kSpProton;
524     break;
525   case 321:
526     if(FillQAHistos){
527       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
528         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
529         TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",kSpKaon,ipid));
530         h->Fill(trk->Pt(),fnsigmas[kSpKaon][ipid]);
531       }
532     }
533     return kSpKaon;
534     break;
535   case 211:
536     if(FillQAHistos){
537       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
538         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
539         TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",kSpPion,ipid));
540         h->Fill(trk->Pt(),fnsigmas[kSpPion][ipid]);
541       }
542     }
543     return kSpPion;
544     break;
545   default:
546     return kSpUndefined;
547   }
548
549 }
550
551 //////////////////////////////////////////////////////////////////////////////////////////////////
552
553 void AliHelperPID::CheckTOF(AliVTrack * trk)
554 {
555   //check if the particle has TOF Matching
556   
557   //get the PIDResponse
558   if(fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF,trk)==0)fHasTOFPID=kFALSE;
559   else fHasTOFPID=kTRUE;
560   
561   //in addition to TOF status we look at the pt
562   if(trk->Pt()<fPtTOFPID)fHasTOFPID=kFALSE;
563   
564   if(fRemoveTracksT0Fill)
565     {
566       Int_t startTimeMask = fPIDResponse->GetTOFResponse().GetStartTimeMask(trk->P());
567       if (startTimeMask < 0)fHasTOFPID=kFALSE; 
568     }
569 }
570
571 //////////////////////////////////////////////////////////////////////////////////////////////////
572
573 Double_t AliHelperPID::TOFBetaCalc(AliVTrack *track) const{
574   //TOF beta calculation
575   Double_t tofTime=track->GetTOFsignal();
576   
577   Double_t c=TMath::C()*1.E-9;// m/ns
578   Float_t startTime = fPIDResponse->GetTOFResponse().GetStartTime(((AliVTrack*)track)->P());//in ps
579   Double_t length= fPIDResponse->GetTOFResponse().GetExpectedSignal(track,AliPID::kElectron)*1E-3*c;
580   tofTime -= startTime;      // subtract startTime to the signal
581   Double_t tof= tofTime*1E-3; // ns, average T0 fill subtracted, no info from T0detector         
582   tof=tof*c;
583   return length/tof;
584 }
585
586 //////////////////////////////////////////////////////////////////////////////////////////////////
587
588 Long64_t AliHelperPID::Merge(TCollection* list)
589 {
590   // Merging interface.
591   // Returns the number of merged objects (including this).
592   Printf("Merging");
593   if (!list)
594     return 0;
595   if (list->IsEmpty())
596     return 1;
597   TIterator* iter = list->MakeIterator();
598   TObject* obj;
599   // collections of TList with output histos
600   TList collections;
601   Int_t count = 0;
602   while ((obj = iter->Next())) {
603     AliHelperPID* entry = dynamic_cast<AliHelperPID*> (obj);
604     if (entry == 0) 
605       continue;
606     TList * outputlist = entry->GetOutputList();      
607     collections.Add(outputlist);
608     count++;
609   }
610   fOutputList->Merge(&collections);
611   delete iter;
612   Printf("OK");
613   return count+1;
614 }