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
 221   for(Int_t i=0; i< kNdets ; i++ )   {
 222     for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){
 223       if(fhResp[i][iP])delete fhResp[i][iP];
 224       if(fhProb[i][iP])delete fhProb[i][iP];
 225     }
 226   }
 227
 228   for(Int_t j=0; j< AliPID::kSPECIES; j++){
 229     if(fhCombResp[j])delete fhCombResp[j];
 230     if(fhCombProb[j])delete fhCombProb[j];
 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   if(dets.Contains("ALL")) for(Int_t i=0; i< kNdets ; i++) fDets[i]=kTRUE;
 248 }
 249 //__________________________________
 250 void AliCFTrackCutPid::SetPriors(Double_t r[AliPID::kSPECIES])
 251 {
 252   //
 253   // Sets the a priori concentrations
 254   //
 255   for(Int_t i=0; i<AliPID::kSPECIES; i++) fPriors[i]=r[i];
 256 }
 257 //__________________________________
 258 void AliCFTrackCutPid::SetPriorFunctions(TF1 *func[AliPID::kSPECIES])
 259 {
 260   //
 261   // Sets the momentu dependent a priori concentrations
 262   //
 263
 264   for(Int_t i=0; i<AliPID::kSPECIES; i++) fPriorsFunc[i]=func[i];
 265   fIsPpriors = kTRUE;
 266 }
 267 //____________________________________________
 268 void AliCFTrackCutPid::SetANDstatus(TString dets)
 269 {
 270   //
 271   //Sets the detectors to be in AND for the combined PID
 272   //
 273   if(dets.Contains("ITS") && fDets[kITS])     {fDetsInAnd[kITS]=kTRUE;}
 274   if(dets.Contains("TPC") && fDets[kTPC])     {fDetsInAnd[kTPC]=kTRUE;}
 275   if(dets.Contains("TRD") && fDets[kTRD])     {fDetsInAnd[kTRD]=kTRUE;}
 276   if(dets.Contains("TOF") && fDets[kTOF])     {fDetsInAnd[kTOF]=kTRUE;}
 277   if(dets.Contains("HMPID") && fDets[kHMPID]) {fDetsInAnd[kHMPID]=kTRUE;}
 278
 279   fIsDetAND = kTRUE;
 280 }
 281 //
 282 void AliCFTrackCutPid::SetDetectorProbabilityRestriction(TString det, Int_t iPart, Double_t upperprob)
 283 {
 284   //
 285   // Sets the detector, the particle and the probability
 286   // limit.
 287   //
 288
 289   if(det.Contains("ITS")) fDetRestr = kITS;
 290   if(det.Contains("TPC")) fDetRestr = kTPC;
 291   if(det.Contains("TRD")) fDetRestr = kTRD;
 292   if(det.Contains("TOF")) fDetRestr = kTOF;
 293   if(det.Contains("HMPID")) fDetRestr = kHMPID;
 294   fiPartRestr = iPart;
 295   fDetProbRestr = upperprob;
 296 }
 297 //__________________________________
 298 void AliCFTrackCutPid::TrackInfo(const AliESDtrack *pTrk, ULong_t status[kNdets+1],Double_t pid[kNdets+1][AliPID::kSPECIES]) const
 299 {
 300   //
 301   // Loads the responses of the XXX chosen detectors for the pTrk
 302   // and the corresponding trk status. The final trk status is also loaded.
 303   //
 304   if(fDets[kITS]) {
 305     pTrk->GetITSpid(pid[kITS]);
 306     status[kITS]=AliESDtrack::kITSpid;
 307   }
 308   if(fDets[kTPC]) {
 309     pTrk->GetTPCpid(pid[kTPC]);
 310     status[kTPC]=AliESDtrack::kTPCpid;
 311   }
 312   if(fDets[kTRD]) {
 313     pTrk->GetTRDpid(pid[kTRD]);
 314     status[kTRD]=AliESDtrack::kTRDpid;
 315   }
 316   if(fDets[kTOF]) {
 317     pTrk->GetTOFpid(pid[kTOF]);
 318     status[kTOF]=AliESDtrack::kTOFpid;
 319   }
 320   if(fDets[kHMPID]) {
 321     pTrk->GetHMPIDpid(pid[kHMPID]);
 322     status[kHMPID]=AliESDtrack::kHMPIDpid;
 323   }
 324   if(fDetRestr>=0){
 325     if(fDetRestr == kITS) pTrk->GetITSpid(pid[kITS]);
 326     if(fDetRestr == kTPC) pTrk->GetITSpid(pid[kTPC]);
 327     if(fDetRestr == kTRD) pTrk->GetITSpid(pid[kTRD]);
 328     if(fDetRestr == kTOF) pTrk->GetITSpid(pid[kTOF]);
 329     if(fDetRestr == kHMPID) pTrk->GetITSpid(pid[kHMPID]);
 330   }
 331
 332   status[kNdets]=pTrk->GetStatus();
 333   pTrk->GetESDpid(pid[kNdets]);
 334 }
 335 //__________________________________
 336 void AliCFTrackCutPid::SetPPriors(AliESDtrack *pTrk)
 337 {
 338   //
 339   //sets the mommentum dependent a priori concentrations
 340   //
 341
 342   for(Int_t i=0; i< AliPID::kSPECIES; i++) {
 343     if(pTrk->P()>fPriorsFunc[i]->GetXmin() && pTrk->P() < fPriorsFunc[i]->GetXmax()) fPriors[i]=fPriorsFunc[i]->Eval(pTrk->P());
 344     else {AliInfo("the track momentum is not in the function range. Priors are equal") fPriors[i] = 0.2;}
 345   }
 346 }
 347 //______________________________________
 348 ULong_t AliCFTrackCutPid::StatusForAND(ULong_t status[kNdets+1]) const
 349 {
 350   //
 351   //In case of AND of more detectors the AND-detector status combination.
 352   //is calculated and also returned
 353   //
 354
 355   ULong_t andstatus=0;
 356   for(Int_t i=0; i< kNdets; i++) {
 357     if(!fDetsInAnd[i]) continue;
 358     andstatus = andstatus | status[i];
 359     AliDebug(1,Form("trk status %lu  %i AND-status  combination: %lu",status[i],i,andstatus));
 360   }
 361   return andstatus;
 362 }
 363 //_______________________________________
 364 Int_t AliCFTrackCutPid::GetID(ULong_t status[kNdets+1],Double_t pid[kNdets+1][AliPID::kSPECIES]) const
 365 {
 366   // Identifies the track if its probability is higher than the cut
 367   // value. The default value is fCut=0.2, therefore the most probable
 368   // one is identified by default. Here all the checks on how to identify
 369   // the track are performed (single detector or combined PID,..., detector
 370   // restriction probability
 371   // Returns:   integer corresponding to the identified particle
 372
 373   Int_t iPart=-1;
 374
 375   if(!fgIsComb){
 376     Bool_t isDet=kFALSE;
 377     for(Int_t i=0; i<kNdets; i++){
 378       if(!fDets[i]) continue;
 379       isDet=kTRUE;
 380       AliDebug(1,Form("trk status %lu  %i-det-pid status %lu   -> combination: %lu",status[kNdets],i,status[i],status[kNdets]&status[i]));
 381       if(!(status[kNdets]&status[i])){
 382         iPart=-10;
 383         AliDebug(1,Form("detector %i -> pid trk status not ok",i));
 384       }
 385       else {
 386         AliDebug(1,Form("resp   : %f  %f  %f  %f  %f",pid[i][0],pid[i][1],pid[i][2],pid[i][3],pid[i][4]));
 387         if(fIsQAOn) iPart = IdentifyQA(pid[i],i);
 388         else iPart = Identify(pid[i]);
 389       }
 390     }
 391     if(!isDet){
 392       AliDebug(1,Form("  !!!        No detector selected, the ESD-pid response is considered"));
 393       iPart = Identify(pid[kNdets]);
 394     }
 395   }else{
 396     Double_t calcprob[5];
 397     CombPID(status,pid,calcprob);
 398     iPart = Identify(calcprob);
 399   }
 400
 401
 402   AliDebug(1,Form("selected particle: %i",iPart));
 403
 404   if(iPart >=0 && fiPartRestr>=0) {
 405     AliPID restr(pid[fDetRestr]);
 406     restr.SetPriors(fPriors);
 407     AliDebug(1,Form("setted upper limit: %f  det %i : probability %f ",fDetProbRestr,fDetRestr,restr.GetProbability((AliPID::EParticleType)fiPartRestr)));
 408     if(restr.GetProbability((AliPID::EParticleType)fiPartRestr) > fDetProbRestr) {
 409       iPart = kDetRestr;
 410       AliDebug(1,"\n\n the detector restrictions refused the ID \n\n");
 411     }
 412   }//cross checks with one detector probability
 413
 414   AliDebug(1,Form("after the check the selected particle is %i",iPart));
 415
 416   return iPart;
 417 }
 418 //__________________________________
 419 Bool_t AliCFTrackCutPid::Check(const Double_t *p, Int_t iPsel, Double_t minDiff) const
 420 {
 421   // Checks if there are no equal values and if the valus
 422   // difference between the selected particle and the others i
 423   // is higher than  a lower limit.
 424   // Returns:  kTRUE= is acceptable
 425
 426   AliDebug(2,Form("input particle: %i",iPsel));
 427   Bool_t ck=kTRUE;
 428
 429   if(iPsel<0) ck=kFALSE;
 430
 431   else {
 432     for(Int_t j=0; j< AliPID::kSPECIES; j++) {
 433       if(j!=iPsel && TMath::Abs(p[j]-p[iPsel])<minDiff) ck=kFALSE;
 434     }
 435     if(!ck) AliDebug(1,"the values are too close ");
 436   }
 437   return ck;
 438 }
 439 //___________________________________________
 440 Int_t AliCFTrackCutPid::Identify(Double_t pid[AliPID::kSPECIES]) const
 441 {
 442   //
 443   // The identification is actually performed here with possible
 444   // checks on the det responses and/or probabilities
 445   //
 446   Int_t iPart = -1;
 447
 448   AliPID getpid(pid,kTRUE);
 449
 450   if(fgIsComb) {
 451     Double_t priors[5]={0.2,0.2,0.2,0.2,0.2};
 452     getpid.SetPriors(priors);
 453   }
 454   else getpid.SetPriors(fPriors);
 455
 456
 457   AliPID::EParticleType sel = getpid.GetMostProbable();
 458   Double_t probability[AliPID::kSPECIES];
 459   for(Int_t iP=0; iP<AliPID::kSPECIES; iP++) probability[iP] = getpid.GetProbability((AliPID::EParticleType)iP);
 460
 461
 462   if(getpid.GetProbability(sel,fPriors)>fCut) iPart= (Int_t)sel;
 463   AliDebug(2,Form("calc response   : %f  %f  %f  %f  %f",pid[0],pid[1],pid[2],pid[3],pid[4]));
 464   AliDebug(2,Form("probabilities   : %f  %f  %f  %f  %f",probability[0],probability[1],probability[2],probability[3],probability[4]));
 465
 466   if(fCheckResponse && !Check(pid,iPart, fMinDiffResponse)) iPart=kCheckResp;
 467   if(fCheckSelection && !Check(probability,iPart,fMinDiffProbability)) iPart=kCheckProb;
 468
 469   return iPart;
 470 }
 471 //___________________________________________
 472 Int_t AliCFTrackCutPid::IdentifyQA(const Double_t pid[AliPID::kSPECIES], Int_t idets) const
 473 {
 474   //
 475   // The same as Identify, but with the QA histo filling
 476   //
 477   //
 478
 479   Int_t iPart = -1;
 480   AliDebug(1,Form("resp   : %f  %f  %f  %f  %f",pid[0],pid[1],pid[2],pid[3],pid[4]));
 481
 482   AliPID getpid(pid,kTRUE);
 483
 484   if(fgIsComb) {
 485     Double_t priors[5]={0.2,0.2,0.2,0.2,0.2};
 486     getpid.SetPriors(priors);
 487   }
 488   else getpid.SetPriors(fPriors);
 489
 490   AliPID::EParticleType sel = getpid.GetMostProbable();
 491   Double_t probability[AliPID::kSPECIES];
 492   for(Int_t iP=0; iP<AliPID::kSPECIES; iP++) {
 493     probability[iP] = getpid.GetProbability((AliPID::EParticleType)iP);
 494     fhProb[idets][iP]->Fill(probability[iP]);
 495   }
 496
 497   AliPID toresp(pid,kTRUE); Double_t qapriors[5]={0.2,0.2,0.2,0.2,0.2};
 498   toresp.SetPriors(qapriors);
 499   for(Int_t iPr=0; iPr<AliPID::kSPECIES; iPr++) fhResp[idets][iPr]->Fill(toresp.GetProbability((AliPID::EParticleType)iPr));
 500
 501   if(getpid.GetProbability(sel,fPriors)>fCut) iPart= (Int_t)sel;
 502   AliDebug(1,Form("resp   : %f  %f  %f  %f  %f",pid[0],pid[1],pid[2],pid[3],pid[4]));
 503   AliDebug(1,Form("probab : %f  %f  %f  %f  %f",probability[0],probability[1],probability[2],probability[3],probability[4]));
 504   if(fCheckResponse && !Check(pid,iPart, fMinDiffResponse)) iPart=kCheckResp;
 505   if(fCheckSelection && !Check(probability,iPart,fMinDiffProbability)) iPart=kCheckProb;
 506   return iPart;
 507 }
 508 //___________________________________________
 509 Bool_t AliCFTrackCutPid::IsSelected(TObject *track){
 510   //
 511   //  method for the pid-cut selction
 512   //
 513
 514   if (!track) return kFALSE ;
 515   TString className(track->ClassName());
 516   if (className.CompareTo("AliESDtrack") != 0) {
 517     AliError("obj must point to a AliESDtrack ");
 518     return kFALSE ;
 519   }
 520
 521   AliESDtrack *esdTrack = (AliESDtrack *)track;
 522   ULong_t status[kNdets+1]={0,0,0,0,0,0};
 523   Double_t pid[kNdets+1][AliPID::kSPECIES];
 524   TrackInfo(esdTrack,status,pid);
 525   if(fIsPpriors) SetPPriors(esdTrack);
 526   if(GetID(status,pid)==fgParticleType) return kTRUE;
 527   else return kFALSE;
 528 }
 529 //__________________________________
 530 void  AliCFTrackCutPid::CombPID(ULong_t status[kNdets+1],Double_t pid[kNdets+1][AliPID::kSPECIES],Double_t *combpid) const
 531 {
 532   //
 533   // Calculates the combined PID according to the chosen detectors.
 534   // and provides the array of probabilities
 535   //
 536
 537   Bool_t isdet=kFALSE;
 538   Double_t sum=0.;
 539   Double_t prod[AliPID::kSPECIES]={1.,1.,1.,1.,1.};
 540   Double_t comb[AliPID::kSPECIES]={0.,0.,0.,0.,0.};
 541   Double_t priors[AliPID::kSPECIES]={0.2,0.2,0.2,0.2,0.2};
 542
 543   ULong_t andstatus =0;
 544   if(fIsDetAND) {
 545     andstatus = StatusForAND(status);
 546     AliDebug(1,Form("AND combination %lu",andstatus));
 547   }
 548   for(Int_t j=0; j<AliPID::kSPECIES; j++){
 549     for(Int_t i=0; i< kNdets; i++){
 550       if(!fDets[i]) continue;
 551       if(status[kNdets]&status[i]) {
 552         if(fIsDetAND) {
 553           ULong_t checkstatus = status[kNdets]&andstatus;
 554           if(checkstatus != andstatus) continue;
 555           else {
 556             prod[j]*=pid[i][j];
 557             isdet = kTRUE;
 558             AliDebug(1,Form("-----> trk status %lu   and status %lu -> trk-ANDdetector status combination %lu",status[kNdets],andstatus,status[kNdets]&andstatus));
 559             AliDebug(1,Form("In det %i ->  particle %i response is %f",i,j,pid[i][j]));
 560           }
 561         }
 562         else {
 563           prod[j]*=pid[i][j];
 564           isdet=kTRUE;
 565           AliDebug(2,Form("In det %i ->  particle %i response is %f",i,j,pid[i][j]));
 566         }
 567       }//combined mode
 568     }//loop on dets
 569   }//loop on species
 570
 571   if(fIsQAOn) {
 572     for(Int_t iqa =0; iqa < kNdets; iqa++){
 573       if(!fDets[iqa]) continue;
 574       AliPID normresp(pid[iqa]);
 575       normresp.SetPriors(priors);
 576       for(Int_t ip =0; ip< AliPID::kSPECIES; ip++){
 577         if(!fhResp[iqa][ip]) {AliDebug(1,Form("no pointer to the histo fhResp%i%i, check if pidcut->Init() was called",iqa,ip));}
 578         else fhResp[iqa][ip]->Fill(normresp.GetProbability((AliPID::EParticleType)ip));
 579       }//loop on part
 580     }//loop on dets
 581   }//if qa
 582
 583
 584   Double_t add = 0; for(Int_t isumm=0; isumm<5; isumm++) add+=prod[isumm];
 585   if(add>0) AliDebug(1,Form("calculated comb response: %f %f %f %f %f",prod[0]/add,prod[1]/add,prod[2]/add,prod[3]/add,prod[4]/add));
 586             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]));
 587
 588   if(!isdet) {
 589     AliWarning("\n !! No detector found for the combined pid --> responses are from the ESDpid !! \n");
 590     Double_t sumesdpid=0;
 591     for(Int_t nn=0; nn<AliPID::kSPECIES; nn++) {
 592       combpid[nn]=pid[kNdets][nn];
 593       sumesdpid+=fPriors[nn]*combpid[nn];
 594       if(fIsQAOn) {
 595         if(!fhCombResp[nn]) AliDebug(1,Form("\n no fhCombResp[%i], check if pidcut->Init() was called",nn));
 596         else fhCombResp[nn]->Fill(combpid[nn]);
 597       }//fIsQAOn
 598     }//loop on species
 599     if(sumesdpid > 0) for(Int_t ih = 0; ih < AliPID::kSPECIES; ih++) fhCombProb[ih]->Fill(fPriors[ih]*combpid[ih]/sumesdpid);
 600     else AliDebug(1,"priors or ESDpid are zero, please check them");
 601   }// end no det
 602
 603   else{
 604     for(Int_t k=0; k<AliPID::kSPECIES; k++){
 605       if(fIsQAOn) {
 606         if(!fhCombResp[k]) AliDebug(1,Form("\n no fhCombResp[%i], check if pidcut->Init() was called \n",k));
 607         else fhCombResp[k]->Fill(prod[k]);
 608       }
 609       AliDebug(1,Form("species: %i priors %f and combined prod %f",k,fPriors[k],prod[k]));
 610       comb[k]=fPriors[k]*prod[k];
 611       AliDebug(1,Form("combined probability %i  %f",k,comb[k]));
 612       sum+=comb[k];
 613     }
 614
 615     if(sum == 0) {
 616       AliDebug(1,"Check the detector responses or the priors, their combined products are zero");
 617       return;
 618     }
 619     for(Int_t n=0; n<AliPID::kSPECIES; n++) {
 620       combpid[n]=fPriors[n]*prod[n]/sum;
 621       if(fIsQAOn) {
 622         if(!fhCombProb[n]) AliDebug(1,Form("no fhCombRespi[%i] defined, check if pidcut->Init() was called",n));
 623         else fhCombProb[n]->Fill(combpid[n]);
 624       }
 625     }
 626   }//end else
 627 }
 628 //__________________________________________
 629 //
 630 //QA part
 631 //_________________________________________
 632 void AliCFTrackCutPid::InitialiseHisto()
 633 {
 634   //
 635   //QA histo initialization
 636   //
 637   for(Int_t iP=0; iP<AliPID::kSPECIES; iP++){
 638     fhCombResp[iP]=0x0;
 639     fhCombProb[iP]=0x0;
 640     for(Int_t iDet =0; iDet<kNdets; iDet++){
 641       fhResp[iDet][iP]=0x0;
 642       fhProb[iDet][iP]=0x0;
 643     }
 644   }
 645 }
 646 //______________________________________________
 647 void AliCFTrackCutPid::Init()
 648 {
 649   //
 650   // initialises QA histograms
 651   //
 652
 653   if(fIsQAOn) DefineHistograms();
 654 }
 655 //_________________________________________________
 656 void AliCFTrackCutPid::DefineHistograms()
 657 {
 658   //
 659   //QA histo booking
 660   //
 661   char *detect[kNdets]={"ITS","TPC","TRD","TOF","HMPID"};
 662   char *partic[AliPID::kSPECIES]={"electron","muon","pion","kaon","proton"};
 663
 664    for(Int_t iDet =0; iDet< kNdets; iDet++)
 665     {
 666       if(!fDets[iDet]) continue;
 667       for(Int_t iP =0; iP < AliPID::kSPECIES; iP++){
 668         fhResp[iDet][iP] = new TH1F(Form("rDet%iPart%i",iDet,iP),Form("%s %s response    ",detect[iDet],partic[iP]),fNbins,fXmin,fXmax);
 669         fhProb[iDet][iP] = new TH1F(Form("pDet%iPart%i",iDet,iP),Form("%s %s probability ",detect[iDet],partic[iP]),fNbins,fXmin,fXmax);
 670       }
 671     }
 672
 673
 674 if(fgIsComb)
 675    {
 676     for(Int_t iPart =0; iPart < AliPID::kSPECIES; iPart++)
 677       {
 678         fhCombResp[iPart] = new TH1F(Form("rCombPart%i",iPart),Form(" %s combined response    ",partic[iPart]),fNbins,fXmin,fXmax);
 679         fhCombProb[iPart] = new TH1F(Form("pCombPart%i",iPart),Form("%s combined probability ",partic[iPart]),fNbins,fXmin,fXmax);
 680       }
 681   }
 682 }
 683 //___________________________________________________
 684
 685 void AliCFTrackCutPid::AddQAHistograms(TList *qalist) const
 686 {
 687   //
 688   // adds QA histograms in a TList
 689   //
 690   if(!fIsQAOn) return;
 691
 692   if(fgIsComb){
 693     for(Int_t iPart =0; iPart<AliPID::kSPECIES; iPart++){
 694       qalist->Add(fhCombResp[iPart]);
 695       qalist->Add(fhCombProb[iPart]);
 696     }
 697   }
 698
 699   for(Int_t iDet=0; iDet<kNdets; iDet++){
 700     if(!fDets[iDet]) continue;
 701     for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){
 702       qalist->Add(fhResp[iDet][iP]);
 703       qalist->Add(fhProb[iDet][iP]);
 704     }
 705   }
 706
 707 }