CORRFW/AliCFTrackCutPid.cxx

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