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