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 "AliAnalysisManager.h"
  35 #include "AliInputEventHandler.h"
  36
  37 using namespace AliHelperPIDNameSpace;
  38 using namespace std;
  39
  40 const char * kPIDTypeName[]={"TPC","TOF","TPC-TOF"} ;
  41 const char * kDetectorName[]={"TPC","TOF"} ;
  42 const char * kParticleSpeciesName[]={"Pions","Kaons","Protons","Undefined"} ;
  43
  44 ClassImp(AliHelperPID)
  45
  46 AliHelperPID::AliHelperPID() : TNamed("HelperPID", "PID object"),fisMC(0), fPIDType(kNSigmaTPCTOF), fNSigmaPID(3), fPIDResponse(0),fOutputList(0),fRequestTOFPID(1),fRemoveTracksT0Fill(0),fUseExclusiveNSigma(0),fPtTOFPID(.6),fHasTOFPID(0){
  47
  48   for(Int_t ipart=0;ipart<kNSpecies;ipart++)
  49     for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++)
  50       fnsigmas[ipart][ipid]=999.;
  51
  52   for(Int_t ipart=0;ipart<kNSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;
  53
  54   fOutputList = new TList;
  55   fOutputList->SetOwner();
  56   fOutputList->SetName("OutputList");
  57
  58   //nsigma plot
  59   for(Int_t ipart=0;ipart<kNSpecies;ipart++){
  60     for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
  61       Double_t miny=-30;
  62       Double_t maxy=30;
  63       if(ipid==kNSigmaTPCTOF){miny=0;maxy=50;}
  64       TH2F *fHistoNSigma=new TH2F(Form("NSigma_%d_%d",ipart,ipid),Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,200,miny,maxy);
  65       fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
  66       fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
  67       fOutputList->Add(fHistoNSigma);
  68     }
  69   }
  70
  71   //nsigmaRec plot
  72   for(Int_t ipart=0;ipart<kNSpecies;ipart++){
  73     for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
  74       Double_t miny=-10;
  75       Double_t maxy=10;
  76       if(ipid==kNSigmaTPCTOF){miny=0;maxy=20;}
  77       TH2F *fHistoNSigma=new TH2F(Form("NSigmaRec_%d_%d",ipart,ipid),
  78                                   Form("n#sigma for reconstructed %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,100,miny,maxy);
  79       fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
  80       fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
  81       fOutputList->Add(fHistoNSigma);
  82     }
  83   }
  84
  85   //nsigmaDC plot
  86   for(Int_t ipart=0;ipart<kNSpecies;ipart++){
  87     for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
  88       Double_t miny=-10;
  89       Double_t maxy=10;
  90       if(ipid==kNSigmaTPCTOF){miny=0;maxy=20;}
  91       TH2F *fHistoNSigma=new TH2F(Form("NSigmaDC_%d_%d",ipart,ipid),
  92                                   Form("n#sigma for double counting %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,100,miny,maxy);
  93       fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
  94       fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
  95       fOutputList->Add(fHistoNSigma);
  96     }
  97   }
  98
  99   //nsigmaMC plot
 100   for(Int_t ipart=0;ipart<kNSpecies;ipart++){
 101     for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 102       Double_t miny=-30;
 103       Double_t maxy=30;
 104       if(ipid==kNSigmaTPCTOF){miny=0;maxy=50;}
 105       TH2F *fHistoNSigma=new TH2F(Form("NSigmaMC_%d_%d",ipart,ipid),
 106                                   Form("n#sigma for MC %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,100,miny,maxy);
 107       fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
 108       fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
 109       fOutputList->Add(fHistoNSigma);
 110     }
 111   }
 112
 113   //PID signal plot
 114   for(Int_t idet=0;idet<kNDetectors;idet++){
 115     Double_t maxy=1000;
 116     if(idet==kTOF)maxy=200000;
 117     TH2F *fHistoPID=new TH2F(Form("PID_%d",idet),Form("%s signal",kDetectorName[idet]),200,0,10,200,-maxy,maxy);
 118     fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
 119     fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
 120     fOutputList->Add(fHistoPID);
 121   }
 122 }
 123
 124 //////////////////////////////////////////////////////////////////////////////////////////////////
 125
 126 TH2F* AliHelperPID::GetHistogram2D(const char * name){
 127   // returns histo named name
 128   return (TH2F*) fOutputList->FindObject(name);
 129 }
 130
 131 //////////////////////////////////////////////////////////////////////////////////////////////////
 132
 133 Int_t AliHelperPID::GetParticleSpecies(AliVTrack * trk, Bool_t FIllQAHistos){
 134   //return the specie according to the minimum nsigma value
 135   //no double counting, this has to be evaluated using CheckDoubleCounting()
 136   //Printf("fPtTOFPID %.1f, fRequestTOFPID %d, fNSigmaPID %.1f, fPIDType %d",fPtTOFPID,fRequestTOFPID,fNSigmaPID,fPIDType);
 137
 138   CalculateNSigmas(trk,FIllQAHistos);//fill the data member fnsigmas with the nsigmas value [ipart][iPID]
 139
 140   if(fUseExclusiveNSigma){
 141     Bool_t *HasDC;
 142     HasDC=GetDoubleCounting(trk,kFALSE);
 143     for(Int_t ipart=0;ipart<kNSpecies;ipart++){
 144       if(HasDC[ipart]==kTRUE)  return kSpUndefined;
 145     }
 146     return FindMinNSigma(trk,FIllQAHistos);//NSigmaRec distr filled here
 147   }
 148   else return FindMinNSigma(trk,FIllQAHistos);//NSigmaRec distr filled here
 149
 150 }
 151
 152 //////////////////////////////////////////////////////////////////////////////////////////////////
 153
 154 void AliHelperPID::CalculateNSigmas(AliVTrack * trk, Bool_t FIllQAHistos){
 155   //defines data member fnsigmas
 156   if(!fPIDResponse) {
 157     AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
 158     AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
 159     fPIDResponse = inputHandler->GetPIDResponse();
 160   }
 161   if(!fPIDResponse) {
 162     AliFatal("Cannot get pid response");
 163   }
 164
 165   // Compute nsigma for each hypthesis
 166   AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(trk);
 167   // --- TPC
 168   Double_t nsigmaTPCkProton = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton);
 169   Double_t nsigmaTPCkKaon   = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon);
 170   Double_t nsigmaTPCkPion   = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion);
 171   // --- TOF
 172   Double_t nsigmaTOFkProton=999.,nsigmaTOFkKaon=999.,nsigmaTOFkPion=999.;
 173   Double_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;
 174
 175   CheckTOF(trk);
 176
 177   if(fHasTOFPID && trk->Pt()>fPtTOFPID){//use TOF information
 178     nsigmaTOFkProton = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton);
 179     nsigmaTOFkKaon   = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon);
 180     nsigmaTOFkPion   = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion);
 181     // --- combined
 182     nsigmaTPCTOFkProton = TMath::Sqrt((nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton)/2.);
 183     nsigmaTPCTOFkKaon   = TMath::Sqrt((nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon)/2.);
 184     nsigmaTPCTOFkPion   = TMath::Sqrt((nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion)/2.);
 185   }else{
 186     // --- combined
 187     // if TOF is missing and below fPtTOFPID only the TPC information is used
 188     nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
 189     nsigmaTPCTOFkKaon   = TMath::Abs(nsigmaTPCkKaon);
 190     nsigmaTPCTOFkPion   = TMath::Abs(nsigmaTPCkPion);
 191   }
 192
 193   //set data member fnsigmas
 194   fnsigmas[kSpPion][kNSigmaTPC]=nsigmaTPCkPion;
 195   fnsigmas[kSpKaon][kNSigmaTPC]=nsigmaTPCkKaon;
 196   fnsigmas[kSpProton][kNSigmaTPC]=nsigmaTPCkProton;
 197   fnsigmas[kSpPion][kNSigmaTOF]=nsigmaTOFkPion;
 198   fnsigmas[kSpKaon][kNSigmaTOF]=nsigmaTOFkKaon;
 199   fnsigmas[kSpProton][kNSigmaTOF]=nsigmaTOFkProton;
 200   fnsigmas[kSpPion][kNSigmaTPCTOF]=nsigmaTPCTOFkPion;
 201   fnsigmas[kSpKaon][kNSigmaTPCTOF]=nsigmaTPCTOFkKaon;
 202   fnsigmas[kSpProton][kNSigmaTPCTOF]=nsigmaTPCTOFkProton;
 203
 204   if(FIllQAHistos){
 205     //Fill NSigma SeparationPlot
 206     for(Int_t ipart=0;ipart<kNSpecies;ipart++){
 207       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 208         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && trk->Pt()<fPtTOFPID)continue;//not filling TOF and combined if no TOF PID
 209         TH2F *h=GetHistogram2D(Form("NSigma_%d_%d",ipart,ipid));
 210         h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
 211       }
 212     }
 213     //Fill PID signal plot
 214     for(Int_t idet=0;idet<kNDetectors;idet++){
 215       TH2F *h=GetHistogram2D(Form("PID_%d",idet));
 216       if(idet==kTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
 217       if(idet==kTOF && fHasTOFPID)h->Fill(trk->P(),trk->GetTOFsignal()*trk->Charge());
 218     }
 219   }
 220 }
 221
 222 //////////////////////////////////////////////////////////////////////////////////////////////////
 223
 224 Int_t AliHelperPID::FindMinNSigma(AliVTrack * trk,Bool_t FillQAHistos){
 225
 226   CheckTOF(trk);
 227   if(fRequestTOFPID && (!fHasTOFPID) && trk->Pt()>fPtTOFPID)return kSpUndefined;
 228
 229   //get the identity of the particle with the minimum Nsigma
 230   Double_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
 231   switch (fPIDType){
 232   case kNSigmaTPC:
 233     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPC]);
 234     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPC])  ;
 235     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPC])  ;
 236     break;
 237   case kNSigmaTOF:
 238     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTOF]);
 239     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTOF])  ;
 240     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTOF])  ;
 241     break;
 242   case kNSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
 243     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
 244     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF])  ;
 245     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF])  ;
 246     break;
 247   }
 248
 249   // guess the particle based on the smaller nsigma (within fNSigmaPID)
 250   if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return kSpUndefined;//if is the default value for the three
 251
 252   if( ( nsigmaKaon   < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon   < fNSigmaPID)){
 253     if(FillQAHistos){
 254       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 255         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
 256         TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",kSpKaon,ipid));
 257         h->Fill(trk->Pt(),fnsigmas[kSpKaon][ipid]);
 258       }
 259     }
 260     return kSpKaon;
 261   }
 262   if( ( nsigmaPion   < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion   < fNSigmaPID)){
 263     if(FillQAHistos){
 264       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 265         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
 266         TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",kSpPion,ipid));
 267         h->Fill(trk->Pt(),fnsigmas[kSpPion][ipid]);
 268       }
 269     }
 270     return kSpPion;
 271   }
 272   if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)){
 273     if(FillQAHistos){
 274       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 275         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
 276         TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",kSpProton,ipid));
 277         h->Fill(trk->Pt(),fnsigmas[kSpProton][ipid]);
 278       }
 279     }
 280     return kSpProton;
 281   }
 282   // else, return undefined
 283   return kSpUndefined;
 284 }
 285
 286 //////////////////////////////////////////////////////////////////////////////////////////////////
 287
 288 Bool_t* AliHelperPID::GetDoubleCounting(AliVTrack * trk,Bool_t FIllQAHistos){
 289   //if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE
 290   //fill DC histos
 291   for(Int_t ipart=0;ipart<=kNSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track
 292
 293   Int_t MinNSigma=FindMinNSigma(trk,kFALSE);//not filling the NSigmaRec histos
 294
 295   CheckTOF(trk);
 296
 297   if(MinNSigma==kSpUndefined)return fHasDoubleCounting;//in case of undefined no Double counting
 298
 299   Double_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
 300   switch (fPIDType) {
 301   case kNSigmaTPC:
 302     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPC]);
 303     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPC])  ;
 304     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPC])  ;
 305     break;
 306   case kNSigmaTOF:
 307     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTOF]);
 308     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTOF])  ;
 309     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTOF])  ;
 310     break;
 311   case kNSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
 312     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
 313     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF])  ;
 314     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF])  ;
 315     break;
 316   }
 317   if(nsigmaPion<fNSigmaPID && MinNSigma!=kSpPion)fHasDoubleCounting[kSpPion]=kTRUE;
 318   if(nsigmaKaon<fNSigmaPID && MinNSigma!=kSpKaon)fHasDoubleCounting[kSpKaon]=kTRUE;
 319   if(nsigmaProton<fNSigmaPID && MinNSigma!=kSpProton)fHasDoubleCounting[kSpProton]=kTRUE;
 320
 321   if(FIllQAHistos){
 322     //fill NSigma distr for double counting
 323     for(Int_t ipart=0;ipart<kNSpecies;ipart++){
 324       if(fHasDoubleCounting[ipart]){
 325         for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 326           if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
 327           TH2F *h=GetHistogram2D(Form("NSigmaDC_%d_%d",ipart,ipid));
 328           h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
 329         }
 330       }
 331     }
 332   }
 333   return fHasDoubleCounting;
 334 }
 335
 336 //////////////////////////////////////////////////////////////////////////////////////////////////
 337
 338 Int_t AliHelperPID::GetMCParticleSpecie(AliVEvent* event, AliVTrack * trk, Bool_t FillQAHistos){
 339   //return the specie according to the MC truth
 340   CheckTOF(trk);
 341
 342   if(!fisMC)AliFatal("Error: AliHelperPID::GetMCParticleSpecie called on data\n");
 343
 344   Int_t code=999;
 345   Bool_t isAOD=kFALSE;
 346   Bool_t isESD=kFALSE;
 347   TString nameoftrack(event->ClassName());
 348   if(!nameoftrack.CompareTo("AliESDEvent"))isESD=kTRUE;
 349   else if(!nameoftrack.CompareTo("AliAODEvent"))isAOD=kTRUE;
 350   else AliFatal("Not processing AODs nor ESDs") ;
 351
 352   if(isAOD){
 353     TClonesArray *arrayMC = 0;
 354     arrayMC = (TClonesArray*) event->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 355     if (!arrayMC)AliFatal("Error: MC particles branch not found!\n");
 356     AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(trk->GetLabel()));
 357     if (!partMC){
 358       AliError("Cannot get MC particle");
 359       return kSpUndefined;
 360     }
 361     code=partMC->GetPdgCode();
 362   }
 363   if(isESD){
 364     AliStack* stack=0x0;
 365     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 366     if (!eventHandler) AliFatal("ERROR: Could not retrieve MC event handler");
 367     AliMCEvent* mcEvent = eventHandler->MCEvent();
 368     if (!mcEvent)AliFatal("ERROR: Could not retrieve MC event");
 369     stack = mcEvent->Stack();
 370     if (!stack) AliFatal("ERROR: stack not available\n");
 371     TParticle *part=0;
 372     part = (TParticle*)stack->Particle(TMath::Abs(trk->GetLabel()));
 373     if (!part){
 374       AliError("Cannot get MC particle");
 375       return kSpUndefined;
 376     }
 377     code = part->GetPdgCode();
 378   }
 379   switch(TMath::Abs(code)){
 380   case 2212:
 381     if(FillQAHistos){
 382       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 383         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
 384         TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",kSpProton,ipid));
 385         h->Fill(trk->Pt(),fnsigmas[kSpProton][ipid]);
 386       }
 387     }
 388     return kSpProton;
 389     break;
 390   case 321:
 391     if(FillQAHistos){
 392       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 393         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
 394         TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",kSpKaon,ipid));
 395         h->Fill(trk->Pt(),fnsigmas[kSpKaon][ipid]);
 396       }
 397     }
 398     return kSpKaon;
 399     break;
 400   case 211:
 401     if(FillQAHistos){
 402       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 403         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
 404         TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",kSpPion,ipid));
 405         h->Fill(trk->Pt(),fnsigmas[kSpPion][ipid]);
 406       }
 407     }
 408     return kSpPion;
 409     break;
 410   default:
 411     return kSpUndefined;
 412   }
 413
 414 }
 415
 416 //////////////////////////////////////////////////////////////////////////////////////////////////
 417
 418 void AliHelperPID::CheckTOF(AliVTrack * trk)
 419 {
 420   //check if the particle has TOF Matching
 421   UInt_t status;
 422   status=trk->GetStatus();
 423   if((status&AliVTrack::kTOFout)==0 || (status&AliVTrack::kTIME)==0)fHasTOFPID=kFALSE;
 424   else fHasTOFPID=kTRUE;
 425   //in addition to KTOFout and kTIME we look at the pt
 426   if(trk->Pt()<fPtTOFPID)fHasTOFPID=kFALSE;
 427
 428   if(fRemoveTracksT0Fill)
 429     {
 430       //get the PIDResponse
 431       if(!fPIDResponse) {
 432         AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
 433         AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
 434         fPIDResponse = inputHandler->GetPIDResponse();
 435       }
 436       if(!fPIDResponse) {
 437         AliFatal("Cannot get pid response");
 438       }
 439       Int_t startTimeMask = fPIDResponse->GetTOFResponse().GetStartTimeMask(trk->P());
 440       if (startTimeMask < 0)fHasTOFPID=kFALSE;
 441     }
 442 }
 443
 444 //////////////////////////////////////////////////////////////////////////////////////////////////
 445
 446 Long64_t AliHelperPID::Merge(TCollection* list)
 447 {
 448   // Merging interface.
 449   // Returns the number of merged objects (including this).
 450   Printf("Merging");
 451   if (!list)
 452     return 0;
 453   if (list->IsEmpty())
 454     return 1;
 455   TIterator* iter = list->MakeIterator();
 456   TObject* obj;
 457   // collections of TList with output histos
 458   TList collections;
 459   Int_t count = 0;
 460   while ((obj = iter->Next())) {
 461     AliHelperPID* entry = dynamic_cast<AliHelperPID*> (obj);
 462     if (entry == 0)
 463       continue;
 464     TList * outputlist = entry->GetOutputList();
 465     collections.Add(outputlist);
 466     count++;
 467   }
 468   fOutputList->Merge(&collections);
 469   delete iter;
 470   Printf("OK");
 471   return count+1;
 472 }