]>
Commit | Line | Data |
---|---|---|
563113d0 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | //This class is intended to provide a selection on | |
17 | //the PID for single charged particles as electrons, muons | |
18 | //pions, kaons and protons. The user is supposed to set only one particle | |
19 | //to look at. The class at the moment uses one single ESD track. | |
20 | //The identification is done in Identify(), whereas the GetID() checks | |
21 | //the track status or if the combined PID should be applied. | |
22 | //The Correction Framework follows the Annalysis framework. | |
23 | //The main method of this class is the IsSelected function which returns | |
24 | //true (false) if the ESD track is (not) identified as the previously | |
25 | //setted particle. | |
26 | // | |
27 | // | |
28 | //usage: | |
29 | // | |
30 | // -----ID by one detector------ | |
31 | //AliCFTrackCutPid *pCut = new AliCFTrackCutPid("pion_sel","pion_sel"); | |
32 | //Double_t priors[5]={...}; | |
33 | //pCut->SetPartycleType(AliPID::kPion,kFALSE) | |
34 | //pCut->SetDetectors("DET"); ^^^^^^ | |
35 | //pCut->SetPriors(priors); | |
36 | // | |
37 | // -----ID combining more detectors---- | |
38 | //AliCFTrackCutPid *pCut = new AliCFTrackCutPid("pion_sel","pion_sel"); | |
39 | //Double_t priors[5]={...}; | |
40 | //pCut->SetPartycleType(AliPID::kPion,kTRUE) | |
41 | //pCut->SetDetectors("DET1 DET2 .."); ^^^^^ | |
42 | //pCut->SetPriors(priors) | |
43 | // | |
44 | //Comments: | |
45 | //The user can choose not to identify a track if one residual beteween | |
46 | //the identified particle probability and one among all the other | |
47 | //probabilties is smaller than a predefined value. | |
48 | //The same can be done for the detector responses. | |
49 | //This happens by setting: | |
50 | // | |
51 | //pCut->SetMinDiffProb(0.005,kTRUE) //for probabilities | |
52 | // | |
53 | //pCut->SetMinDiffResp(0.005,kTRUE) //for det responses | |
54 | // | |
55 | //Annalisa.Mastroserio@ba.infn.it | |
56 | // | |
57 | ||
58 | ||
f9658e11 | 59 | |
563113d0 | 60 | #include "AliCFTrackCutPid.h" |
61 | #include "AliLog.h" | |
62 | #include <TMath.h> | |
63 | #include <TList.h> | |
64 | ||
65 | ClassImp(AliCFTrackCutPid) | |
66 | ||
67 | //__________________________________________________________________________________ | |
68 | // CUT ON TRACK PID | |
69 | //__________________________________________________________________________________ | |
70 | AliCFTrackCutPid::AliCFTrackCutPid() : | |
71 | AliCFCutBase(), | |
72 | fCut(0.2), | |
73 | fMinDiffResponse(0.0001), | |
74 | fMinDiffProbability(0.001), | |
75 | fgParticleType(10), | |
76 | fgIsComb(kTRUE), | |
bc6176fb | 77 | fgIsAOD(kFALSE), |
563113d0 | 78 | fCheckResponse(kFALSE), |
79 | fCheckSelection(kTRUE), | |
80 | fIsPpriors(kFALSE), | |
81 | fIsDetAND(kFALSE), | |
82 | fXmin(-0.005), | |
83 | fXmax(1.005), | |
84 | fNbins(101), | |
85 | fDetRestr(-1), | |
86 | fiPartRestr(-1), | |
39f34d57 | 87 | fDetProbRestr(1), |
88 | fProbThreshold(0.) | |
563113d0 | 89 | { |
90 | // | |
91 | //Default constructor | |
92 | // | |
00a38d07 | 93 | for(Int_t j=0; j< AliPID::kSPECIES; j++) { |
563113d0 | 94 | fPriors[j]=0.2; |
ae8c412c | 95 | } |
96 | for(Int_t j=0; j< AliPID::kSPECIES; j++) { | |
563113d0 | 97 | fPriorsFunc[j]=0x0; |
98 | } | |
563113d0 | 99 | for(Int_t jDet=0; jDet< kNdets; jDet++) { |
100 | fDets[jDet]=kFALSE; | |
101 | fDetsInAnd[jDet]=kFALSE; | |
102 | } | |
103 | ||
104 | InitialiseHisto(); | |
105 | } | |
106 | //______________________________ | |
107 | AliCFTrackCutPid::AliCFTrackCutPid(const Char_t* name, const Char_t* title) : | |
108 | AliCFCutBase(name,title), | |
109 | fCut(0.2), | |
110 | fMinDiffResponse(0.0001), | |
111 | fMinDiffProbability(0.001), | |
112 | fgParticleType(10), | |
113 | fgIsComb(kTRUE), | |
bc6176fb | 114 | fgIsAOD(kFALSE), |
563113d0 | 115 | fCheckResponse(kFALSE), |
116 | fCheckSelection(kTRUE), | |
117 | fIsPpriors(kFALSE), | |
118 | fIsDetAND(kFALSE), | |
119 | fXmin(-0.005), | |
120 | fXmax(1.005), | |
121 | fNbins(101), | |
122 | fDetRestr(-1), | |
123 | fiPartRestr(-1), | |
39f34d57 | 124 | fDetProbRestr(1), |
125 | fProbThreshold(0.) | |
563113d0 | 126 | { |
127 | // | |
128 | //Constructor | |
129 | // | |
00a38d07 | 130 | for(Int_t j=0; j< AliPID::kSPECIES; j++) { |
563113d0 | 131 | fPriors[j]=0.2; |
ae8c412c | 132 | } |
133 | for(Int_t j=0; j< AliPID::kSPECIES; j++) { | |
563113d0 | 134 | fPriorsFunc[j]=0x0; |
135 | } | |
563113d0 | 136 | for(Int_t jDet=0; jDet< kNdets; jDet++) { |
137 | fDets[jDet]=kFALSE; | |
138 | fDetsInAnd[jDet]=kFALSE; | |
139 | } | |
140 | ||
141 | InitialiseHisto(); | |
142 | } | |
143 | //______________________________ | |
144 | AliCFTrackCutPid::AliCFTrackCutPid(const AliCFTrackCutPid& c) : | |
145 | AliCFCutBase(c), | |
146 | fCut(c.fCut), | |
147 | fMinDiffResponse(c.fMinDiffResponse), | |
148 | fMinDiffProbability(c.fMinDiffProbability), | |
149 | fgParticleType(c.fgParticleType), | |
150 | fgIsComb(c.fgIsComb), | |
bc6176fb | 151 | fgIsAOD(c.fgIsAOD), |
563113d0 | 152 | fCheckResponse(c.fCheckResponse), |
153 | fCheckSelection(c.fCheckSelection), | |
154 | fIsPpriors(c.fIsPpriors), | |
155 | fIsDetAND(c.fIsDetAND), | |
156 | fXmin(c.fXmin), | |
157 | fXmax(c.fXmax), | |
158 | fNbins(c.fNbins), | |
159 | fDetRestr(c.fDetRestr), | |
160 | fiPartRestr(c.fiPartRestr), | |
39f34d57 | 161 | fDetProbRestr(c.fDetProbRestr), |
162 | fProbThreshold(c.fProbThreshold) | |
563113d0 | 163 | { |
164 | // | |
165 | //Copy constructor | |
166 | // | |
167 | for(Int_t i=0; i< kNdets ; i++ ) { | |
168 | fDets[i]=c.fDets[i]; | |
169 | fDetsInAnd[i]=c.fDetsInAnd[i]; | |
170 | for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){ | |
171 | fhResp[i][iP]=c.fhResp[i][iP]; | |
172 | fhProb[i][iP]=c.fhProb[i][iP]; | |
173 | } | |
174 | } | |
00a38d07 | 175 | for(Int_t j=0; j< AliPID::kSPECIES; j++){ |
563113d0 | 176 | fPriors[j]=c.fPriors[j]; |
ae8c412c | 177 | } |
178 | for(Int_t j=0; j< AliPID::kSPECIES; j++){ | |
563113d0 | 179 | fPriorsFunc[j]=c.fPriorsFunc[j]; |
180 | fhCombResp[j]=c.fhCombResp[j]; | |
181 | fhCombProb[j]=c.fhCombProb[j]; | |
182 | } | |
183 | } | |
184 | //______________________________ | |
185 | AliCFTrackCutPid& AliCFTrackCutPid::operator=(const AliCFTrackCutPid& c) | |
186 | { | |
187 | // | |
188 | // Assignment operator | |
189 | // | |
190 | if (this != &c) { | |
191 | AliCFCutBase::operator=(c) ; | |
192 | this->fCut=c.fCut; | |
193 | this->fMinDiffResponse=c.fMinDiffResponse; | |
bc6176fb | 194 | this->fMinDiffProbability=c.fMinDiffProbability; |
563113d0 | 195 | this->fgParticleType=c.fgParticleType; |
196 | this->fgIsComb=c.fgIsComb; | |
bc6176fb | 197 | this->fgIsAOD=c.fgIsAOD; |
563113d0 | 198 | this->fCheckResponse=c.fCheckResponse; |
199 | this->fCheckSelection=c.fCheckSelection; | |
200 | this->fIsPpriors=c.fIsPpriors; | |
201 | this->fIsDetAND=c.fIsDetAND; | |
202 | this->fXmin=c.fXmin; | |
203 | this->fXmax=c.fXmax; | |
204 | this->fNbins=c.fNbins; | |
205 | this->fDetRestr=c.fDetRestr; | |
206 | this->fiPartRestr=c.fiPartRestr; | |
207 | this->fDetProbRestr=c.fDetProbRestr; | |
39f34d57 | 208 | this->fProbThreshold=c.fProbThreshold; |
563113d0 | 209 | |
210 | for(Int_t i=0; i< kNdets ; i++ ) { | |
211 | this->fDets[i]=c.fDets[i]; | |
212 | for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){ | |
213 | this->fhResp[i][iP]=c.fhResp[i][iP]; | |
214 | this->fhProb[i][iP]=c.fhProb[i][iP]; | |
215 | } | |
216 | } | |
217 | ||
00a38d07 | 218 | for(Int_t j=0; j< AliPID::kSPECIES; j++){ |
563113d0 | 219 | this->fPriors[j]=c.fPriors[j]; |
ae8c412c | 220 | } |
221 | for(Int_t j=0; j< AliPID::kSPECIES; j++){ | |
563113d0 | 222 | this->fhCombResp[j]=c.fhCombResp[j]; |
223 | this->fhCombProb[j]=c.fhCombProb[j]; | |
224 | this-> fPriorsFunc[j]=c.fPriorsFunc[j]; | |
225 | ||
226 | } | |
227 | } | |
228 | return *this ; | |
229 | } | |
230 | //__________________________________________________________________________________ | |
231 | AliCFTrackCutPid::~AliCFTrackCutPid() { | |
232 | // | |
233 | //dtor | |
234 | // | |
7f0c5c1e | 235 | |
563113d0 | 236 | for(Int_t i=0; i< kNdets ; i++ ) { |
237 | for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){ | |
238 | if(fhResp[i][iP])delete fhResp[i][iP]; | |
239 | if(fhProb[i][iP])delete fhProb[i][iP]; | |
240 | } | |
241 | } | |
242 | ||
243 | for(Int_t j=0; j< AliPID::kSPECIES; j++){ | |
244 | if(fhCombResp[j])delete fhCombResp[j]; | |
245 | if(fhCombProb[j])delete fhCombProb[j]; | |
246 | ||
247 | } | |
563113d0 | 248 | } |
249 | //__________________________________ | |
250 | void AliCFTrackCutPid::SetDetectors(TString dets) | |
251 | { | |
252 | // | |
253 | // The string of detectors is translated into | |
254 | // the respective booelan data members | |
255 | ||
256 | if(dets.Contains("ITS")) {fDets[kITS]=kTRUE;} | |
257 | if(dets.Contains("TPC")) {fDets[kTPC]=kTRUE;} | |
258 | if(dets.Contains("TRD")) {fDets[kTRD]=kTRUE;} | |
259 | if(dets.Contains("TOF")) {fDets[kTOF]=kTRUE;} | |
260 | if(dets.Contains("HMPID")) {fDets[kHMPID]=kTRUE;} | |
7f0c5c1e | 261 | |
262 | if(dets.Contains("ALL")) for(Int_t i=0; i< kNdets ; i++) fDets[i]=kTRUE; | |
563113d0 | 263 | } |
264 | //__________________________________ | |
265 | void AliCFTrackCutPid::SetPriors(Double_t r[AliPID::kSPECIES]) | |
266 | { | |
267 | // | |
268 | // Sets the a priori concentrations | |
269 | // | |
270 | for(Int_t i=0; i<AliPID::kSPECIES; i++) fPriors[i]=r[i]; | |
271 | } | |
272 | //__________________________________ | |
273 | void AliCFTrackCutPid::SetPriorFunctions(TF1 *func[AliPID::kSPECIES]) | |
274 | { | |
275 | // | |
276 | // Sets the momentu dependent a priori concentrations | |
277 | // | |
278 | ||
279 | for(Int_t i=0; i<AliPID::kSPECIES; i++) fPriorsFunc[i]=func[i]; | |
280 | fIsPpriors = kTRUE; | |
281 | } | |
282 | //____________________________________________ | |
283 | void AliCFTrackCutPid::SetANDstatus(TString dets) | |
284 | { | |
285 | // | |
286 | //Sets the detectors to be in AND for the combined PID | |
287 | // | |
288 | if(dets.Contains("ITS") && fDets[kITS]) {fDetsInAnd[kITS]=kTRUE;} | |
289 | if(dets.Contains("TPC") && fDets[kTPC]) {fDetsInAnd[kTPC]=kTRUE;} | |
290 | if(dets.Contains("TRD") && fDets[kTRD]) {fDetsInAnd[kTRD]=kTRUE;} | |
291 | if(dets.Contains("TOF") && fDets[kTOF]) {fDetsInAnd[kTOF]=kTRUE;} | |
292 | if(dets.Contains("HMPID") && fDets[kHMPID]) {fDetsInAnd[kHMPID]=kTRUE;} | |
293 | ||
294 | fIsDetAND = kTRUE; | |
295 | } | |
296 | // | |
297 | void AliCFTrackCutPid::SetDetectorProbabilityRestriction(TString det, Int_t iPart, Double_t upperprob) | |
298 | { | |
299 | // | |
300 | // Sets the detector, the particle and the probability | |
301 | // limit. | |
302 | // | |
303 | ||
304 | if(det.Contains("ITS")) fDetRestr = kITS; | |
305 | if(det.Contains("TPC")) fDetRestr = kTPC; | |
306 | if(det.Contains("TRD")) fDetRestr = kTRD; | |
307 | if(det.Contains("TOF")) fDetRestr = kTOF; | |
308 | if(det.Contains("HMPID")) fDetRestr = kHMPID; | |
309 | fiPartRestr = iPart; | |
310 | fDetProbRestr = upperprob; | |
311 | } | |
312 | //__________________________________ | |
313 | void AliCFTrackCutPid::TrackInfo(const AliESDtrack *pTrk, ULong_t status[kNdets+1],Double_t pid[kNdets+1][AliPID::kSPECIES]) const | |
314 | { | |
315 | // | |
316 | // Loads the responses of the XXX chosen detectors for the pTrk | |
317 | // and the corresponding trk status. The final trk status is also loaded. | |
318 | // | |
319 | if(fDets[kITS]) { | |
320 | pTrk->GetITSpid(pid[kITS]); | |
321 | status[kITS]=AliESDtrack::kITSpid; | |
322 | } | |
323 | if(fDets[kTPC]) { | |
324 | pTrk->GetTPCpid(pid[kTPC]); | |
325 | status[kTPC]=AliESDtrack::kTPCpid; | |
326 | } | |
327 | if(fDets[kTRD]) { | |
328 | pTrk->GetTRDpid(pid[kTRD]); | |
329 | status[kTRD]=AliESDtrack::kTRDpid; | |
330 | } | |
331 | if(fDets[kTOF]) { | |
332 | pTrk->GetTOFpid(pid[kTOF]); | |
333 | status[kTOF]=AliESDtrack::kTOFpid; | |
334 | } | |
335 | if(fDets[kHMPID]) { | |
336 | pTrk->GetHMPIDpid(pid[kHMPID]); | |
337 | status[kHMPID]=AliESDtrack::kHMPIDpid; | |
338 | } | |
339 | if(fDetRestr>=0){ | |
340 | if(fDetRestr == kITS) pTrk->GetITSpid(pid[kITS]); | |
341 | if(fDetRestr == kTPC) pTrk->GetITSpid(pid[kTPC]); | |
342 | if(fDetRestr == kTRD) pTrk->GetITSpid(pid[kTRD]); | |
343 | if(fDetRestr == kTOF) pTrk->GetITSpid(pid[kTOF]); | |
344 | if(fDetRestr == kHMPID) pTrk->GetITSpid(pid[kHMPID]); | |
345 | } | |
346 | ||
347 | status[kNdets]=pTrk->GetStatus(); | |
348 | pTrk->GetESDpid(pid[kNdets]); | |
349 | } | |
350 | //__________________________________ | |
351 | void AliCFTrackCutPid::SetPPriors(AliESDtrack *pTrk) | |
352 | { | |
353 | // | |
354 | //sets the mommentum dependent a priori concentrations | |
355 | // | |
356 | ||
357 | for(Int_t i=0; i< AliPID::kSPECIES; i++) { | |
358 | if(pTrk->P()>fPriorsFunc[i]->GetXmin() && pTrk->P() < fPriorsFunc[i]->GetXmax()) fPriors[i]=fPriorsFunc[i]->Eval(pTrk->P()); | |
13242232 | 359 | else {AliInfo("the track momentum is not in the function range. Priors are equal"); fPriors[i] = 0.2;} |
563113d0 | 360 | } |
361 | } | |
362 | //______________________________________ | |
363 | ULong_t AliCFTrackCutPid::StatusForAND(ULong_t status[kNdets+1]) const | |
364 | { | |
365 | // | |
366 | //In case of AND of more detectors the AND-detector status combination. | |
367 | //is calculated and also returned | |
368 | // | |
369 | ||
370 | ULong_t andstatus=0; | |
371 | for(Int_t i=0; i< kNdets; i++) { | |
372 | if(!fDetsInAnd[i]) continue; | |
373 | andstatus = andstatus | status[i]; | |
374 | AliDebug(1,Form("trk status %lu %i AND-status combination: %lu",status[i],i,andstatus)); | |
375 | } | |
376 | return andstatus; | |
377 | } | |
378 | //_______________________________________ | |
379 | Int_t AliCFTrackCutPid::GetID(ULong_t status[kNdets+1],Double_t pid[kNdets+1][AliPID::kSPECIES]) const | |
380 | { | |
381 | // Identifies the track if its probability is higher than the cut | |
382 | // value. The default value is fCut=0.2, therefore the most probable | |
383 | // one is identified by default. Here all the checks on how to identify | |
384 | // the track are performed (single detector or combined PID,..., detector | |
385 | // restriction probability | |
386 | // Returns: integer corresponding to the identified particle | |
387 | ||
388 | Int_t iPart=-1; | |
389 | ||
390 | if(!fgIsComb){ | |
391 | Bool_t isDet=kFALSE; | |
392 | for(Int_t i=0; i<kNdets; i++){ | |
393 | if(!fDets[i]) continue; | |
394 | isDet=kTRUE; | |
395 | AliDebug(1,Form("trk status %lu %i-det-pid status %lu -> combination: %lu",status[kNdets],i,status[i],status[kNdets]&status[i])); | |
396 | if(!(status[kNdets]&status[i])){ | |
397 | iPart=-10; | |
398 | AliDebug(1,Form("detector %i -> pid trk status not ok",i)); | |
399 | } | |
400 | else { | |
401 | AliDebug(1,Form("resp : %f %f %f %f %f",pid[i][0],pid[i][1],pid[i][2],pid[i][3],pid[i][4])); | |
402 | if(fIsQAOn) iPart = IdentifyQA(pid[i],i); | |
403 | else iPart = Identify(pid[i]); | |
404 | } | |
405 | } | |
406 | if(!isDet){ | |
407 | AliDebug(1,Form(" !!! No detector selected, the ESD-pid response is considered")); | |
408 | iPart = Identify(pid[kNdets]); | |
409 | } | |
410 | }else{ | |
411 | Double_t calcprob[5]; | |
f9658e11 | 412 | CombPID(status,pid,calcprob); |
563113d0 | 413 | iPart = Identify(calcprob); |
f9658e11 | 414 | |
563113d0 | 415 | } |
416 | ||
417 | ||
418 | AliDebug(1,Form("selected particle: %i",iPart)); | |
419 | ||
420 | if(iPart >=0 && fiPartRestr>=0) { | |
421 | AliPID restr(pid[fDetRestr]); | |
422 | restr.SetPriors(fPriors); | |
423 | AliDebug(1,Form("setted upper limit: %f det %i : probability %f ",fDetProbRestr,fDetRestr,restr.GetProbability((AliPID::EParticleType)fiPartRestr))); | |
424 | if(restr.GetProbability((AliPID::EParticleType)fiPartRestr) > fDetProbRestr) { | |
425 | iPart = kDetRestr; | |
426 | AliDebug(1,"\n\n the detector restrictions refused the ID \n\n"); | |
427 | } | |
428 | }//cross checks with one detector probability | |
429 | ||
430 | AliDebug(1,Form("after the check the selected particle is %i",iPart)); | |
431 | ||
432 | return iPart; | |
433 | } | |
bc6176fb | 434 | //_________________________________________________________________________________ |
435 | Int_t AliCFTrackCutPid::GetAODID(AliAODTrack *aodtrack) const | |
436 | { | |
437 | // | |
438 | // Identifies the AOD Track using the combined pid responses | |
439 | // | |
440 | ||
441 | Double_t combpid[AliPID::kSPECIES]; | |
f9658e11 | 442 | for(Int_t i=0; i< AliPID::kSPECIES; i++) { |
443 | combpid[i]= aodtrack->PID()[i]; | |
444 | if(!fhCombResp[i]) AliDebug(1,Form("\n no fhCombResp[%i], check if pidcut->Init() was called",i)); | |
445 | else fhCombResp[i]->Fill(combpid[i]); | |
446 | } | |
bc6176fb | 447 | return Identify(combpid); |
448 | } | |
563113d0 | 449 | //__________________________________ |
450 | Bool_t AliCFTrackCutPid::Check(const Double_t *p, Int_t iPsel, Double_t minDiff) const | |
451 | { | |
452 | // Checks if there are no equal values and if the valus | |
453 | // difference between the selected particle and the others i | |
454 | // is higher than a lower limit. | |
455 | // Returns: kTRUE= is acceptable | |
456 | ||
7f0c5c1e | 457 | AliDebug(2,Form("input particle: %i",iPsel)); |
563113d0 | 458 | Bool_t ck=kTRUE; |
459 | ||
460 | if(iPsel<0) ck=kFALSE; | |
461 | ||
462 | else { | |
463 | for(Int_t j=0; j< AliPID::kSPECIES; j++) { | |
464 | if(j!=iPsel && TMath::Abs(p[j]-p[iPsel])<minDiff) ck=kFALSE; | |
465 | } | |
466 | if(!ck) AliDebug(1,"the values are too close "); | |
467 | } | |
468 | return ck; | |
469 | } | |
470 | //___________________________________________ | |
471 | Int_t AliCFTrackCutPid::Identify(Double_t pid[AliPID::kSPECIES]) const | |
472 | { | |
473 | // | |
474 | // The identification is actually performed here with possible | |
475 | // checks on the det responses and/or probabilities | |
476 | // | |
bc6176fb | 477 | |
563113d0 | 478 | Int_t iPart = -1; |
bc6176fb | 479 | |
480 | AliDebug(2,Form("calc response bef: %f %f %f %f %f",pid[0],pid[1],pid[2],pid[3],pid[4])); | |
f9658e11 | 481 | AliDebug(2,Form("priors : %f %f %f %f %f",fPriors[0],fPriors[1],fPriors[2],fPriors[3],fPriors[4])); |
482 | ||
563113d0 | 483 | AliPID getpid(pid,kTRUE); |
bc6176fb | 484 | getpid.SetPriors(fPriors); |
563113d0 | 485 | |
bc6176fb | 486 | Double_t probability[AliPID::kSPECIES]={0.,0.,0.,0.,0.}; |
487 | for(Int_t iP=0; iP<AliPID::kSPECIES; iP++) { | |
488 | probability[iP] = getpid.GetProbability((AliPID::EParticleType)iP); | |
489 | AliDebug(2,Form("prob %i %f",iP, probability[iP])); | |
f9658e11 | 490 | if(fIsQAOn) fhCombProb[iP]->Fill(probability[iP]); |
bc6176fb | 491 | } |
563113d0 | 492 | |
39f34d57 | 493 | |
494 | if (fProbThreshold > 0.) { | |
495 | if (probability[fgParticleType] >= fProbThreshold) iPart=fgParticleType; | |
496 | } | |
497 | else { | |
498 | AliPID::EParticleType sel = getpid.GetMostProbable(); | |
499 | if(getpid.GetProbability(sel,fPriors)>fCut) iPart= (Int_t)sel; | |
500 | AliDebug(2,Form("probabilities : %f %f %f %f %f",probability[0],probability[1],probability[2],probability[3],probability[4])); | |
501 | } | |
7f0c5c1e | 502 | |
563113d0 | 503 | if(fCheckResponse && !Check(pid,iPart, fMinDiffResponse)) iPart=kCheckResp; |
504 | if(fCheckSelection && !Check(probability,iPart,fMinDiffProbability)) iPart=kCheckProb; | |
7f0c5c1e | 505 | |
563113d0 | 506 | return iPart; |
507 | } | |
508 | //___________________________________________ | |
509 | Int_t AliCFTrackCutPid::IdentifyQA(const Double_t pid[AliPID::kSPECIES], Int_t idets) const | |
510 | { | |
511 | // | |
512 | // The same as Identify, but with the QA histo filling | |
513 | // | |
514 | // | |
515 | ||
516 | Int_t iPart = -1; | |
517 | AliDebug(1,Form("resp : %f %f %f %f %f",pid[0],pid[1],pid[2],pid[3],pid[4])); | |
518 | ||
519 | AliPID getpid(pid,kTRUE); | |
f9658e11 | 520 | getpid.SetPriors(fPriors); |
521 | ||
563113d0 | 522 | AliPID::EParticleType sel = getpid.GetMostProbable(); |
523 | Double_t probability[AliPID::kSPECIES]; | |
524 | for(Int_t iP=0; iP<AliPID::kSPECIES; iP++) { | |
525 | probability[iP] = getpid.GetProbability((AliPID::EParticleType)iP); | |
526 | fhProb[idets][iP]->Fill(probability[iP]); | |
527 | } | |
528 | ||
ae8c412c | 529 | AliPID toresp(pid,kTRUE); |
530 | Double_t qapriors[10]={0.2,0.2,0.2,0.2,0.2,0,0,0,0,0}; | |
563113d0 | 531 | toresp.SetPriors(qapriors); |
532 | for(Int_t iPr=0; iPr<AliPID::kSPECIES; iPr++) fhResp[idets][iPr]->Fill(toresp.GetProbability((AliPID::EParticleType)iPr)); | |
533 | ||
534 | if(getpid.GetProbability(sel,fPriors)>fCut) iPart= (Int_t)sel; | |
535 | AliDebug(1,Form("resp : %f %f %f %f %f",pid[0],pid[1],pid[2],pid[3],pid[4])); | |
536 | AliDebug(1,Form("probab : %f %f %f %f %f",probability[0],probability[1],probability[2],probability[3],probability[4])); | |
537 | if(fCheckResponse && !Check(pid,iPart, fMinDiffResponse)) iPart=kCheckResp; | |
538 | if(fCheckSelection && !Check(probability,iPart,fMinDiffProbability)) iPart=kCheckProb; | |
539 | return iPart; | |
540 | } | |
541 | //___________________________________________ | |
9eeae5d5 | 542 | Bool_t AliCFTrackCutPid::IsSelected(TObject *track){ |
563113d0 | 543 | // |
544 | // method for the pid-cut selction | |
545 | // | |
bc6176fb | 546 | Bool_t sel = kFALSE; |
563113d0 | 547 | |
548 | if (!track) return kFALSE ; | |
549 | TString className(track->ClassName()); | |
bc6176fb | 550 | if (className.CompareTo("AliESDtrack") == 0) { |
f9658e11 | 551 | AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track); |
ae8c412c | 552 | if (!esdTrack) return kFALSE; |
563113d0 | 553 | ULong_t status[kNdets+1]={0,0,0,0,0,0}; |
554 | Double_t pid[kNdets+1][AliPID::kSPECIES]; | |
555 | TrackInfo(esdTrack,status,pid); | |
556 | if(fIsPpriors) SetPPriors(esdTrack); | |
bc6176fb | 557 | if(GetID(status,pid)==fgParticleType) sel = kTRUE; |
558 | } | |
559 | ||
560 | if (className.CompareTo("AliAODTrack") == 0) { | |
f9658e11 | 561 | AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track); |
ae8c412c | 562 | if (!aodtrack) return kFALSE ; |
bc6176fb | 563 | if(GetAODID(aodtrack) == fgParticleType) sel = kTRUE; |
564 | } | |
565 | ||
566 | return sel; | |
567 | ||
563113d0 | 568 | } |
569 | //__________________________________ | |
570 | void AliCFTrackCutPid::CombPID(ULong_t status[kNdets+1],Double_t pid[kNdets+1][AliPID::kSPECIES],Double_t *combpid) const | |
571 | { | |
572 | // | |
573 | // Calculates the combined PID according to the chosen detectors. | |
574 | // and provides the array of probabilities | |
575 | // | |
576 | ||
577 | Bool_t isdet=kFALSE; | |
563113d0 | 578 | Double_t prod[AliPID::kSPECIES]={1.,1.,1.,1.,1.}; |
7f0c5c1e | 579 | |
563113d0 | 580 | ULong_t andstatus =0; |
581 | if(fIsDetAND) { | |
582 | andstatus = StatusForAND(status); | |
583 | AliDebug(1,Form("AND combination %lu",andstatus)); | |
584 | } | |
f9658e11 | 585 | |
586 | //Products of single detector responses | |
563113d0 | 587 | for(Int_t j=0; j<AliPID::kSPECIES; j++){ |
588 | for(Int_t i=0; i< kNdets; i++){ | |
589 | if(!fDets[i]) continue; | |
590 | if(status[kNdets]&status[i]) { | |
591 | if(fIsDetAND) { | |
592 | ULong_t checkstatus = status[kNdets]&andstatus; | |
593 | if(checkstatus != andstatus) continue; | |
594 | else { | |
595 | prod[j]*=pid[i][j]; | |
596 | isdet = kTRUE; | |
597 | AliDebug(1,Form("-----> trk status %lu and status %lu -> trk-ANDdetector status combination %lu",status[kNdets],andstatus,status[kNdets]&andstatus)); | |
7f0c5c1e | 598 | AliDebug(1,Form("In det %i -> particle %i response is %f",i,j,pid[i][j])); |
563113d0 | 599 | } |
f9658e11 | 600 | } else { |
563113d0 | 601 | prod[j]*=pid[i][j]; |
602 | isdet=kTRUE; | |
7f0c5c1e | 603 | AliDebug(2,Form("In det %i -> particle %i response is %f",i,j,pid[i][j])); |
f9658e11 | 604 | |
605 | if(fIsQAOn){ | |
606 | if(!fhResp[i][j]) {AliDebug(1,Form("no pointer to the histo fhResp%i%i, check if pidcut->Init() was called",i,j));} | |
607 | else fhResp[i][j]->Fill(pid[i][j]); | |
608 | ||
609 | if(!fhProb[i][j]) {AliDebug(1,Form("no pointer to the histo fhProb%i%i, check if pidcut->Init() was called",i,j));} | |
610 | else { | |
611 | AliPID detprob(pid[i],kTRUE); | |
612 | detprob.SetPriors(fPriors); | |
613 | fhProb[i][j]->Fill(detprob.GetProbability((AliPID::EParticleType)j)); | |
614 | } | |
615 | } | |
563113d0 | 616 | } |
617 | }//combined mode | |
618 | }//loop on dets | |
619 | }//loop on species | |
7f0c5c1e | 620 | |
f9658e11 | 621 | |
622 | //no detectors found, then go to ESD pid... | |
563113d0 | 623 | if(!isdet) { |
7f0c5c1e | 624 | AliWarning("\n !! No detector found for the combined pid --> responses are from the ESDpid !! \n"); |
625 | Double_t sumesdpid=0; | |
f9658e11 | 626 | for(Int_t nn=0; nn<AliPID::kSPECIES; nn++) sumesdpid+=pid[kNdets][nn]; |
627 | if(sumesdpid<=0) { | |
628 | AliDebug(1,"priors or ESDpid are inconsistent, please check them"); | |
629 | return; | |
630 | } else { | |
631 | for(Int_t k=0; k<AliPID::kSPECIES; k++){ | |
632 | combpid[k] = pid[kNdets][k]/sumesdpid; | |
633 | if(fIsQAOn) { | |
634 | if(!fhCombResp[k]) AliDebug(1,Form("\n no fhCombResp[%i], check if pidcut->Init() was called",k)); | |
635 | else fhCombResp[k]->Fill(combpid[k]); | |
636 | } | |
637 | }//loop on species | |
638 | } | |
639 | return; | |
640 | } | |
563113d0 | 641 | |
f9658e11 | 642 | Double_t add = 0; for(Int_t isumm=0; isumm<5; isumm++) add+=prod[isumm]; |
643 | if(add>0) { | |
644 | for(Int_t ip =0; ip < AliPID::kSPECIES; ip++) { | |
645 | combpid[ip] = prod[ip]/add; | |
646 | if(fIsQAOn) { | |
647 | if(!fhCombResp[ip]) AliDebug(1,Form("\n no fhCombResp[%i], check if pidcut->Init() was called",ip)); | |
648 | else fhCombResp[ip]->Fill(combpid[ip]); | |
649 | ||
650 | } | |
651 | } | |
652 | AliDebug(1,Form("calculated comb response: %f %f %f %f %f",combpid[0],combpid[1],combpid[2],combpid[3],combpid[4])); | |
653 | } else { | |
654 | AliDebug(1,"single detector responses are inconsistent, please check them...."); | |
655 | return; | |
656 | } | |
657 | AliDebug(1,Form("the ESDpid response: %f %f %f %f %f",pid[kNdets][0],pid[kNdets][1],pid[kNdets][2],pid[kNdets][3],pid[kNdets][4])); | |
658 | ||
563113d0 | 659 | } |
660 | //__________________________________________ | |
661 | // | |
662 | //QA part | |
663 | //_________________________________________ | |
664 | void AliCFTrackCutPid::InitialiseHisto() | |
665 | { | |
666 | // | |
667 | //QA histo initialization | |
668 | // | |
669 | for(Int_t iP=0; iP<AliPID::kSPECIES; iP++){ | |
670 | fhCombResp[iP]=0x0; | |
671 | fhCombProb[iP]=0x0; | |
672 | for(Int_t iDet =0; iDet<kNdets; iDet++){ | |
673 | fhResp[iDet][iP]=0x0; | |
674 | fhProb[iDet][iP]=0x0; | |
675 | } | |
676 | } | |
677 | } | |
678 | //______________________________________________ | |
563113d0 | 679 | void AliCFTrackCutPid::DefineHistograms() |
680 | { | |
681 | // | |
682 | //QA histo booking | |
683 | // | |
7f0c5c1e | 684 | |
bc6176fb | 685 | if(fgIsAOD){ |
1f2d22de | 686 | const char *partic[AliPID::kSPECIES]={"electron","muon","pion","kaon","proton"}; |
bc6176fb | 687 | |
688 | for(Int_t iPart =0; iPart < AliPID::kSPECIES; iPart++) | |
563113d0 | 689 | { |
ce2a1393 | 690 | fhCombResp[iPart] = new TH1F(Form("%s_rCombPart%i",GetName(),iPart),Form(" %s combined response (AODTrack) ",partic[iPart]),fNbins,fXmin,fXmax); |
f9658e11 | 691 | fhCombResp[iPart]->SetXTitle(Form(" %s combined response ",partic[iPart])); |
692 | fhCombResp[iPart]->SetYTitle("entries"); | |
ce2a1393 | 693 | AliDebug(1,Form( "%s is booked!!",fhCombResp[iPart]->GetName())); |
694 | fhCombProb[iPart] = new TH1F(Form("%s_pCombPart%i",GetName(),iPart),Form("%s combined probability (AODTrack) ",partic[iPart]),fNbins,fXmin,fXmax); | |
695 | fhCombProb[iPart]->SetXTitle(Form(" %s combined probability ",partic[iPart])); | |
f9658e11 | 696 | fhCombProb[iPart]->SetYTitle("entries"); |
ce2a1393 | 697 | AliDebug(1,Form( "%s is booked!!",fhCombProb[iPart]->GetName())); |
bc6176fb | 698 | } |
699 | } | |
700 | ||
701 | ||
702 | else { | |
1f2d22de | 703 | const char *detect[kNdets]={"ITS","TPC","TRD","TOF","HMPID"}; |
704 | const char *partic[AliPID::kSPECIES]={"electron","muon","pion","kaon","proton"}; | |
bc6176fb | 705 | |
706 | for(Int_t iDet =0; iDet< kNdets; iDet++) | |
707 | { | |
708 | if(!fDets[iDet]) continue; | |
709 | for(Int_t iP =0; iP < AliPID::kSPECIES; iP++){ | |
ce2a1393 | 710 | fhResp[iDet][iP] = new TH1F(Form("%s_rDet%iPart%i",GetName(),iDet,iP),Form("%s response for %s ",detect[iDet],partic[iP]),fNbins,fXmin,fXmax); |
f9658e11 | 711 | fhResp[iDet][iP]->SetXTitle(Form(" %s response ",partic[iP])); |
712 | fhResp[iDet][iP]->SetYTitle("entries"); | |
ce2a1393 | 713 | fhProb[iDet][iP] = new TH1F(Form("%s_pDet%iPart%i",GetName(),iDet,iP),Form("%s calculated probability for %s",detect[iDet],partic[iP]),fNbins,fXmin,fXmax); |
714 | fhProb[iDet][iP]->SetXTitle(Form(" %s probability ",partic[iP])); | |
f9658e11 | 715 | fhProb[iDet][iP]->SetYTitle("entries"); |
bc6176fb | 716 | } |
717 | } | |
718 | ||
719 | ||
720 | if(fgIsComb) | |
721 | { | |
722 | for(Int_t iPart =0; iPart < AliPID::kSPECIES; iPart++) | |
723 | { | |
ce2a1393 | 724 | fhCombResp[iPart] = new TH1F(Form("%s_rCombPart%i",GetName(),iPart),Form(" %s combined response ",partic[iPart]),fNbins,fXmin,fXmax); |
f9658e11 | 725 | fhCombResp[iPart]->SetXTitle(Form(" %s response ",partic[iPart])); |
726 | fhCombResp[iPart]->SetYTitle("entries"); | |
ce2a1393 | 727 | AliDebug(1,Form( "%s is booked!!",fhCombResp[iPart]->GetName())); |
728 | fhCombProb[iPart] = new TH1F(Form("%s_pCombPart%i",GetName(),iPart),Form("%s combined probability ",partic[iPart]),fNbins,fXmin,fXmax); | |
f9658e11 | 729 | fhCombProb[iPart]->SetXTitle(Form(" %s response ",partic[iPart])); |
730 | fhCombProb[iPart]->SetYTitle("entries"); | |
ce2a1393 | 731 | AliDebug(1,Form( "%s is booked!!",fhCombProb[iPart]->GetName())); |
bc6176fb | 732 | } |
733 | } | |
734 | } | |
735 | ||
563113d0 | 736 | } |
737 | //___________________________________________________ | |
738 | ||
107a3100 | 739 | void AliCFTrackCutPid::AddQAHistograms(TList *qalist) |
563113d0 | 740 | { |
741 | // | |
742 | // adds QA histograms in a TList | |
743 | // | |
744 | if(!fIsQAOn) return; | |
107a3100 | 745 | DefineHistograms(); |
563113d0 | 746 | |
f9658e11 | 747 | if(fgIsComb || fgIsAOD){ |
563113d0 | 748 | for(Int_t iPart =0; iPart<AliPID::kSPECIES; iPart++){ |
749 | qalist->Add(fhCombResp[iPart]); | |
750 | qalist->Add(fhCombProb[iPart]); | |
751 | } | |
752 | } | |
753 | ||
754 | for(Int_t iDet=0; iDet<kNdets; iDet++){ | |
7f0c5c1e | 755 | if(!fDets[iDet]) continue; |
563113d0 | 756 | for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){ |
7f0c5c1e | 757 | qalist->Add(fhResp[iDet][iP]); |
758 | qalist->Add(fhProb[iDet][iP]); | |
563113d0 | 759 | } |
760 | } | |
7f0c5c1e | 761 | |
563113d0 | 762 | } |
f9658e11 | 763 |