1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
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);
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)
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:
51 //pCut->SetMinDiffProb(0.005,kTRUE) //for probabilities
53 //pCut->SetMinDiffResp(0.005,kTRUE) //for det responses
55 //Annalisa.Mastroserio@ba.infn.it
59 #include "AliCFTrackCutPid.h"
64 ClassImp(AliCFTrackCutPid)
66 //__________________________________________________________________________________
68 //__________________________________________________________________________________
69 AliCFTrackCutPid::AliCFTrackCutPid() :
72 fMinDiffResponse(0.0001),
73 fMinDiffProbability(0.001),
76 fCheckResponse(kFALSE),
77 fCheckSelection(kTRUE),
90 for(Int_t j=0; j< AliPID::kSPECIES; j++) {
95 for(Int_t jDet=0; jDet< kNdets; jDet++) {
97 fDetsInAnd[jDet]=kFALSE;
102 //______________________________
103 AliCFTrackCutPid::AliCFTrackCutPid(const Char_t* name, const Char_t* title) :
104 AliCFCutBase(name,title),
106 fMinDiffResponse(0.0001),
107 fMinDiffProbability(0.001),
110 fCheckResponse(kFALSE),
111 fCheckSelection(kTRUE),
124 for(Int_t j=0; j< AliPID::kSPECIES; j++) {
129 for(Int_t jDet=0; jDet< kNdets; jDet++) {
131 fDetsInAnd[jDet]=kFALSE;
136 //______________________________
137 AliCFTrackCutPid::AliCFTrackCutPid(const AliCFTrackCutPid& c) :
140 fMinDiffResponse(c.fMinDiffResponse),
141 fMinDiffProbability(c.fMinDiffProbability),
142 fgParticleType(c.fgParticleType),
143 fgIsComb(c.fgIsComb),
144 fCheckResponse(c.fCheckResponse),
145 fCheckSelection(c.fCheckSelection),
146 fIsPpriors(c.fIsPpriors),
147 fIsDetAND(c.fIsDetAND),
151 fDetRestr(c.fDetRestr),
152 fiPartRestr(c.fiPartRestr),
153 fDetProbRestr(c.fDetProbRestr)
158 for(Int_t i=0; i< kNdets ; i++ ) {
160 fDetsInAnd[i]=c.fDetsInAnd[i];
161 for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){
162 fhResp[i][iP]=c.fhResp[i][iP];
163 fhProb[i][iP]=c.fhProb[i][iP];
166 for(Int_t j=0; j< AliPID::kSPECIES; j++){
167 fPriors[j]=c.fPriors[j];
168 fPriorsFunc[j]=c.fPriorsFunc[j];
169 fhCombResp[j]=c.fhCombResp[j];
170 fhCombProb[j]=c.fhCombProb[j];
173 //______________________________
174 AliCFTrackCutPid& AliCFTrackCutPid::operator=(const AliCFTrackCutPid& c)
177 // Assignment operator
180 AliCFCutBase::operator=(c) ;
182 this->fMinDiffResponse=c.fMinDiffResponse;
183 this->fMinDiffProbability=c.fMinDiffProbability;
184 this->fgParticleType=c.fgParticleType;
185 this->fgIsComb=c.fgIsComb;
186 this->fCheckResponse=c.fCheckResponse;
187 this->fCheckSelection=c.fCheckSelection;
188 this->fIsPpriors=c.fIsPpriors;
189 this->fIsDetAND=c.fIsDetAND;
192 this->fNbins=c.fNbins;
193 this->fDetRestr=c.fDetRestr;
194 this->fiPartRestr=c.fiPartRestr;
195 this->fDetProbRestr=c.fDetProbRestr;
197 for(Int_t i=0; i< kNdets ; i++ ) {
198 this->fDets[i]=c.fDets[i];
199 for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){
200 this->fhResp[i][iP]=c.fhResp[i][iP];
201 this->fhProb[i][iP]=c.fhProb[i][iP];
205 for(Int_t j=0; j< AliPID::kSPECIES; j++){
206 this->fPriors[j]=c.fPriors[j];
207 this->fhCombResp[j]=c.fhCombResp[j];
208 this->fhCombProb[j]=c.fhCombProb[j];
209 this-> fPriorsFunc[j]=c.fPriorsFunc[j];
215 //__________________________________________________________________________________
216 AliCFTrackCutPid::~AliCFTrackCutPid() {
220 for(Int_t i=0; i< kNdets ; i++ ) {
221 for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){
222 if(fhResp[i][iP])delete fhResp[i][iP];
223 if(fhProb[i][iP])delete fhProb[i][iP];
227 for(Int_t j=0; j< AliPID::kSPECIES; j++){
228 if(fhCombResp[j])delete fhCombResp[j];
229 if(fhCombProb[j])delete fhCombProb[j];
234 //__________________________________
235 void AliCFTrackCutPid::SetDetectors(TString dets)
238 // The string of detectors is translated into
239 // the respective booelan data members
241 if(dets.Contains("ITS")) {fDets[kITS]=kTRUE;}
242 if(dets.Contains("TPC")) {fDets[kTPC]=kTRUE;}
243 if(dets.Contains("TRD")) {fDets[kTRD]=kTRUE;}
244 if(dets.Contains("TOF")) {fDets[kTOF]=kTRUE;}
245 if(dets.Contains("HMPID")) {fDets[kHMPID]=kTRUE;}
247 //__________________________________
248 void AliCFTrackCutPid::SetPriors(Double_t r[AliPID::kSPECIES])
251 // Sets the a priori concentrations
253 for(Int_t i=0; i<AliPID::kSPECIES; i++) fPriors[i]=r[i];
255 //__________________________________
256 void AliCFTrackCutPid::SetPriorFunctions(TF1 *func[AliPID::kSPECIES])
259 // Sets the momentu dependent a priori concentrations
262 for(Int_t i=0; i<AliPID::kSPECIES; i++) fPriorsFunc[i]=func[i];
265 //____________________________________________
266 void AliCFTrackCutPid::SetANDstatus(TString dets)
269 //Sets the detectors to be in AND for the combined PID
271 if(dets.Contains("ITS") && fDets[kITS]) {fDetsInAnd[kITS]=kTRUE;}
272 if(dets.Contains("TPC") && fDets[kTPC]) {fDetsInAnd[kTPC]=kTRUE;}
273 if(dets.Contains("TRD") && fDets[kTRD]) {fDetsInAnd[kTRD]=kTRUE;}
274 if(dets.Contains("TOF") && fDets[kTOF]) {fDetsInAnd[kTOF]=kTRUE;}
275 if(dets.Contains("HMPID") && fDets[kHMPID]) {fDetsInAnd[kHMPID]=kTRUE;}
280 void AliCFTrackCutPid::SetDetectorProbabilityRestriction(TString det, Int_t iPart, Double_t upperprob)
283 // Sets the detector, the particle and the probability
287 if(det.Contains("ITS")) fDetRestr = kITS;
288 if(det.Contains("TPC")) fDetRestr = kTPC;
289 if(det.Contains("TRD")) fDetRestr = kTRD;
290 if(det.Contains("TOF")) fDetRestr = kTOF;
291 if(det.Contains("HMPID")) fDetRestr = kHMPID;
293 fDetProbRestr = upperprob;
295 //__________________________________
296 void AliCFTrackCutPid::TrackInfo(const AliESDtrack *pTrk, ULong_t status[kNdets+1],Double_t pid[kNdets+1][AliPID::kSPECIES]) const
299 // Loads the responses of the XXX chosen detectors for the pTrk
300 // and the corresponding trk status. The final trk status is also loaded.
303 pTrk->GetITSpid(pid[kITS]);
304 status[kITS]=AliESDtrack::kITSpid;
307 pTrk->GetTPCpid(pid[kTPC]);
308 status[kTPC]=AliESDtrack::kTPCpid;
311 pTrk->GetTRDpid(pid[kTRD]);
312 status[kTRD]=AliESDtrack::kTRDpid;
315 pTrk->GetTOFpid(pid[kTOF]);
316 status[kTOF]=AliESDtrack::kTOFpid;
319 pTrk->GetHMPIDpid(pid[kHMPID]);
320 status[kHMPID]=AliESDtrack::kHMPIDpid;
323 if(fDetRestr == kITS) pTrk->GetITSpid(pid[kITS]);
324 if(fDetRestr == kTPC) pTrk->GetITSpid(pid[kTPC]);
325 if(fDetRestr == kTRD) pTrk->GetITSpid(pid[kTRD]);
326 if(fDetRestr == kTOF) pTrk->GetITSpid(pid[kTOF]);
327 if(fDetRestr == kHMPID) pTrk->GetITSpid(pid[kHMPID]);
330 status[kNdets]=pTrk->GetStatus();
331 pTrk->GetESDpid(pid[kNdets]);
333 //__________________________________
334 void AliCFTrackCutPid::SetPPriors(AliESDtrack *pTrk)
337 //sets the mommentum dependent a priori concentrations
340 for(Int_t i=0; i< AliPID::kSPECIES; i++) {
341 if(pTrk->P()>fPriorsFunc[i]->GetXmin() && pTrk->P() < fPriorsFunc[i]->GetXmax()) fPriors[i]=fPriorsFunc[i]->Eval(pTrk->P());
342 else {AliInfo("the track momentum is not in the function range. Priors are equal") fPriors[i] = 0.2;}
345 //______________________________________
346 ULong_t AliCFTrackCutPid::StatusForAND(ULong_t status[kNdets+1]) const
349 //In case of AND of more detectors the AND-detector status combination.
350 //is calculated and also returned
354 for(Int_t i=0; i< kNdets; i++) {
355 if(!fDetsInAnd[i]) continue;
356 andstatus = andstatus | status[i];
357 AliDebug(1,Form("trk status %lu %i AND-status combination: %lu",status[i],i,andstatus));
361 //_______________________________________
362 Int_t AliCFTrackCutPid::GetID(ULong_t status[kNdets+1],Double_t pid[kNdets+1][AliPID::kSPECIES]) const
364 // Identifies the track if its probability is higher than the cut
365 // value. The default value is fCut=0.2, therefore the most probable
366 // one is identified by default. Here all the checks on how to identify
367 // the track are performed (single detector or combined PID,..., detector
368 // restriction probability
369 // Returns: integer corresponding to the identified particle
375 for(Int_t i=0; i<kNdets; i++){
376 if(!fDets[i]) continue;
378 AliDebug(1,Form("trk status %lu %i-det-pid status %lu -> combination: %lu",status[kNdets],i,status[i],status[kNdets]&status[i]));
379 if(!(status[kNdets]&status[i])){
381 AliDebug(1,Form("detector %i -> pid trk status not ok",i));
384 AliDebug(1,Form("resp : %f %f %f %f %f",pid[i][0],pid[i][1],pid[i][2],pid[i][3],pid[i][4]));
385 if(fIsQAOn) iPart = IdentifyQA(pid[i],i);
386 else iPart = Identify(pid[i]);
390 AliDebug(1,Form(" !!! No detector selected, the ESD-pid response is considered"));
391 iPart = Identify(pid[kNdets]);
394 Double_t calcprob[5];
395 CombPID(status,pid,calcprob);
396 iPart = Identify(calcprob);
400 AliDebug(1,Form("selected particle: %i",iPart));
402 if(iPart >=0 && fiPartRestr>=0) {
403 AliPID restr(pid[fDetRestr]);
404 restr.SetPriors(fPriors);
405 AliDebug(1,Form("setted upper limit: %f det %i : probability %f ",fDetProbRestr,fDetRestr,restr.GetProbability((AliPID::EParticleType)fiPartRestr)));
406 if(restr.GetProbability((AliPID::EParticleType)fiPartRestr) > fDetProbRestr) {
408 AliDebug(1,"\n\n the detector restrictions refused the ID \n\n");
410 }//cross checks with one detector probability
412 AliDebug(1,Form("after the check the selected particle is %i",iPart));
416 //__________________________________
417 Bool_t AliCFTrackCutPid::Check(const Double_t *p, Int_t iPsel, Double_t minDiff) const
419 // Checks if there are no equal values and if the valus
420 // difference between the selected particle and the others i
421 // is higher than a lower limit.
422 // Returns: kTRUE= is acceptable
424 AliDebug(1,Form("input particle: %i",iPsel));
427 if(iPsel<0) ck=kFALSE;
430 for(Int_t j=0; j< AliPID::kSPECIES; j++) {
431 if(j!=iPsel && TMath::Abs(p[j]-p[iPsel])<minDiff) ck=kFALSE;
433 if(!ck) AliDebug(1,"the values are too close ");
437 //___________________________________________
438 Int_t AliCFTrackCutPid::Identify(Double_t pid[AliPID::kSPECIES]) const
441 // The identification is actually performed here with possible
442 // checks on the det responses and/or probabilities
446 AliPID getpid(pid,kTRUE);
449 Double_t priors[5]={0.2,0.2,0.2,0.2,0.2};
450 getpid.SetPriors(priors);
452 else getpid.SetPriors(fPriors);
455 AliPID::EParticleType sel = getpid.GetMostProbable();
456 Double_t probability[AliPID::kSPECIES];
457 for(Int_t iP=0; iP<AliPID::kSPECIES; iP++) probability[iP] = getpid.GetProbability((AliPID::EParticleType)iP);
460 if(getpid.GetProbability(sel,fPriors)>fCut) iPart= (Int_t)sel;
461 AliDebug(1,Form("resp : %f %f %f %f %f",pid[0],pid[1],pid[2],pid[3],pid[4]));
462 AliDebug(1,Form("probab : %f %f %f %f %f",probability[0],probability[1],probability[2],probability[3],probability[4]));
463 if(fCheckResponse && !Check(pid,iPart, fMinDiffResponse)) iPart=kCheckResp;
464 if(fCheckSelection && !Check(probability,iPart,fMinDiffProbability)) iPart=kCheckProb;
467 //___________________________________________
468 Int_t AliCFTrackCutPid::IdentifyQA(const Double_t pid[AliPID::kSPECIES], Int_t idets) const
471 // The same as Identify, but with the QA histo filling
476 AliDebug(1,Form("resp : %f %f %f %f %f",pid[0],pid[1],pid[2],pid[3],pid[4]));
478 AliPID getpid(pid,kTRUE);
481 Double_t priors[5]={0.2,0.2,0.2,0.2,0.2};
482 getpid.SetPriors(priors);
484 else getpid.SetPriors(fPriors);
486 AliPID::EParticleType sel = getpid.GetMostProbable();
487 Double_t probability[AliPID::kSPECIES];
488 for(Int_t iP=0; iP<AliPID::kSPECIES; iP++) {
489 probability[iP] = getpid.GetProbability((AliPID::EParticleType)iP);
490 fhProb[idets][iP]->Fill(probability[iP]);
493 AliPID toresp(pid,kTRUE); Double_t qapriors[5]={0.2,0.2,0.2,0.2,0.2};
494 toresp.SetPriors(qapriors);
495 for(Int_t iPr=0; iPr<AliPID::kSPECIES; iPr++) fhResp[idets][iPr]->Fill(toresp.GetProbability((AliPID::EParticleType)iPr));
497 if(getpid.GetProbability(sel,fPriors)>fCut) iPart= (Int_t)sel;
498 AliDebug(1,Form("resp : %f %f %f %f %f",pid[0],pid[1],pid[2],pid[3],pid[4]));
499 AliDebug(1,Form("probab : %f %f %f %f %f",probability[0],probability[1],probability[2],probability[3],probability[4]));
500 if(fCheckResponse && !Check(pid,iPart, fMinDiffResponse)) iPart=kCheckResp;
501 if(fCheckSelection && !Check(probability,iPart,fMinDiffProbability)) iPart=kCheckProb;
504 //___________________________________________
505 Bool_t AliCFTrackCutPid::IsSelected(TObject *track){
507 // method for the pid-cut selction
510 if (!track) return kFALSE ;
511 TString className(track->ClassName());
512 if (className.CompareTo("AliESDtrack") != 0) {
513 AliError("obj must point to a AliESDtrack ");
517 AliESDtrack *esdTrack = (AliESDtrack *)track;
518 ULong_t status[kNdets+1]={0,0,0,0,0,0};
519 Double_t pid[kNdets+1][AliPID::kSPECIES];
520 TrackInfo(esdTrack,status,pid);
521 if(fIsPpriors) SetPPriors(esdTrack);
522 if(GetID(status,pid)==fgParticleType) return kTRUE;
525 //__________________________________
526 void AliCFTrackCutPid::CombPID(ULong_t status[kNdets+1],Double_t pid[kNdets+1][AliPID::kSPECIES],Double_t *combpid) const
529 // Calculates the combined PID according to the chosen detectors.
530 // and provides the array of probabilities
535 Double_t prod[AliPID::kSPECIES]={1.,1.,1.,1.,1.};
536 Double_t comb[AliPID::kSPECIES]={0.,0.,0.,0.,0.};
537 Double_t priors[AliPID::kSPECIES]={0.2,0.2,0.2,0.2,0.2};
538 ULong_t andstatus =0;
540 andstatus = StatusForAND(status);
541 AliDebug(1,Form("AND combination %lu",andstatus));
543 for(Int_t j=0; j<AliPID::kSPECIES; j++){
544 for(Int_t i=0; i< kNdets; i++){
545 if(!fDets[i]) continue;
546 if(status[kNdets]&status[i]) {
548 ULong_t checkstatus = status[kNdets]&andstatus;
549 if(checkstatus != andstatus) continue;
553 AliDebug(1,Form("-----> trk status %lu and status %lu -> trk-ANDdetector status combination %lu",status[kNdets],andstatus,status[kNdets]&andstatus));
554 AliDebug(1,Form("In det loop %i -> particle %i response is %f",i,j,pid[i][j]));
560 AliDebug(1,Form("In det loop %i -> particle %i response is %f",i,j,pid[i][j]));
567 for(Int_t iqa =0; iqa < kNdets; iqa++){
568 if(!fDets[iqa]) continue;
569 AliPID normresp(pid[iqa]);
570 normresp.SetPriors(priors);
571 for(Int_t ip =0; ip< AliPID::kSPECIES; ip++){
572 if(!fhResp[iqa][ip]) {AliDebug(1,Form("no pointer to the histo fhResp%i%i, check if pidcut->Init() was called",iqa,ip));}
573 else fhResp[iqa][ip]->Fill(normresp.GetProbability((AliPID::EParticleType)ip));
579 AliDebug(1,"No proper status for the combined pid -> probabilities are set to zero");
580 for(Int_t nn=0; nn<AliPID::kSPECIES; nn++) {combpid[nn]=0;}
584 for(Int_t k=0; k<AliPID::kSPECIES; k++){
586 if(!fhCombResp[k]) AliDebug(1,Form("no fhCombResp[%i], check if pidcut->Init() was called",k));
587 else fhCombResp[k]->Fill(prod[k]);
589 AliDebug(1,Form("species %i priors %f and prod %f",k,fPriors[k],prod[k]));
590 comb[k]=fPriors[k]*prod[k];
591 AliDebug(1,Form("comb %i %f",k,comb[k]));
596 AliDebug(1,"Check the detector responses or the priors, their combined products are zero");
599 for(Int_t n=0; n<AliPID::kSPECIES; n++) {
600 combpid[n]=fPriors[n]*prod[n]/sum;
602 if(!fhCombProb[n]) Printf("no fhCombResp[%i], check if pidcut->Init() was called",n);
603 fhCombProb[n]->Fill(combpid[n]);
608 //__________________________________________
611 //_________________________________________
612 void AliCFTrackCutPid::InitialiseHisto()
615 //QA histo initialization
617 for(Int_t iP=0; iP<AliPID::kSPECIES; iP++){
620 for(Int_t iDet =0; iDet<kNdets; iDet++){
621 fhResp[iDet][iP]=0x0;
622 fhProb[iDet][iP]=0x0;
626 //______________________________________________
627 void AliCFTrackCutPid::Init()
630 // initialises QA histograms
633 if(fIsQAOn) DefineHistograms();
635 //_________________________________________________
636 void AliCFTrackCutPid::DefineHistograms()
641 char *detect[5]={"ITS","TPC","TRD","TOF","HMPID"};
642 char *partic[5]={"electron","muon","pion","kaon","proton"};
644 for(Int_t iDet =0; iDet< kNdets; iDet++)
646 if(!fDets[iDet]) continue;
647 for(Int_t iP =0; iP < AliPID::kSPECIES; iP++){
648 fhResp[iDet][iP] = new TH1F(Form("rDet%iPart%i",iDet,iP),Form("%s %s response ",detect[iDet],partic[iP]),fNbins,fXmin,fXmax);
649 fhProb[iDet][iP] = new TH1F(Form("pDet%iPart%i",iDet,iP),Form("%s %s probability ",detect[iDet],partic[iP]),fNbins,fXmin,fXmax);
654 for(Int_t iPart =0; iPart < AliPID::kSPECIES; iPart++)
656 fhCombResp[iPart] = new TH1F(Form("rCombPart%i",iPart),Form(" %s combined response ",partic[iPart]),fNbins,fXmin,fXmax);
657 fhCombProb[iPart] = new TH1F(Form("pCombPart%i",iPart),Form("%s combined probability ",partic[iPart]),fNbins,fXmin,fXmax);
661 //___________________________________________________
663 void AliCFTrackCutPid::AddQAHistograms(TList *qalist) const
666 // adds QA histograms in a TList
671 for(Int_t iPart =0; iPart<AliPID::kSPECIES; iPart++){
672 qalist->Add(fhCombResp[iPart]);
673 qalist->Add(fhCombProb[iPart]);
677 for(Int_t iDet=0; iDet<kNdets; iDet++){
678 for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){
679 if(!fgIsComb)qalist->Add(fhResp[iDet][iP]);
680 if(!fgIsComb)qalist->Add(fhProb[iDet][iP]);