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