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