PWG/Tools/AliHelperPID.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 //-----------------------------------------------------------------
  18 //         AliHelperPID class
  19 //-----------------------------------------------------------------
  20
  21 #include "AliHelperPID.h"
  22 #include "AliVEvent.h"
  23 #include "AliAODEvent.h"
  24 #include "AliMCEvent.h"
  25 #include "AliMCEventHandler.h"
  26 #include "TH1F.h"
  27 #include "TH2F.h"
  28 #include "TList.h"
  29 #include "AliStack.h"
  30 #include "AliVTrack.h"
  31 #include "TParticle.h"
  32 #include "AliAODMCParticle.h"
  33 #include "AliPIDResponse.h"
  34 #include "AliPIDCombined.h"
  35 #include "AliAnalysisManager.h"
  36 #include "AliInputEventHandler.h"
  37
  38 using namespace AliHelperPIDNameSpace;
  39 using namespace std;
  40
  41 const char * kPIDTypeName[]={"TPC","TOF","TPC-TOF"} ;
  42 const char * kDetectorName[]={"ITS","TPC","TOF"} ;
  43 const char * kParticleSpeciesName[]={"Pions","Kaons","Protons","Undefined"} ;
  44
  45 ClassImp(AliHelperPID)
  46
  47 AliHelperPID::AliHelperPID() : TNamed("HelperPID", "PID object"),fisMC(0), fPIDType(kNSigmaTPCTOF), fNSigmaPID(3), fBayesCut(0.8), fPIDResponse(0x0), fPIDCombined(0x0),fOutputList(0x0),fRequestTOFPID(1),fRemoveTracksT0Fill(0),fUseExclusiveNSigma(0),fPtTOFPID(.6),fHasTOFPID(0){
  48
  49   for(Int_t ipart=0;ipart<kNSpecies;ipart++)
  50     for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++)
  51       fnsigmas[ipart][ipid]=999.;
  52
  53   for(Int_t ipart=0;ipart<kNSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;
  54
  55   fOutputList = new TList;
  56   fOutputList->SetOwner();
  57   fOutputList->SetName("OutputList");
  58
  59   //nsigma plot
  60   for(Int_t ipart=0;ipart<kNSpecies;ipart++){
  61     for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
  62       Double_t miny=-30;
  63       Double_t maxy=30;
  64       if(ipid==kNSigmaTPCTOF){miny=0;maxy=50;}
  65       TH2F *fHistoNSigma=new TH2F(Form("NSigma_%d_%d",ipart,ipid),Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
  66       fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
  67       fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
  68       fOutputList->Add(fHistoNSigma);
  69     }
  70   }
  71
  72   //nsigmaRec plot
  73   for(Int_t ipart=0;ipart<kNSpecies;ipart++){
  74     for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
  75       Double_t miny=-10;
  76       Double_t maxy=10;
  77       if(ipid==kNSigmaTPCTOF){miny=0;maxy=20;}
  78       TH2F *fHistoNSigma=new TH2F(Form("NSigmaRec_%d_%d",ipart,ipid),
  79                                   Form("n#sigma for reconstructed %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
  80       fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
  81       fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
  82       fOutputList->Add(fHistoNSigma);
  83     }
  84   }
  85
  86   //BayesRec plot
  87   for(Int_t ipart=0;ipart<kNSpecies;ipart++){
  88     Double_t miny=0.;
  89     Double_t maxy=1;
  90     TH2F *fHistoBayes=new TH2F(Form("BayesRec_%d",ipart),
  91                                Form("probability for reconstructed %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
  92     fHistoBayes->GetXaxis()->SetTitle("P_{T} (GeV / c)");
  93     fHistoBayes->GetYaxis()->SetTitle(Form("Bayes prob %s",kParticleSpeciesName[ipart]));
  94     fOutputList->Add(fHistoBayes);
  95   }
  96
  97   //nsigmaDC plot
  98   for(Int_t ipart=0;ipart<kNSpecies;ipart++){
  99     for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 100       Double_t miny=-10;
 101       Double_t maxy=10;
 102       if(ipid==kNSigmaTPCTOF){miny=0;maxy=20;}
 103       TH2F *fHistoNSigma=new TH2F(Form("NSigmaDC_%d_%d",ipart,ipid),
 104                                   Form("n#sigma for double counting %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
 105       fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
 106       fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
 107       fOutputList->Add(fHistoNSigma);
 108     }
 109   }
 110
 111   //nsigmaMC plot
 112   for(Int_t ipart=0;ipart<kNSpecies;ipart++){
 113     for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 114       Double_t miny=-30;
 115       Double_t maxy=30;
 116       if(ipid==kNSigmaTPCTOF){miny=0;maxy=50;}
 117       TH2F *fHistoNSigma=new TH2F(Form("NSigmaMC_%d_%d",ipart,ipid),
 118                                   Form("n#sigma for MC %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
 119       fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
 120       fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
 121       fOutputList->Add(fHistoNSigma);
 122     }
 123   }
 124
 125   //PID signal plot
 126   for(Int_t idet=0;idet<kNDetectors;idet++){
 127     for(Int_t ipart=0;ipart<kNSpecies;ipart++){
 128       Double_t maxy=500;
 129       if(idet==kTOF)maxy=1.1;
 130       TH2F *fHistoPID=new TH2F(Form("PID_%d_%d",idet,ipart),Form("%s signal - %s",kDetectorName[idet],kParticleSpeciesName[ipart]),200,0,10,500,-maxy,maxy);
 131       fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
 132       fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
 133       fOutputList->Add(fHistoPID);
 134     }
 135   }
 136   //PID signal plot, before PID cut
 137   for(Int_t idet=0;idet<kNDetectors;idet++){
 138     Double_t maxy=500;
 139     if(idet==kTOF)maxy=1.1;
 140     TH2F *fHistoPID=new TH2F(Form("PIDAll_%d",idet),Form("%s signal",kDetectorName[idet]),200,0,10,500,-maxy,maxy);
 141     fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
 142     fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
 143     fOutputList->Add(fHistoPID);
 144   }
 145 }
 146
 147 //////////////////////////////////////////////////////////////////////////////////////////////////
 148
 149 TH2F* AliHelperPID::GetHistogram2D(const char * name){
 150   // returns histo named name
 151   return (TH2F*) fOutputList->FindObject(name);
 152 }
 153
 154 //////////////////////////////////////////////////////////////////////////////////////////////////
 155
 156 Int_t AliHelperPID::GetParticleSpecies(AliVTrack * trk, Bool_t FIllQAHistos){
 157   //return the specie according to the minimum nsigma value
 158   //no double counting, this has to be evaluated using CheckDoubleCounting()
 159
 160   Int_t ID=kSpUndefined;
 161
 162   //get the PID response
 163   if(!fPIDResponse) {
 164     AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
 165     AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
 166     fPIDResponse = inputHandler->GetPIDResponse();
 167   }
 168   if(!fPIDResponse) {
 169     AliFatal("Cannot get pid response");
 170   }
 171
 172   //calculate nsigmas (used also by the bayesian)
 173   CalculateNSigmas(trk,FIllQAHistos);//fill the data member fnsigmas with the nsigmas value [ipart][iPID]
 174
 175   //Do PID
 176   if(fPIDType==kBayes){//use bayesianPID
 177
 178     if(!fPIDCombined) {
 179       // ------- setup PIDCombined
 180       fPIDCombined=new AliPIDCombined;
 181       fPIDCombined->SetDefaultTPCPriors();
 182       fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
 183     }
 184     if(!fPIDCombined) {
 185       AliFatal("PIDCombined object not found");
 186     }
 187
 188     ID = GetIDBayes(trk,FIllQAHistos);
 189
 190   }else{ //use nsigma PID
 191
 192     ID=FindMinNSigma(trk,FIllQAHistos);
 193
 194     if(fUseExclusiveNSigma){ //if one particle has double counting and exclusive nsigma is requested ID = kSpUndefined
 195       Bool_t *HasDC;
 196       HasDC=GetDoubleCounting(trk,FIllQAHistos);
 197       for(Int_t ipart=0;ipart<kNSpecies;ipart++){
 198         if(HasDC[ipart]==kTRUE)  ID = kSpUndefined;
 199       }
 200     }
 201   }
 202
 203   //Fill PID signal plot
 204   if(ID != kSpUndefined){
 205     for(Int_t idet=0;idet<kNDetectors;idet++){
 206       TH2F *h=GetHistogram2D(Form("PID_%d_%d",idet,ID));
 207       if(idet==kITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
 208       if(idet==kTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
 209       if(idet==kTOF && fHasTOFPID)h->Fill(trk->P(),TOFBetaCalc(trk)*trk->Charge());
 210     }
 211   }
 212   //Fill PID signal plot without cuts
 213   for(Int_t idet=0;idet<kNDetectors;idet++){
 214     TH2F *h=GetHistogram2D(Form("PIDAll_%d",idet));
 215     if(idet==kITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
 216     if(idet==kTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
 217     if(idet==kTOF && fHasTOFPID)h->Fill(trk->P(),TOFBetaCalc(trk)*trk->Charge());
 218   }
 219
 220   return ID;
 221 }
 222
 223 //////////////////////////////////////////////////////////////////////////////////////////////////
 224
 225 Int_t AliHelperPID::GetParticleSpecies(AliVParticle * part) {
 226   // return PID according to MC truth
 227   switch(TMath::Abs(part->PdgCode())){
 228   case 2212:
 229     return kSpProton;
 230     break;
 231   case 321:
 232     return kSpKaon;
 233     break;
 234   case 211:
 235     return kSpPion;
 236     break;
 237   default:
 238     return kSpUndefined;
 239   }
 240 }
 241
 242 //////////////////////////////////////////////////////////////////////////////////////////////////
 243
 244 Int_t AliHelperPID::GetIDBayes(AliVTrack * trk, Bool_t FIllQAHistos){
 245
 246   Bool_t *IDs=GetAllCompatibleIdentitiesNSigma(trk,FIllQAHistos);
 247
 248   Double_t probBayes[AliPID::kSPECIES];
 249
 250   UInt_t detUsed= 0;
 251   CheckTOF(trk);
 252   if(fHasTOFPID && trk->Pt()>fPtTOFPID){//use TOF information
 253     detUsed = CalcPIDCombined(trk, fPIDResponse, AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC, probBayes);
 254     if(detUsed != (AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC))return kSpUndefined;//check that TPC and TOF are used
 255   }else{
 256     detUsed = CalcPIDCombined(trk, fPIDResponse,AliPIDResponse::kDetTPC, probBayes);
 257     if(detUsed != AliPIDResponse::kDetTPC)return kSpUndefined;//check that TPC is used
 258   }
 259
 260   //the probability has to be normalized to one, we check it
 261   Double_t sump=0.;
 262   for(Int_t ipart=0;ipart<AliPID::kSPECIES;ipart++)sump+=probBayes[ipart];
 263   if(sump<.99 && sump>1.01){//FIXME precision problem in the sum, workaround
 264     AliFatal("Bayesian probability not normalized to one");
 265   }
 266
 267   //probabilities are normalized to one, if the cut is above .5 there is no problem
 268   if(probBayes[AliPID::kPion]>fBayesCut && IDs[kSpPion]==1){
 269     TH2F *h=GetHistogram2D(Form("BayesRec_%d",kSpPion));
 270     h->Fill(trk->Pt(),probBayes[AliPID::kPion]);
 271     return kSpPion;
 272   }
 273   else if(probBayes[AliPID::kKaon]>fBayesCut && IDs[kSpKaon]==1){
 274     TH2F *h=GetHistogram2D(Form("BayesRec_%d",kSpKaon));
 275     h->Fill(trk->Pt(),probBayes[AliPID::kKaon]);
 276     return kSpKaon;
 277   }
 278   else if(probBayes[AliPID::kProton]>fBayesCut && IDs[kSpProton]==1){
 279     TH2F *h=GetHistogram2D(Form("BayesRec_%d",kSpProton));
 280     h->Fill(trk->Pt(),probBayes[AliPID::kProton]);
 281     return kSpProton;
 282   }
 283   else{
 284     return kSpUndefined;
 285   }
 286 }
 287
 288 //////////////////////////////////////////////////////////////////////////////////////////////////
 289
 290 UInt_t AliHelperPID::CalcPIDCombined(const AliVTrack *track,const AliPIDResponse *PIDResponse, Int_t detMask, Double_t* prob) const{
 291   //
 292   // Bayesian PID calculation
 293   //
 294   for(Int_t i=0;i<AliPID::kSPECIES;i++)
 295     {
 296       prob[i]=0.;
 297     }
 298   fPIDCombined->SetDetectorMask(detMask);
 299
 300   return fPIDCombined->ComputeProbabilities((AliVTrack*)track, PIDResponse, prob);
 301 }
 302
 303 //////////////////////////////////////////////////////////////////////////////////////////////////
 304
 305 void AliHelperPID::CalculateNSigmas(AliVTrack * trk, Bool_t FIllQAHistos){
 306   //defines data member fnsigmas
 307
 308   // Compute nsigma for each hypthesis
 309   AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(trk);
 310   // --- TPC
 311   Double_t nsigmaTPCkProton = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton);
 312   Double_t nsigmaTPCkKaon   = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon);
 313   Double_t nsigmaTPCkPion   = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion);
 314   // --- TOF
 315   Double_t nsigmaTOFkProton=999.,nsigmaTOFkKaon=999.,nsigmaTOFkPion=999.;
 316   Double_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;
 317
 318   CheckTOF(trk);
 319
 320   if(fHasTOFPID && trk->Pt()>fPtTOFPID){//use TOF information
 321     nsigmaTOFkProton = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton);
 322     nsigmaTOFkKaon   = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon);
 323     nsigmaTOFkPion   = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion);
 324     Double_t d2Proton=nsigmaTPCkProton * nsigmaTPCkProton + nsigmaTOFkProton * nsigmaTOFkProton;
 325     Double_t d2Kaon=nsigmaTPCkKaon * nsigmaTPCkKaon + nsigmaTOFkKaon * nsigmaTOFkKaon;
 326     Double_t d2Pion=nsigmaTPCkPion * nsigmaTPCkPion + nsigmaTOFkPion * nsigmaTOFkPion;
 327     //commented, this is done in analogy with AliESDTrackCuts, nsigma combind according to the probability
 328     // --- combined
 329     // -----------------------------------
 330     // How to get to a n-sigma cut?
 331     //
 332     // The accumulated statistics from 0 to d is
 333     //
 334     // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 335     // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 336     // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
 337     //
 338     // work around precision problem
 339     // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 340     //if(TMath::Exp(- d2Proton / 2) > 1e-15)nsigmaTPCTOFkProton  =  TMath::Sqrt(2)*TMath::ErfInverse(1 - TMath::Exp(- d2Proton / 2));
 341     //if(TMath::Exp(- d2Kaon / 2) > 1e-15)nsigmaTPCTOFkKaon    = TMath::Sqrt(2)*TMath::ErfInverse(1 - TMath::Exp(- d2Kaon / 2));
 342     //if(TMath::Exp(- d2Pion / 2) > 1e-15)nsigmaTPCTOFkPion     = TMath::Sqrt(2)*TMath::ErfInverse(1 - TMath::Exp(- d2Pion / 2));
 343
 344     //used for the 2PC PID paper (circular cut)
 345     nsigmaTPCTOFkProton  =  TMath::Sqrt(d2Proton);
 346     nsigmaTPCTOFkKaon    =  TMath::Sqrt(d2Kaon);
 347     nsigmaTPCTOFkPion    =  TMath::Sqrt(d2Pion);
 348   }else{
 349     // --- combined
 350     // if TOF is missing and below fPtTOFPID only the TPC information is used
 351     nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
 352     nsigmaTPCTOFkKaon   = TMath::Abs(nsigmaTPCkKaon);
 353     nsigmaTPCTOFkPion   = TMath::Abs(nsigmaTPCkPion);
 354   }
 355
 356   //set data member fnsigmas
 357   fnsigmas[kSpPion][kNSigmaTPC]=nsigmaTPCkPion;
 358   fnsigmas[kSpKaon][kNSigmaTPC]=nsigmaTPCkKaon;
 359   fnsigmas[kSpProton][kNSigmaTPC]=nsigmaTPCkProton;
 360   fnsigmas[kSpPion][kNSigmaTOF]=nsigmaTOFkPion;
 361   fnsigmas[kSpKaon][kNSigmaTOF]=nsigmaTOFkKaon;
 362   fnsigmas[kSpProton][kNSigmaTOF]=nsigmaTOFkProton;
 363   fnsigmas[kSpPion][kNSigmaTPCTOF]=nsigmaTPCTOFkPion;
 364   fnsigmas[kSpKaon][kNSigmaTPCTOF]=nsigmaTPCTOFkKaon;
 365   fnsigmas[kSpProton][kNSigmaTPCTOF]=nsigmaTPCTOFkProton;
 366
 367   if(FIllQAHistos){
 368     //Fill NSigma SeparationPlot
 369     for(Int_t ipart=0;ipart<kNSpecies;ipart++){
 370       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 371         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && trk->Pt()<fPtTOFPID)continue;//not filling TOF and combined if no TOF PID
 372         TH2F *h=GetHistogram2D(Form("NSigma_%d_%d",ipart,ipid));
 373         h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
 374       }
 375     }
 376   }
 377 }
 378
 379 //////////////////////////////////////////////////////////////////////////////////////////////////
 380
 381 Int_t AliHelperPID::FindMinNSigma(AliVTrack * trk,Bool_t FillQAHistos){
 382
 383   CheckTOF(trk);
 384   if(fRequestTOFPID && (!fHasTOFPID) && trk->Pt()>fPtTOFPID)return kSpUndefined;
 385
 386   //get the identity of the particle with the minimum Nsigma
 387   Double_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
 388   switch (fPIDType){
 389   case kNSigmaTPC:
 390     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPC]);
 391     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPC])  ;
 392     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPC])  ;
 393     break;
 394   case kNSigmaTOF:
 395     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTOF]);
 396     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTOF])  ;
 397     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTOF])  ;
 398     break;
 399   case kNSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
 400     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
 401     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF])  ;
 402     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF])  ;
 403     break;
 404   case kBayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
 405     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
 406     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF])  ;
 407     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF])  ;
 408     break;
 409   }
 410
 411   // guess the particle based on the smaller nsigma (within fNSigmaPID)
 412   if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return kSpUndefined;//if is the default value for the three
 413
 414   if( ( nsigmaKaon   < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon   < fNSigmaPID)){
 415     if(FillQAHistos){
 416       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 417         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
 418         TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",kSpKaon,ipid));
 419         h->Fill(trk->Pt(),fnsigmas[kSpKaon][ipid]);
 420       }
 421     }
 422     return kSpKaon;
 423   }
 424   if( ( nsigmaPion   < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion   < fNSigmaPID)){
 425     if(FillQAHistos){
 426       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 427         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
 428         TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",kSpPion,ipid));
 429         h->Fill(trk->Pt(),fnsigmas[kSpPion][ipid]);
 430       }
 431     }
 432     return kSpPion;
 433   }
 434   if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)){
 435     if(FillQAHistos){
 436       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 437         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
 438         TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",kSpProton,ipid));
 439         h->Fill(trk->Pt(),fnsigmas[kSpProton][ipid]);
 440       }
 441     }
 442     return kSpProton;
 443   }
 444   // else, return undefined
 445   return kSpUndefined;
 446 }
 447
 448 //////////////////////////////////////////////////////////////////////////////////////////////////
 449
 450 Bool_t* AliHelperPID::GetDoubleCounting(AliVTrack * trk,Bool_t FIllQAHistos){
 451   //if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE
 452   //fill DC histos
 453   for(Int_t ipart=0;ipart<kNSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track
 454
 455   Int_t MinNSigma=FindMinNSigma(trk,kFALSE);//not filling the NSigmaRec histos
 456
 457   CheckTOF(trk);
 458
 459   if(MinNSigma==kSpUndefined)return fHasDoubleCounting;//in case of undefined no Double counting
 460
 461   Double_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
 462   switch (fPIDType) {
 463   case kNSigmaTPC:
 464     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPC]);
 465     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPC])  ;
 466     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPC])  ;
 467     break;
 468   case kNSigmaTOF:
 469     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTOF]);
 470     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTOF])  ;
 471     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTOF])  ;
 472     break;
 473   case kNSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
 474     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
 475     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF])  ;
 476     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF])  ;
 477     break;
 478   case kBayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
 479     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
 480     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF])  ;
 481     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF])  ;
 482     break;
 483   }
 484   if(nsigmaPion<fNSigmaPID && MinNSigma!=kSpPion)fHasDoubleCounting[kSpPion]=kTRUE;
 485   if(nsigmaKaon<fNSigmaPID && MinNSigma!=kSpKaon)fHasDoubleCounting[kSpKaon]=kTRUE;
 486   if(nsigmaProton<fNSigmaPID && MinNSigma!=kSpProton)fHasDoubleCounting[kSpProton]=kTRUE;
 487
 488   if(FIllQAHistos){
 489     //fill NSigma distr for double counting
 490     for(Int_t ipart=0;ipart<kNSpecies;ipart++){
 491       if(fHasDoubleCounting[ipart]){
 492         for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 493           if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
 494           TH2F *h=GetHistogram2D(Form("NSigmaDC_%d_%d",ipart,ipid));
 495           h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
 496         }
 497       }
 498     }
 499   }
 500   return fHasDoubleCounting;
 501 }
 502
 503 //////////////////////////////////////////////////////////////////////////////////////////////////
 504
 505 Bool_t* AliHelperPID::GetAllCompatibleIdentitiesNSigma(AliVTrack * trk,Bool_t FIllQAHistos){
 506
 507   Bool_t *IDs=GetDoubleCounting(trk,FIllQAHistos);
 508   IDs[FindMinNSigma(trk,FIllQAHistos)]=kTRUE;
 509   return IDs;
 510
 511 }
 512
 513 //////////////////////////////////////////////////////////////////////////////////////////////////
 514
 515 Int_t AliHelperPID::GetMCParticleSpecie(AliVEvent* event, AliVTrack * trk, Bool_t FillQAHistos){
 516   //return the specie according to the MC truth
 517   CheckTOF(trk);
 518
 519   if(!fisMC)AliFatal("Error: AliHelperPID::GetMCParticleSpecie called on data\n");
 520
 521   Int_t code=999;
 522   Bool_t isAOD=kFALSE;
 523   Bool_t isESD=kFALSE;
 524   TString nameoftrack(event->ClassName());
 525   if(!nameoftrack.CompareTo("AliESDEvent"))isESD=kTRUE;
 526   else if(!nameoftrack.CompareTo("AliAODEvent"))isAOD=kTRUE;
 527   else AliFatal("Not processing AODs nor ESDs") ;
 528
 529   if(isAOD){
 530     TClonesArray *arrayMC = 0;
 531     arrayMC = (TClonesArray*) event->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 532     if (!arrayMC)AliFatal("Error: MC particles branch not found!\n");
 533     AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(trk->GetLabel()));
 534     if (!partMC){
 535       AliError("Cannot get MC particle");
 536       return kSpUndefined;
 537     }
 538     code=partMC->GetPdgCode();
 539   }
 540   if(isESD){
 541     AliStack* stack=0x0;
 542     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 543     if (!eventHandler) AliFatal("ERROR: Could not retrieve MC event handler");
 544     AliMCEvent* mcEvent = eventHandler->MCEvent();
 545     if (!mcEvent)AliFatal("ERROR: Could not retrieve MC event");
 546     stack = mcEvent->Stack();
 547     if (!stack) AliFatal("ERROR: stack not available\n");
 548     TParticle *part=0;
 549     part = (TParticle*)stack->Particle(TMath::Abs(trk->GetLabel()));
 550     if (!part){
 551       AliError("Cannot get MC particle");
 552       return kSpUndefined;
 553     }
 554     code = part->GetPdgCode();
 555   }
 556   switch(TMath::Abs(code)){
 557   case 2212:
 558     if(FillQAHistos){
 559       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 560         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
 561         TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",kSpProton,ipid));
 562         h->Fill(trk->Pt(),fnsigmas[kSpProton][ipid]);
 563       }
 564     }
 565     return kSpProton;
 566     break;
 567   case 321:
 568     if(FillQAHistos){
 569       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 570         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
 571         TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",kSpKaon,ipid));
 572         h->Fill(trk->Pt(),fnsigmas[kSpKaon][ipid]);
 573       }
 574     }
 575     return kSpKaon;
 576     break;
 577   case 211:
 578     if(FillQAHistos){
 579       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 580         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
 581         TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",kSpPion,ipid));
 582         h->Fill(trk->Pt(),fnsigmas[kSpPion][ipid]);
 583       }
 584     }
 585     return kSpPion;
 586     break;
 587   default:
 588     return kSpUndefined;
 589   }
 590
 591 }
 592
 593 //////////////////////////////////////////////////////////////////////////////////////////////////
 594
 595 void AliHelperPID::CheckTOF(AliVTrack * trk)
 596 {
 597   //check if the particle has TOF Matching
 598
 599   //get the PIDResponse
 600   if(fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF,trk)==0)fHasTOFPID=kFALSE;
 601   else fHasTOFPID=kTRUE;
 602
 603   //in addition to TOF status we look at the pt
 604   if(trk->Pt()<fPtTOFPID)fHasTOFPID=kFALSE;
 605
 606   if(fRemoveTracksT0Fill)
 607     {
 608       Int_t startTimeMask = fPIDResponse->GetTOFResponse().GetStartTimeMask(trk->P());
 609       if (startTimeMask < 0)fHasTOFPID=kFALSE;
 610     }
 611 }
 612
 613 //////////////////////////////////////////////////////////////////////////////////////////////////
 614
 615 Double_t AliHelperPID::TOFBetaCalc(AliVTrack *track) const{
 616   //TOF beta calculation
 617   Double_t tofTime=track->GetTOFsignal();
 618
 619   Double_t c=TMath::C()*1.E-9;// m/ns
 620   Float_t startTime = fPIDResponse->GetTOFResponse().GetStartTime(((AliVTrack*)track)->P());//in ps
 621   Double_t length= fPIDResponse->GetTOFResponse().GetExpectedSignal(track,AliPID::kElectron)*1E-3*c;
 622   tofTime -= startTime;      // subtract startTime to the signal
 623   Double_t tof= tofTime*1E-3; // ns, average T0 fill subtracted, no info from T0detector
 624   tof=tof*c;
 625   return length/tof;
 626 }
 627
 628 //////////////////////////////////////////////////////////////////////////////////////////////////
 629
 630 Double_t AliHelperPID::GetMass(AliHelperParticleSpecies_t id) const{
 631   //return Mass according to AliHelperParticleSpecies_t. If undefined return -999.
 632   Double_t mass=-999.;
 633   if (id == kSpProton) { mass=9.38271999999999995e-01; }
 634   if (id == kSpKaon)   { mass=4.93676999999999977e-01; }
 635   if (id == kSpPion)    { mass=1.39570000000000000e-01; }
 636   return mass;
 637 }
 638
 639 //////////////////////////////////////////////////////////////////////////////////////////////////
 640
 641 Long64_t AliHelperPID::Merge(TCollection* list)
 642 {
 643   // Merging interface.
 644   // Returns the number of merged objects (including this).
 645   Printf("Merging");
 646   if (!list)
 647     return 0;
 648   if (list->IsEmpty())
 649     return 1;
 650   TIterator* iter = list->MakeIterator();
 651   TObject* obj;
 652   // collections of TList with output histos
 653   TList collections;
 654   Int_t count = 0;
 655   while ((obj = iter->Next())) {
 656     AliHelperPID* entry = dynamic_cast<AliHelperPID*> (obj);
 657     if (entry == 0)
 658       continue;
 659     TList * outputlist = entry->GetOutputList();
 660     collections.Add(outputlist);
 661     count++;
 662   }
 663   fOutputList->Merge(&collections);
 664   delete iter;
 665   Printf("OK");
 666   return count+1;
 667 }