]>
Commit | Line | Data |
---|---|---|
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 | ||
38 | using namespace AliHelperPIDNameSpace; | |
39 | using namespace std; | |
40 | ||
41 | const char * kPIDTypeName[]={"TPC","TOF","TPC-TOF"} ; | |
d2062f6b | 42 | const char * kDetectorName[]={"ITS","TPC","TOF"} ; |
32e49607 | 43 | const char * kParticleSpeciesName[]={"Pions","Kaons","Protons","Undefined"} ; |
44 | ||
45 | ClassImp(AliHelperPID) | |
46 | ||
d2062f6b | 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){ |
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 | } | |
145 | ||
32e49607 | 146 | } |
147 | ||
148 | ////////////////////////////////////////////////////////////////////////////////////////////////// | |
149 | ||
150 | TH2F* AliHelperPID::GetHistogram2D(const char * name){ | |
151 | // returns histo named name | |
152 | return (TH2F*) fOutputList->FindObject(name); | |
153 | } | |
154 | ||
155 | ////////////////////////////////////////////////////////////////////////////////////////////////// | |
156 | ||
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() | |
32e49607 | 160 | |
d2062f6b | 161 | Int_t ID=kSpUndefined; |
162 | ||
163 | //get the PID response | |
164 | if(!fPIDResponse) { | |
165 | AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager(); | |
166 | AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler()); | |
167 | fPIDResponse = inputHandler->GetPIDResponse(); | |
168 | } | |
169 | if(!fPIDResponse) { | |
170 | AliFatal("Cannot get pid response"); | |
171 | } | |
172 | ||
173 | //calculate nsigmas (used also by the bayesian) | |
32e49607 | 174 | CalculateNSigmas(trk,FIllQAHistos);//fill the data member fnsigmas with the nsigmas value [ipart][iPID] |
175 | ||
d2062f6b | 176 | //Do PID |
177 | if(fPIDType==kBayes){//use bayesianPID | |
178 | ||
179 | if(!fPIDCombined) { | |
180 | AliFatal("PIDCombined object has to be set in the steering macro"); | |
181 | } | |
182 | ||
183 | ID = GetIDBayes(trk,FIllQAHistos); | |
184 | ||
185 | }else{ //use nsigma PID | |
186 | ||
187 | ID=FindMinNSigma(trk,FIllQAHistos); | |
188 | ||
189 | if(fUseExclusiveNSigma){ //if one particle has double counting and exclusive nsigma is requested ID = kSpUndefined | |
190 | Bool_t *HasDC; | |
191 | HasDC=GetDoubleCounting(trk,FIllQAHistos); | |
192 | for(Int_t ipart=0;ipart<kNSpecies;ipart++){ | |
193 | if(HasDC[ipart]==kTRUE) ID = kSpUndefined; | |
194 | } | |
5838bc91 | 195 | } |
196 | } | |
32e49607 | 197 | |
d2062f6b | 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()); | |
205 | } | |
206 | } | |
214084ec | 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()); | |
213 | } | |
214 | ||
d2062f6b | 215 | return ID; |
32e49607 | 216 | } |
217 | ||
218 | ////////////////////////////////////////////////////////////////////////////////////////////////// | |
219 | ||
7d3a3d9b | 220 | Int_t AliHelperPID::GetParticleSpecies(AliVParticle * part) { |
221 | // return PID according to MC truth | |
222 | switch(TMath::Abs(part->PdgCode())){ | |
223 | case 2212: | |
224 | return kSpProton; | |
225 | break; | |
226 | case 321: | |
227 | return kSpKaon; | |
228 | break; | |
229 | case 211: | |
230 | return kSpPion; | |
231 | break; | |
232 | default: | |
233 | return kSpUndefined; | |
234 | } | |
235 | } | |
236 | ||
237 | ////////////////////////////////////////////////////////////////////////////////////////////////// | |
238 | ||
d2062f6b | 239 | Int_t AliHelperPID::GetIDBayes(AliVTrack * trk, Bool_t FIllQAHistos){ |
240 | ||
241 | Bool_t *IDs=GetAllCompatibleIdentitiesNSigma(trk,FIllQAHistos); | |
242 | ||
243 | Double_t probBayes[AliPID::kSPECIES]; | |
244 | ||
245 | UInt_t detUsed= 0; | |
246 | CheckTOF(trk); | |
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 | |
250 | }else{ | |
251 | detUsed = CalcPIDCombined(trk, fPIDResponse,AliPIDResponse::kDetTPC, probBayes); | |
252 | if(detUsed != AliPIDResponse::kDetTPC)return kSpUndefined;//check that TPC is used | |
32e49607 | 253 | } |
d2062f6b | 254 | |
255 | //the probability has to be normalized to one, we check it | |
256 | Double_t sump=0.; | |
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"); | |
260 | } | |
261 | ||
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]); | |
266 | return kSpPion; | |
32e49607 | 267 | } |
d2062f6b | 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]); | |
271 | return kSpKaon; | |
272 | } | |
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]); | |
276 | return kSpProton; | |
277 | } | |
278 | else{ | |
279 | return kSpUndefined; | |
280 | } | |
281 | } | |
282 | ||
283 | ////////////////////////////////////////////////////////////////////////////////////////////////// | |
284 | ||
285 | UInt_t AliHelperPID::CalcPIDCombined(const AliVTrack *track,const AliPIDResponse *PIDResponse, Int_t detMask, Double_t* prob) const{ | |
286 | // | |
287 | // Bayesian PID calculation | |
288 | // | |
289 | for(Int_t i=0;i<AliPID::kSPECIES;i++) | |
290 | { | |
291 | prob[i]=0.; | |
292 | } | |
293 | fPIDCombined->SetDetectorMask(detMask); | |
294 | ||
295 | return fPIDCombined->ComputeProbabilities((AliVTrack*)track, PIDResponse, prob); | |
296 | } | |
297 | ||
298 | ////////////////////////////////////////////////////////////////////////////////////////////////// | |
299 | ||
300 | void AliHelperPID::CalculateNSigmas(AliVTrack * trk, Bool_t FIllQAHistos){ | |
301 | //defines data member fnsigmas | |
32e49607 | 302 | |
303 | // Compute nsigma for each hypthesis | |
304 | AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(trk); | |
305 | // --- TPC | |
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); | |
309 | // --- TOF | |
310 | Double_t nsigmaTOFkProton=999.,nsigmaTOFkKaon=999.,nsigmaTOFkPion=999.; | |
311 | Double_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.; | |
312 | ||
313 | CheckTOF(trk); | |
314 | ||
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); | |
aed2199f | 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 | |
32e49607 | 323 | // --- combined |
aed2199f | 324 | // ----------------------------------- |
325 | // How to get to a n-sigma cut? | |
326 | // | |
327 | // The accumulated statistics from 0 to d is | |
328 | // | |
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) | |
332 | // | |
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)); | |
338 | ||
339 | //used for the 2PC PID paper (circular cut) | |
340 | nsigmaTPCTOFkProton = TMath::Sqrt(d2Proton); | |
341 | nsigmaTPCTOFkKaon = TMath::Sqrt(d2Kaon); | |
342 | nsigmaTPCTOFkPion = TMath::Sqrt(d2Pion); | |
32e49607 | 343 | }else{ |
344 | // --- combined | |
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); | |
349 | } | |
350 | ||
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; | |
361 | ||
362 | if(FIllQAHistos){ | |
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]); | |
369 | } | |
370 | } | |
32e49607 | 371 | } |
372 | } | |
373 | ||
374 | ////////////////////////////////////////////////////////////////////////////////////////////////// | |
375 | ||
376 | Int_t AliHelperPID::FindMinNSigma(AliVTrack * trk,Bool_t FillQAHistos){ | |
377 | ||
378 | CheckTOF(trk); | |
379 | if(fRequestTOFPID && (!fHasTOFPID) && trk->Pt()>fPtTOFPID)return kSpUndefined; | |
380 | ||
381 | //get the identity of the particle with the minimum Nsigma | |
382 | Double_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.; | |
383 | switch (fPIDType){ | |
384 | case kNSigmaTPC: | |
385 | nsigmaProton = TMath::Abs(fnsigmas[kSpProton][kNSigmaTPC]); | |
386 | nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPC]) ; | |
387 | nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTPC]) ; | |
388 | break; | |
389 | case kNSigmaTOF: | |
390 | nsigmaProton = TMath::Abs(fnsigmas[kSpProton][kNSigmaTOF]); | |
391 | nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTOF]) ; | |
392 | nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTOF]) ; | |
393 | break; | |
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]) ; | |
398 | break; | |
d2062f6b | 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]) ; | |
403 | break; | |
32e49607 | 404 | } |
405 | ||
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 | |
408 | ||
409 | if( ( nsigmaKaon < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon < fNSigmaPID)){ | |
410 | if(FillQAHistos){ | |
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]); | |
415 | } | |
416 | } | |
417 | return kSpKaon; | |
418 | } | |
419 | if( ( nsigmaPion < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion < fNSigmaPID)){ | |
420 | if(FillQAHistos){ | |
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]); | |
425 | } | |
426 | } | |
427 | return kSpPion; | |
428 | } | |
429 | if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)){ | |
430 | if(FillQAHistos){ | |
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]); | |
435 | } | |
436 | } | |
437 | return kSpProton; | |
438 | } | |
439 | // else, return undefined | |
440 | return kSpUndefined; | |
441 | } | |
442 | ||
443 | ////////////////////////////////////////////////////////////////////////////////////////////////// | |
444 | ||
445 | Bool_t* AliHelperPID::GetDoubleCounting(AliVTrack * trk,Bool_t FIllQAHistos){ | |
446 | //if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE | |
447 | //fill DC histos | |
8ccf185c | 448 | for(Int_t ipart=0;ipart<kNSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track |
32e49607 | 449 | |
450 | Int_t MinNSigma=FindMinNSigma(trk,kFALSE);//not filling the NSigmaRec histos | |
451 | ||
452 | CheckTOF(trk); | |
453 | ||
454 | if(MinNSigma==kSpUndefined)return fHasDoubleCounting;//in case of undefined no Double counting | |
455 | ||
456 | Double_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.; | |
457 | switch (fPIDType) { | |
458 | case kNSigmaTPC: | |
459 | nsigmaProton = TMath::Abs(fnsigmas[kSpProton][kNSigmaTPC]); | |
460 | nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPC]) ; | |
461 | nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTPC]) ; | |
462 | break; | |
463 | case kNSigmaTOF: | |
464 | nsigmaProton = TMath::Abs(fnsigmas[kSpProton][kNSigmaTOF]); | |
465 | nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTOF]) ; | |
466 | nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTOF]) ; | |
467 | break; | |
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]) ; | |
472 | break; | |
d2062f6b | 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]) ; | |
477 | break; | |
32e49607 | 478 | } |
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; | |
482 | ||
483 | if(FIllQAHistos){ | |
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]); | |
491 | } | |
492 | } | |
493 | } | |
494 | } | |
495 | return fHasDoubleCounting; | |
496 | } | |
497 | ||
498 | ////////////////////////////////////////////////////////////////////////////////////////////////// | |
499 | ||
d2062f6b | 500 | Bool_t* AliHelperPID::GetAllCompatibleIdentitiesNSigma(AliVTrack * trk,Bool_t FIllQAHistos){ |
501 | ||
502 | Bool_t *IDs=GetDoubleCounting(trk,FIllQAHistos); | |
503 | IDs[FindMinNSigma(trk,FIllQAHistos)]=kTRUE; | |
504 | return IDs; | |
505 | ||
506 | } | |
507 | ||
508 | ////////////////////////////////////////////////////////////////////////////////////////////////// | |
509 | ||
32e49607 | 510 | Int_t AliHelperPID::GetMCParticleSpecie(AliVEvent* event, AliVTrack * trk, Bool_t FillQAHistos){ |
511 | //return the specie according to the MC truth | |
512 | CheckTOF(trk); | |
513 | ||
514 | if(!fisMC)AliFatal("Error: AliHelperPID::GetMCParticleSpecie called on data\n"); | |
515 | ||
516 | Int_t code=999; | |
517 | Bool_t isAOD=kFALSE; | |
518 | Bool_t isESD=kFALSE; | |
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") ; | |
523 | ||
524 | if(isAOD){ | |
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())); | |
529 | if (!partMC){ | |
530 | AliError("Cannot get MC particle"); | |
531 | return kSpUndefined; | |
532 | } | |
533 | code=partMC->GetPdgCode(); | |
534 | } | |
535 | if(isESD){ | |
536 | AliStack* stack=0x0; | |
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"); | |
543 | TParticle *part=0; | |
544 | part = (TParticle*)stack->Particle(TMath::Abs(trk->GetLabel())); | |
545 | if (!part){ | |
546 | AliError("Cannot get MC particle"); | |
547 | return kSpUndefined; | |
548 | } | |
549 | code = part->GetPdgCode(); | |
550 | } | |
551 | switch(TMath::Abs(code)){ | |
552 | case 2212: | |
553 | if(FillQAHistos){ | |
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]); | |
558 | } | |
559 | } | |
560 | return kSpProton; | |
561 | break; | |
562 | case 321: | |
563 | if(FillQAHistos){ | |
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]); | |
568 | } | |
569 | } | |
570 | return kSpKaon; | |
571 | break; | |
572 | case 211: | |
573 | if(FillQAHistos){ | |
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]); | |
578 | } | |
579 | } | |
580 | return kSpPion; | |
581 | break; | |
582 | default: | |
583 | return kSpUndefined; | |
584 | } | |
585 | ||
586 | } | |
587 | ||
588 | ////////////////////////////////////////////////////////////////////////////////////////////////// | |
589 | ||
590 | void AliHelperPID::CheckTOF(AliVTrack * trk) | |
591 | { | |
592 | //check if the particle has TOF Matching | |
d2062f6b | 593 | |
594 | //get the PIDResponse | |
595 | if(fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF,trk)==0)fHasTOFPID=kFALSE; | |
32e49607 | 596 | else fHasTOFPID=kTRUE; |
d2062f6b | 597 | |
598 | //in addition to TOF status we look at the pt | |
32e49607 | 599 | if(trk->Pt()<fPtTOFPID)fHasTOFPID=kFALSE; |
b46d749a | 600 | |
601 | if(fRemoveTracksT0Fill) | |
602 | { | |
b46d749a | 603 | Int_t startTimeMask = fPIDResponse->GetTOFResponse().GetStartTimeMask(trk->P()); |
604 | if (startTimeMask < 0)fHasTOFPID=kFALSE; | |
605 | } | |
32e49607 | 606 | } |
607 | ||
608 | ////////////////////////////////////////////////////////////////////////////////////////////////// | |
609 | ||
d2062f6b | 610 | Double_t AliHelperPID::TOFBetaCalc(AliVTrack *track) const{ |
611 | //TOF beta calculation | |
612 | Double_t tofTime=track->GetTOFsignal(); | |
613 | ||
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 | |
619 | tof=tof*c; | |
620 | return length/tof; | |
621 | } | |
622 | ||
623 | ////////////////////////////////////////////////////////////////////////////////////////////////// | |
624 | ||
7d3a3d9b | 625 | Double_t AliHelperPID::GetMass(AliHelperParticleSpecies_t id) const{ |
626 | //return Mass according to AliHelperParticleSpecies_t. If undefined return -999. | |
627 | Double_t mass=-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; } | |
631 | return mass; | |
632 | } | |
633 | ||
634 | ////////////////////////////////////////////////////////////////////////////////////////////////// | |
635 | ||
32e49607 | 636 | Long64_t AliHelperPID::Merge(TCollection* list) |
637 | { | |
638 | // Merging interface. | |
639 | // Returns the number of merged objects (including this). | |
640 | Printf("Merging"); | |
641 | if (!list) | |
642 | return 0; | |
643 | if (list->IsEmpty()) | |
644 | return 1; | |
645 | TIterator* iter = list->MakeIterator(); | |
646 | TObject* obj; | |
647 | // collections of TList with output histos | |
648 | TList collections; | |
649 | Int_t count = 0; | |
650 | while ((obj = iter->Next())) { | |
651 | AliHelperPID* entry = dynamic_cast<AliHelperPID*> (obj); | |
652 | if (entry == 0) | |
653 | continue; | |
654 | TList * outputlist = entry->GetOutputList(); | |
655 | collections.Add(outputlist); | |
656 | count++; | |
657 | } | |
658 | fOutputList->Merge(&collections); | |
659 | delete iter; | |
660 | Printf("OK"); | |
661 | return count+1; | |
662 | } |