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