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