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