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