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   fProbThreshold(0.)
  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 //______________________________
 106 AliCFTrackCutPid::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),
 113   fgIsAOD(kFALSE),
 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),
 123   fDetProbRestr(1),
 124   fProbThreshold(0.)
 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 //______________________________
 142 AliCFTrackCutPid::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),
 149   fgIsAOD(c.fgIsAOD),
 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),
 159   fDetProbRestr(c.fDetProbRestr),
 160   fProbThreshold(c.fProbThreshold)
 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 //______________________________
 181 AliCFTrackCutPid& 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;
 190     this->fMinDiffProbability=c.fMinDiffProbability;
 191     this->fgParticleType=c.fgParticleType;
 192     this->fgIsComb=c.fgIsComb;
 193     this->fgIsAOD=c.fgIsAOD;
 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;
 204     this->fProbThreshold=c.fProbThreshold;
 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 //__________________________________________________________________________________
 225 AliCFTrackCutPid::~AliCFTrackCutPid() {
 226   //
 227   //dtor
 228   //
 229
 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   }
 242 }
 243 //__________________________________
 244 void 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;}
 255
 256   if(dets.Contains("ALL")) for(Int_t i=0; i< kNdets ; i++) fDets[i]=kTRUE;
 257 }
 258 //__________________________________
 259 void 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 //__________________________________
 267 void 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 //____________________________________________
 277 void 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 //
 291 void 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 //__________________________________
 307 void 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 //__________________________________
 345 void 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 //______________________________________
 357 ULong_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 //_______________________________________
 373 Int_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];
 406     CombPID(status,pid,calcprob);
 407     iPart = Identify(calcprob);
 408
 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 }
 428 //_________________________________________________________________________________
 429 Int_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];
 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   }
 441   return Identify(combpid);
 442 }
 443 //__________________________________
 444 Bool_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
 451   AliDebug(2,Form("input particle: %i",iPsel));
 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 //___________________________________________
 465 Int_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   //
 471
 472   Int_t iPart = -1;
 473
 474   AliDebug(2,Form("calc response bef: %f  %f  %f  %f  %f",pid[0],pid[1],pid[2],pid[3],pid[4]));
 475   AliDebug(2,Form("priors           : %f  %f  %f  %f  %f",fPriors[0],fPriors[1],fPriors[2],fPriors[3],fPriors[4]));
 476
 477   AliPID getpid(pid,kTRUE);
 478   getpid.SetPriors(fPriors);
 479
 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]));
 484     if(fIsQAOn) fhCombProb[iP]->Fill(probability[iP]);
 485   }
 486
 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   }
 496
 497   if(fCheckResponse && !Check(pid,iPart, fMinDiffResponse)) iPart=kCheckResp;
 498   if(fCheckSelection && !Check(probability,iPart,fMinDiffProbability)) iPart=kCheckProb;
 499
 500   return iPart;
 501 }
 502 //___________________________________________
 503 Int_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);
 514   getpid.SetPriors(fPriors);
 515
 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 //___________________________________________
 535 Bool_t AliCFTrackCutPid::IsSelected(TObject *track){
 536   //
 537   //  method for the pid-cut selction
 538   //
 539   Bool_t sel = kFALSE;
 540
 541   if (!track) return kFALSE ;
 542   TString className(track->ClassName());
 543   if (className.CompareTo("AliESDtrack") == 0) {
 544   AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
 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);
 549   if(GetID(status,pid)==fgParticleType) sel = kTRUE;
 550   }
 551
 552   if (className.CompareTo("AliAODTrack") == 0) {
 553   AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
 554   if(GetAODID(aodtrack) == fgParticleType) sel = kTRUE;
 555   }
 556
 557  return sel;
 558
 559 }
 560 //__________________________________
 561 void  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;
 569   Double_t prod[AliPID::kSPECIES]={1.,1.,1.,1.,1.};
 570
 571   ULong_t andstatus =0;
 572   if(fIsDetAND) {
 573     andstatus = StatusForAND(status);
 574     AliDebug(1,Form("AND combination %lu",andstatus));
 575   }
 576
 577   //Products of single detector responses
 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));
 589             AliDebug(1,Form("In det %i ->  particle %i response is %f",i,j,pid[i][j]));
 590           }
 591         } else {
 592           prod[j]*=pid[i][j];
 593           isdet=kTRUE;
 594           AliDebug(2,Form("In det %i ->  particle %i response is %f",i,j,pid[i][j]));
 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            }
 607         }
 608       }//combined mode
 609     }//loop on dets
 610   }//loop on species
 611
 612
 613    //no detectors found, then go to ESD pid...
 614   if(!isdet) {
 615     AliWarning("\n !! No detector found for the combined pid --> responses are from the ESDpid !! \n");
 616     Double_t sumesdpid=0;
 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     }
 632
 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
 650 }
 651 //__________________________________________
 652 //
 653 //QA part
 654 //_________________________________________
 655 void 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 //______________________________________________
 670 void AliCFTrackCutPid::DefineHistograms()
 671 {
 672   //
 673   //QA histo booking
 674   //
 675
 676   if(fgIsAOD){
 677      const char *partic[AliPID::kSPECIES]={"electron","muon","pion","kaon","proton"};
 678
 679      for(Int_t iPart =0; iPart < AliPID::kSPECIES; iPart++)
 680       {
 681        fhCombResp[iPart] = new TH1F(Form("%s_rCombPart%i",GetName(),iPart),Form(" %s combined response  (AODTrack)  ",partic[iPart]),fNbins,fXmin,fXmax);
 682        fhCombResp[iPart]->SetXTitle(Form(" %s combined response ",partic[iPart]));
 683        fhCombResp[iPart]->SetYTitle("entries");
 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]));
 687        fhCombProb[iPart]->SetYTitle("entries");
 688        AliDebug(1,Form(  "%s is booked!!",fhCombProb[iPart]->GetName()));
 689      }
 690    }
 691
 692
 693   else {
 694    const char *detect[kNdets]={"ITS","TPC","TRD","TOF","HMPID"};
 695    const char *partic[AliPID::kSPECIES]={"electron","muon","pion","kaon","proton"};
 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++){
 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);
 702         fhResp[iDet][iP]->SetXTitle(Form(" %s response ",partic[iP]));
 703         fhResp[iDet][iP]->SetYTitle("entries");
 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]));
 706         fhProb[iDet][iP]->SetYTitle("entries");
 707        }
 708      }
 709
 710
 711   if(fgIsComb)
 712      {
 713       for(Int_t iPart =0; iPart < AliPID::kSPECIES; iPart++)
 714         {
 715           fhCombResp[iPart] = new TH1F(Form("%s_rCombPart%i",GetName(),iPart),Form(" %s combined response    ",partic[iPart]),fNbins,fXmin,fXmax);
 716           fhCombResp[iPart]->SetXTitle(Form(" %s response ",partic[iPart]));
 717           fhCombResp[iPart]->SetYTitle("entries");
 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);
 720           fhCombProb[iPart]->SetXTitle(Form(" %s response ",partic[iPart]));
 721           fhCombProb[iPart]->SetYTitle("entries");
 722           AliDebug(1,Form(  "%s is booked!!",fhCombProb[iPart]->GetName()));
 723         }
 724      }
 725    }
 726
 727 }
 728 //___________________________________________________
 729
 730 void AliCFTrackCutPid::AddQAHistograms(TList *qalist)
 731 {
 732   //
 733   // adds QA histograms in a TList
 734   //
 735   if(!fIsQAOn) return;
 736   DefineHistograms();
 737
 738   if(fgIsComb || fgIsAOD){
 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++){
 746     if(!fDets[iDet]) continue;
 747     for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){
 748       qalist->Add(fhResp[iDet][iP]);
 749       qalist->Add(fhProb[iDet][iP]);
 750     }
 751   }
 752
 753 }
 754