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),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   }
 147   else return FindMinNSigma(trk,FIllQAHistos);//NSigmaRec distr filled here
 148
 149 }
 150
 151 //////////////////////////////////////////////////////////////////////////////////////////////////
 152
 153 void AliHelperPID::CalculateNSigmas(AliVTrack * trk, Bool_t FIllQAHistos){
 154   //defines data member fnsigmas
 155   if(!fPIDResponse) {
 156     AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
 157     AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
 158     fPIDResponse = inputHandler->GetPIDResponse();
 159   }
 160   if(!fPIDResponse) {
 161     AliFatal("Cannot get pid response");
 162   }
 163
 164   // Compute nsigma for each hypthesis
 165   AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(trk);
 166   // --- TPC
 167   Double_t nsigmaTPCkProton = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton);
 168   Double_t nsigmaTPCkKaon   = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon);
 169   Double_t nsigmaTPCkPion   = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion);
 170   // --- TOF
 171   Double_t nsigmaTOFkProton=999.,nsigmaTOFkKaon=999.,nsigmaTOFkPion=999.;
 172   Double_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;
 173
 174   CheckTOF(trk);
 175
 176   if(fHasTOFPID && trk->Pt()>fPtTOFPID){//use TOF information
 177     nsigmaTOFkProton = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton);
 178     nsigmaTOFkKaon   = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon);
 179     nsigmaTOFkPion   = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion);
 180     // --- combined
 181     nsigmaTPCTOFkProton = TMath::Sqrt((nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton)/2.);
 182     nsigmaTPCTOFkKaon   = TMath::Sqrt((nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon)/2.);
 183     nsigmaTPCTOFkPion   = TMath::Sqrt((nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion)/2.);
 184   }else{
 185     // --- combined
 186     // if TOF is missing and below fPtTOFPID only the TPC information is used
 187     nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
 188     nsigmaTPCTOFkKaon   = TMath::Abs(nsigmaTPCkKaon);
 189     nsigmaTPCTOFkPion   = TMath::Abs(nsigmaTPCkPion);
 190   }
 191
 192   //set data member fnsigmas
 193   fnsigmas[kSpPion][kNSigmaTPC]=nsigmaTPCkPion;
 194   fnsigmas[kSpKaon][kNSigmaTPC]=nsigmaTPCkKaon;
 195   fnsigmas[kSpProton][kNSigmaTPC]=nsigmaTPCkProton;
 196   fnsigmas[kSpPion][kNSigmaTOF]=nsigmaTOFkPion;
 197   fnsigmas[kSpKaon][kNSigmaTOF]=nsigmaTOFkKaon;
 198   fnsigmas[kSpProton][kNSigmaTOF]=nsigmaTOFkProton;
 199   fnsigmas[kSpPion][kNSigmaTPCTOF]=nsigmaTPCTOFkPion;
 200   fnsigmas[kSpKaon][kNSigmaTPCTOF]=nsigmaTPCTOFkKaon;
 201   fnsigmas[kSpProton][kNSigmaTPCTOF]=nsigmaTPCTOFkProton;
 202
 203   if(FIllQAHistos){
 204     //Fill NSigma SeparationPlot
 205     for(Int_t ipart=0;ipart<kNSpecies;ipart++){
 206       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 207         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && trk->Pt()<fPtTOFPID)continue;//not filling TOF and combined if no TOF PID
 208         TH2F *h=GetHistogram2D(Form("NSigma_%d_%d",ipart,ipid));
 209         h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
 210       }
 211     }
 212     //Fill PID signal plot
 213     for(Int_t idet=0;idet<kNDetectors;idet++){
 214       TH2F *h=GetHistogram2D(Form("PID_%d",idet));
 215       if(idet==kTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
 216       if(idet==kTOF && fHasTOFPID)h->Fill(trk->P(),trk->GetTOFsignal()*trk->Charge());
 217     }
 218   }
 219 }
 220
 221 //////////////////////////////////////////////////////////////////////////////////////////////////
 222
 223 Int_t AliHelperPID::FindMinNSigma(AliVTrack * trk,Bool_t FillQAHistos){
 224
 225   CheckTOF(trk);
 226   if(fRequestTOFPID && (!fHasTOFPID) && trk->Pt()>fPtTOFPID)return kSpUndefined;
 227
 228   //get the identity of the particle with the minimum Nsigma
 229   Double_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
 230   switch (fPIDType){
 231   case kNSigmaTPC:
 232     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPC]);
 233     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPC])  ;
 234     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPC])  ;
 235     break;
 236   case kNSigmaTOF:
 237     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTOF]);
 238     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTOF])  ;
 239     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTOF])  ;
 240     break;
 241   case kNSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
 242     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
 243     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF])  ;
 244     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF])  ;
 245     break;
 246   }
 247
 248   // guess the particle based on the smaller nsigma (within fNSigmaPID)
 249   if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return kSpUndefined;//if is the default value for the three
 250
 251   if( ( nsigmaKaon   < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon   < fNSigmaPID)){
 252     if(FillQAHistos){
 253       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 254         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
 255         TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",kSpKaon,ipid));
 256         h->Fill(trk->Pt(),fnsigmas[kSpKaon][ipid]);
 257       }
 258     }
 259     return kSpKaon;
 260   }
 261   if( ( nsigmaPion   < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion   < fNSigmaPID)){
 262     if(FillQAHistos){
 263       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 264         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
 265         TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",kSpPion,ipid));
 266         h->Fill(trk->Pt(),fnsigmas[kSpPion][ipid]);
 267       }
 268     }
 269     return kSpPion;
 270   }
 271   if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)){
 272     if(FillQAHistos){
 273       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 274         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
 275         TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",kSpProton,ipid));
 276         h->Fill(trk->Pt(),fnsigmas[kSpProton][ipid]);
 277       }
 278     }
 279     return kSpProton;
 280   }
 281   // else, return undefined
 282   return kSpUndefined;
 283 }
 284
 285 //////////////////////////////////////////////////////////////////////////////////////////////////
 286
 287 Bool_t* AliHelperPID::GetDoubleCounting(AliVTrack * trk,Bool_t FIllQAHistos){
 288   //if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE
 289   //fill DC histos
 290   for(Int_t ipart=0;ipart<=kNSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track
 291
 292   Int_t MinNSigma=FindMinNSigma(trk,kFALSE);//not filling the NSigmaRec histos
 293
 294   CheckTOF(trk);
 295
 296   if(MinNSigma==kSpUndefined)return fHasDoubleCounting;//in case of undefined no Double counting
 297
 298   Double_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
 299   switch (fPIDType) {
 300   case kNSigmaTPC:
 301     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPC]);
 302     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPC])  ;
 303     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPC])  ;
 304     break;
 305   case kNSigmaTOF:
 306     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTOF]);
 307     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTOF])  ;
 308     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTOF])  ;
 309     break;
 310   case kNSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
 311     nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
 312     nsigmaKaon    =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF])  ;
 313     nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF])  ;
 314     break;
 315   }
 316   if(nsigmaPion<fNSigmaPID && MinNSigma!=kSpPion)fHasDoubleCounting[kSpPion]=kTRUE;
 317   if(nsigmaKaon<fNSigmaPID && MinNSigma!=kSpKaon)fHasDoubleCounting[kSpKaon]=kTRUE;
 318   if(nsigmaProton<fNSigmaPID && MinNSigma!=kSpProton)fHasDoubleCounting[kSpProton]=kTRUE;
 319
 320   if(FIllQAHistos){
 321     //fill NSigma distr for double counting
 322     for(Int_t ipart=0;ipart<kNSpecies;ipart++){
 323       if(fHasDoubleCounting[ipart]){
 324         for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 325           if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
 326           TH2F *h=GetHistogram2D(Form("NSigmaDC_%d_%d",ipart,ipid));
 327           h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
 328         }
 329       }
 330     }
 331   }
 332   return fHasDoubleCounting;
 333 }
 334
 335 //////////////////////////////////////////////////////////////////////////////////////////////////
 336
 337 Int_t AliHelperPID::GetMCParticleSpecie(AliVEvent* event, AliVTrack * trk, Bool_t FillQAHistos){
 338   //return the specie according to the MC truth
 339   CheckTOF(trk);
 340
 341   if(!fisMC)AliFatal("Error: AliHelperPID::GetMCParticleSpecie called on data\n");
 342
 343   Int_t code=999;
 344   Bool_t isAOD=kFALSE;
 345   Bool_t isESD=kFALSE;
 346   TString nameoftrack(event->ClassName());
 347   if(!nameoftrack.CompareTo("AliESDEvent"))isESD=kTRUE;
 348   else if(!nameoftrack.CompareTo("AliAODEvent"))isAOD=kTRUE;
 349   else AliFatal("Not processing AODs nor ESDs") ;
 350
 351   if(isAOD){
 352     TClonesArray *arrayMC = 0;
 353     arrayMC = (TClonesArray*) event->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 354     if (!arrayMC)AliFatal("Error: MC particles branch not found!\n");
 355     AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(trk->GetLabel()));
 356     if (!partMC){
 357       AliError("Cannot get MC particle");
 358       return kSpUndefined;
 359     }
 360     code=partMC->GetPdgCode();
 361   }
 362   if(isESD){
 363     AliStack* stack=0x0;
 364     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 365     if (!eventHandler) AliFatal("ERROR: Could not retrieve MC event handler");
 366     AliMCEvent* mcEvent = eventHandler->MCEvent();
 367     if (!mcEvent)AliFatal("ERROR: Could not retrieve MC event");
 368     stack = mcEvent->Stack();
 369     if (!stack) AliFatal("ERROR: stack not available\n");
 370     TParticle *part=0;
 371     part = (TParticle*)stack->Particle(TMath::Abs(trk->GetLabel()));
 372     if (!part){
 373       AliError("Cannot get MC particle");
 374       return kSpUndefined;
 375     }
 376     code = part->GetPdgCode();
 377   }
 378   switch(TMath::Abs(code)){
 379   case 2212:
 380     if(FillQAHistos){
 381       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 382         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
 383         TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",kSpProton,ipid));
 384         h->Fill(trk->Pt(),fnsigmas[kSpProton][ipid]);
 385       }
 386     }
 387     return kSpProton;
 388     break;
 389   case 321:
 390     if(FillQAHistos){
 391       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 392         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
 393         TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",kSpKaon,ipid));
 394         h->Fill(trk->Pt(),fnsigmas[kSpKaon][ipid]);
 395       }
 396     }
 397     return kSpKaon;
 398     break;
 399   case 211:
 400     if(FillQAHistos){
 401       for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
 402         if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
 403         TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",kSpPion,ipid));
 404         h->Fill(trk->Pt(),fnsigmas[kSpPion][ipid]);
 405       }
 406     }
 407     return kSpPion;
 408     break;
 409   default:
 410     return kSpUndefined;
 411   }
 412
 413 }
 414
 415 //////////////////////////////////////////////////////////////////////////////////////////////////
 416
 417 void AliHelperPID::CheckTOF(AliVTrack * trk)
 418 {
 419   //check if the particle has TOF Matching
 420   UInt_t status;
 421   status=trk->GetStatus();
 422   if((status&AliVTrack::kTOFout)==0 || (status&AliVTrack::kTIME)==0)fHasTOFPID=kFALSE;
 423   else fHasTOFPID=kTRUE;
 424   //in addition to KTOFout and kTIME we look at the pt
 425   if(trk->Pt()<fPtTOFPID)fHasTOFPID=kFALSE;
 426 }
 427
 428 //////////////////////////////////////////////////////////////////////////////////////////////////
 429
 430 Long64_t AliHelperPID::Merge(TCollection* list)
 431 {
 432   // Merging interface.
 433   // Returns the number of merged objects (including this).
 434   Printf("Merging");
 435   if (!list)
 436     return 0;
 437   if (list->IsEmpty())
 438     return 1;
 439   TIterator* iter = list->MakeIterator();
 440   TObject* obj;
 441   // collections of TList with output histos
 442   TList collections;
 443   Int_t count = 0;
 444   while ((obj = iter->Next())) {
 445     AliHelperPID* entry = dynamic_cast<AliHelperPID*> (obj);
 446     if (entry == 0)
 447       continue;
 448     TList * outputlist = entry->GetOutputList();
 449     collections.Add(outputlist);
 450     count++;
 451   }
 452   fOutputList->Merge(&collections);
 453   delete iter;
 454   Printf("OK");
 455   return count+1;
 456 }