]> git.uio.no Git - u/mrichter/AliRoot.git/blame - CORRFW/AliCFTrackCutPid.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
syst. check that reads the MC information for the case of TPC-only
[u/mrichter/AliRoot.git] / CORRFW / AliCFTrackCutPid.cxx
CommitLineData
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
64ClassImp(AliCFTrackCutPid)
65
66//__________________________________________________________________________________
67// CUT ON TRACK PID
68//__________________________________________________________________________________
69AliCFTrackCutPid::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//______________________________
104AliCFTrackCutPid::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//______________________________
139AliCFTrackCutPid::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//______________________________
177AliCFTrackCutPid& 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//__________________________________________________________________________________
220AliCFTrackCutPid::~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//__________________________________
239void 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//__________________________________
254void 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//__________________________________
262void 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//____________________________________________
272void 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//
286void 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//__________________________________
302void 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//__________________________________
340void 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//______________________________________
352ULong_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//_______________________________________
368Int_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//_________________________________________________________________________________
423Int_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//__________________________________
434Bool_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//___________________________________________
455Int_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//___________________________________________
495Int_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//___________________________________________
532Bool_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//__________________________________
558void 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//_________________________________________
665void 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 680void 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 729void 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}