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