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 #include "AliCFTrackCutPid.h"
  60 #include "AliLog.h"
  61 #include <TMath.h>
  62 #include <TList.h>
  63
  64 ClassImp(AliCFTrackCutPid)
  65
  66 //__________________________________________________________________________________
  67 // CUT ON TRACK PID
  68 //__________________________________________________________________________________
  69 AliCFTrackCutPid::AliCFTrackCutPid() :
  70   AliCFCutBase(),
  71   fCut(0.2),
  72   fMinDiffResponse(0.0001),
  73   fMinDiffProbability(0.001),
  74   fgParticleType(10),
  75   fgIsComb(kTRUE),
  76   fCheckResponse(kFALSE),
  77   fCheckSelection(kTRUE),
  78   fIsPpriors(kFALSE),
  79   fIsDetAND(kFALSE),
  80   fXmin(-0.005),
  81   fXmax(1.005),
  82   fNbins(101),
  83   fDetRestr(-1),
  84   fiPartRestr(-1),
  85   fDetProbRestr(1)
  86 {
  87   //
  88   //Default constructor
  89   //
  90   for(Int_t j=0; j< AliPID::kSPECIES; j++) {
  91     fPriors[j]=0.2;
  92     fPriorsFunc[j]=0x0;
  93   }
  94
  95   for(Int_t jDet=0; jDet< kNdets; jDet++)  {
  96     fDets[jDet]=kFALSE;
  97     fDetsInAnd[jDet]=kFALSE;
  98   }
  99
 100   InitialiseHisto();
 101 }
 102 //______________________________
 103 AliCFTrackCutPid::AliCFTrackCutPid(const Char_t* name, const Char_t* title) :
 104   AliCFCutBase(name,title),
 105   fCut(0.2),
 106   fMinDiffResponse(0.0001),
 107   fMinDiffProbability(0.001),
 108   fgParticleType(10),
 109   fgIsComb(kTRUE),
 110   fCheckResponse(kFALSE),
 111   fCheckSelection(kTRUE),
 112   fIsPpriors(kFALSE),
 113   fIsDetAND(kFALSE),
 114   fXmin(-0.005),
 115   fXmax(1.005),
 116   fNbins(101),
 117   fDetRestr(-1),
 118   fiPartRestr(-1),
 119   fDetProbRestr(1)
 120 {
 121   //
 122   //Constructor
 123   //
 124   for(Int_t j=0; j< AliPID::kSPECIES; j++) {
 125     fPriors[j]=0.2;
 126     fPriorsFunc[j]=0x0;
 127   }
 128
 129   for(Int_t jDet=0; jDet< kNdets; jDet++)  {
 130     fDets[jDet]=kFALSE;
 131     fDetsInAnd[jDet]=kFALSE;
 132   }
 133
 134   InitialiseHisto();
 135 }
 136 //______________________________
 137 AliCFTrackCutPid::AliCFTrackCutPid(const AliCFTrackCutPid& c) :
 138   AliCFCutBase(c),
 139   fCut(c.fCut),
 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),
 148   fXmin(c.fXmin),
 149   fXmax(c.fXmax),
 150   fNbins(c.fNbins),
 151   fDetRestr(c.fDetRestr),
 152   fiPartRestr(c.fiPartRestr),
 153   fDetProbRestr(c.fDetProbRestr)
 154 {
 155   //
 156   //Copy constructor
 157   //
 158   for(Int_t i=0; i< kNdets ; i++ ) {
 159     fDets[i]=c.fDets[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];
 164     }
 165   }
 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];
 171   }
 172 }
 173 //______________________________
 174 AliCFTrackCutPid& AliCFTrackCutPid::operator=(const AliCFTrackCutPid& c)
 175 {
 176   //
 177   // Assignment operator
 178   //
 179   if (this != &c) {
 180     AliCFCutBase::operator=(c) ;
 181     this->fCut=c.fCut;
 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;
 190     this->fXmin=c.fXmin;
 191     this->fXmax=c.fXmax;
 192     this->fNbins=c.fNbins;
 193     this->fDetRestr=c.fDetRestr;
 194     this->fiPartRestr=c.fiPartRestr;
 195     this->fDetProbRestr=c.fDetProbRestr;
 196
 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];
 202       }
 203     }
 204
 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];
 210
 211     }
 212   }
 213   return *this ;
 214 }
 215 //__________________________________________________________________________________
 216 AliCFTrackCutPid::~AliCFTrackCutPid() {
 217   //
 218   //dtor
 219   //
 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];
 224     }
 225   }
 226
 227   for(Int_t j=0; j< AliPID::kSPECIES; j++){
 228     if(fhCombResp[j])delete fhCombResp[j];
 229     if(fhCombProb[j])delete fhCombProb[j];
 230
 231   }
 232
 233 }
 234 //__________________________________
 235 void AliCFTrackCutPid::SetDetectors(TString dets)
 236 {
 237   //
 238   // The string of detectors is translated into
 239   // the respective booelan data members
 240
 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;}
 246 }
 247 //__________________________________
 248 void AliCFTrackCutPid::SetPriors(Double_t r[AliPID::kSPECIES])
 249 {
 250   //
 251   // Sets the a priori concentrations
 252   //
 253   for(Int_t i=0; i<AliPID::kSPECIES; i++) fPriors[i]=r[i];
 254 }
 255 //__________________________________
 256 void AliCFTrackCutPid::SetPriorFunctions(TF1 *func[AliPID::kSPECIES])
 257 {
 258   //
 259   // Sets the momentu dependent a priori concentrations
 260   //
 261
 262   for(Int_t i=0; i<AliPID::kSPECIES; i++) fPriorsFunc[i]=func[i];
 263   fIsPpriors = kTRUE;
 264 }
 265 //____________________________________________
 266 void AliCFTrackCutPid::SetANDstatus(TString dets)
 267 {
 268   //
 269   //Sets the detectors to be in AND for the combined PID
 270   //
 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;}
 276
 277   fIsDetAND = kTRUE;
 278 }
 279 //
 280 void AliCFTrackCutPid::SetDetectorProbabilityRestriction(TString det, Int_t iPart, Double_t upperprob)
 281 {
 282   //
 283   // Sets the detector, the particle and the probability
 284   // limit.
 285   //
 286
 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;
 292   fiPartRestr = iPart;
 293   fDetProbRestr = upperprob;
 294 }
 295 //__________________________________
 296 void AliCFTrackCutPid::TrackInfo(const AliESDtrack *pTrk, ULong_t status[kNdets+1],Double_t pid[kNdets+1][AliPID::kSPECIES]) const
 297 {
 298   //
 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.
 301   //
 302   if(fDets[kITS]) {
 303     pTrk->GetITSpid(pid[kITS]);
 304     status[kITS]=AliESDtrack::kITSpid;
 305   }
 306   if(fDets[kTPC]) {
 307     pTrk->GetTPCpid(pid[kTPC]);
 308     status[kTPC]=AliESDtrack::kTPCpid;
 309   }
 310   if(fDets[kTRD]) {
 311     pTrk->GetTRDpid(pid[kTRD]);
 312     status[kTRD]=AliESDtrack::kTRDpid;
 313   }
 314   if(fDets[kTOF]) {
 315     pTrk->GetTOFpid(pid[kTOF]);
 316     status[kTOF]=AliESDtrack::kTOFpid;
 317   }
 318   if(fDets[kHMPID]) {
 319     pTrk->GetHMPIDpid(pid[kHMPID]);
 320     status[kHMPID]=AliESDtrack::kHMPIDpid;
 321   }
 322   if(fDetRestr>=0){
 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]);
 328   }
 329
 330   status[kNdets]=pTrk->GetStatus();
 331   pTrk->GetESDpid(pid[kNdets]);
 332 }
 333 //__________________________________
 334 void AliCFTrackCutPid::SetPPriors(AliESDtrack *pTrk)
 335 {
 336   //
 337   //sets the mommentum dependent a priori concentrations
 338   //
 339
 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;}
 343   }
 344 }
 345 //______________________________________
 346 ULong_t AliCFTrackCutPid::StatusForAND(ULong_t status[kNdets+1]) const
 347 {
 348   //
 349   //In case of AND of more detectors the AND-detector status combination.
 350   //is calculated and also returned
 351   //
 352
 353   ULong_t andstatus=0;
 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));
 358   }
 359   return andstatus;
 360 }
 361 //_______________________________________
 362 Int_t AliCFTrackCutPid::GetID(ULong_t status[kNdets+1],Double_t pid[kNdets+1][AliPID::kSPECIES]) const
 363 {
 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
 370
 371   Int_t iPart=-1;
 372
 373   if(!fgIsComb){
 374     Bool_t isDet=kFALSE;
 375     for(Int_t i=0; i<kNdets; i++){
 376       if(!fDets[i]) continue;
 377       isDet=kTRUE;
 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])){
 380         iPart=-10;
 381         AliDebug(1,Form("detector %i -> pid trk status not ok",i));
 382       }
 383       else {
 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]);
 387       }
 388     }
 389     if(!isDet){
 390       AliDebug(1,Form("  !!!        No detector selected, the ESD-pid response is considered"));
 391       iPart = Identify(pid[kNdets]);
 392     }
 393   }else{
 394     Double_t calcprob[5];
 395     CombPID(status,pid,calcprob);
 396     iPart = Identify(calcprob);
 397   }
 398
 399
 400   AliDebug(1,Form("selected particle: %i",iPart));
 401
 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) {
 407       iPart = kDetRestr;
 408       AliDebug(1,"\n\n the detector restrictions refused the ID \n\n");
 409     }
 410   }//cross checks with one detector probability
 411
 412   AliDebug(1,Form("after the check the selected particle is %i",iPart));
 413
 414   return iPart;
 415 }
 416 //__________________________________
 417 Bool_t AliCFTrackCutPid::Check(const Double_t *p, Int_t iPsel, Double_t minDiff) const
 418 {
 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
 423
 424   AliDebug(1,Form("input particle: %i",iPsel));
 425   Bool_t ck=kTRUE;
 426
 427   if(iPsel<0) ck=kFALSE;
 428
 429   else {
 430     for(Int_t j=0; j< AliPID::kSPECIES; j++) {
 431       if(j!=iPsel && TMath::Abs(p[j]-p[iPsel])<minDiff) ck=kFALSE;
 432     }
 433     if(!ck) AliDebug(1,"the values are too close ");
 434   }
 435   return ck;
 436 }
 437 //___________________________________________
 438 Int_t AliCFTrackCutPid::Identify(Double_t pid[AliPID::kSPECIES]) const
 439 {
 440   //
 441   // The identification is actually performed here with possible
 442   // checks on the det responses and/or probabilities
 443   //
 444   Int_t iPart = -1;
 445
 446   AliPID getpid(pid,kTRUE);
 447
 448   if(fgIsComb) {
 449     Double_t priors[5]={0.2,0.2,0.2,0.2,0.2};
 450     getpid.SetPriors(priors);
 451   }
 452   else getpid.SetPriors(fPriors);
 453
 454
 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);
 458
 459
 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;
 465   return iPart;
 466 }
 467 //___________________________________________
 468 Int_t AliCFTrackCutPid::IdentifyQA(const Double_t pid[AliPID::kSPECIES], Int_t idets) const
 469 {
 470   //
 471   // The same as Identify, but with the QA histo filling
 472   //
 473   //
 474
 475   Int_t iPart = -1;
 476   AliDebug(1,Form("resp   : %f  %f  %f  %f  %f",pid[0],pid[1],pid[2],pid[3],pid[4]));
 477
 478   AliPID getpid(pid,kTRUE);
 479
 480   if(fgIsComb) {
 481     Double_t priors[5]={0.2,0.2,0.2,0.2,0.2};
 482     getpid.SetPriors(priors);
 483   }
 484   else getpid.SetPriors(fPriors);
 485
 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]);
 491   }
 492
 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));
 496
 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;
 502   return iPart;
 503 }
 504 //___________________________________________
 505 Bool_t AliCFTrackCutPid::IsSelected(TObject *track){
 506   //
 507   //  method for the pid-cut selction
 508   //
 509
 510   if (!track) return kFALSE ;
 511   TString className(track->ClassName());
 512   if (className.CompareTo("AliESDtrack") != 0) {
 513     AliError("obj must point to a AliESDtrack ");
 514     return kFALSE ;
 515   }
 516
 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;
 523   else return kFALSE;
 524 }
 525 //__________________________________
 526 void  AliCFTrackCutPid::CombPID(ULong_t status[kNdets+1],Double_t pid[kNdets+1][AliPID::kSPECIES],Double_t *combpid) const
 527 {
 528   //
 529   // Calculates the combined PID according to the chosen detectors.
 530   // and provides the array of probabilities
 531   //
 532
 533   Bool_t isdet=kFALSE;
 534   Double_t sum=0.;
 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;
 539   if(fIsDetAND) {
 540     andstatus = StatusForAND(status);
 541     AliDebug(1,Form("AND combination %lu",andstatus));
 542   }
 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]) {
 547         if(fIsDetAND) {
 548           ULong_t checkstatus = status[kNdets]&andstatus;
 549           if(checkstatus != andstatus) continue;
 550           else {
 551             prod[j]*=pid[i][j];
 552             isdet = kTRUE;
 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]));
 555           }
 556         }
 557         else {
 558           prod[j]*=pid[i][j];
 559           isdet=kTRUE;
 560           AliDebug(1,Form("In det loop %i ->  particle %i response is %f",i,j,pid[i][j]));
 561         }
 562       }//combined mode
 563     }//loop on dets
 564   }//loop on species
 565
 566   if(fIsQAOn) {
 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));
 574       }//loop on part
 575     }//loop on dets
 576   }//if qa
 577
 578   if(!isdet) {
 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;}
 581   }
 582
 583   else{
 584     for(Int_t k=0; k<AliPID::kSPECIES; k++){
 585       if(fIsQAOn) {
 586         if(!fhCombResp[k]) AliDebug(1,Form("no fhCombResp[%i], check if pidcut->Init() was called",k));
 587         else fhCombResp[k]->Fill(prod[k]);
 588       }
 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]));
 592       sum+=comb[k];
 593     }
 594
 595     if(sum == 0) {
 596       AliDebug(1,"Check the detector responses or the priors, their combined products are zero");
 597       return;
 598     }
 599     for(Int_t n=0; n<AliPID::kSPECIES; n++) {
 600       combpid[n]=fPriors[n]*prod[n]/sum;
 601       if(fIsQAOn) {
 602         if(!fhCombProb[n]) Printf("no fhCombResp[%i], check if pidcut->Init() was called",n);
 603         fhCombProb[n]->Fill(combpid[n]);
 604       }
 605     }
 606   }
 607 }
 608 //__________________________________________
 609 //
 610 //QA part
 611 //_________________________________________
 612 void AliCFTrackCutPid::InitialiseHisto()
 613 {
 614   //
 615   //QA histo initialization
 616   //
 617   for(Int_t iP=0; iP<AliPID::kSPECIES; iP++){
 618     fhCombResp[iP]=0x0;
 619     fhCombProb[iP]=0x0;
 620     for(Int_t iDet =0; iDet<kNdets; iDet++){
 621       fhResp[iDet][iP]=0x0;
 622       fhProb[iDet][iP]=0x0;
 623     }
 624   }
 625 }
 626 //______________________________________________
 627 void AliCFTrackCutPid::Init()
 628 {
 629   //
 630   // initialises QA histograms
 631   //
 632
 633   if(fIsQAOn) DefineHistograms();
 634 }
 635 //_________________________________________________
 636 void AliCFTrackCutPid::DefineHistograms()
 637 {
 638   //
 639   //QA histo booking
 640   //
 641   char *detect[5]={"ITS","TPC","TRD","TOF","HMPID"};
 642   char *partic[5]={"electron","muon","pion","kaon","proton"};
 643
 644   for(Int_t iDet =0; iDet< kNdets; iDet++)
 645     {
 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);
 650       }
 651     }
 652
 653   if(fgIsComb){
 654     for(Int_t iPart =0; iPart < AliPID::kSPECIES; iPart++)
 655       {
 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);
 658       }
 659   }
 660 }
 661 //___________________________________________________
 662
 663 void AliCFTrackCutPid::AddQAHistograms(TList *qalist) const
 664 {
 665   //
 666   // adds QA histograms in a TList
 667   //
 668   if(!fIsQAOn) return;
 669
 670   if(fgIsComb){
 671     for(Int_t iPart =0; iPart<AliPID::kSPECIES; iPart++){
 672       qalist->Add(fhCombResp[iPart]);
 673       qalist->Add(fhCombProb[iPart]);
 674     }
 675   }
 676
 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]);
 681     }
 682   }
 683 }