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