[u/mrichter/AliRoot.git] / PWG / Tools / AliHelperPID.cxx


/**************************************************************************
 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//-----------------------------------------------------------------
//         AliHelperPID class
//-----------------------------------------------------------------

#include "AliHelperPID.h"
#include "AliVEvent.h"      
#include "AliAODEvent.h"      
#include "AliMCEvent.h"      
#include "AliMCEventHandler.h"      
#include "TH1F.h"             
#include "TH2F.h"             
#include "TList.h"            
#include "AliStack.h"            
#include "AliVTrack.h"      
#include "TParticle.h"
#include "AliAODMCParticle.h" 
#include "AliPIDResponse.h"   
#include "AliAnalysisManager.h"
#include "AliInputEventHandler.h"

using namespace AliHelperPIDNameSpace;
using namespace std;

const char * kPIDTypeName[]={"TPC","TOF","TPC-TOF"} ;
const char * kDetectorName[]={"TPC","TOF"} ;
const char * kParticleSpeciesName[]={"Pions","Kaons","Protons","Undefined"} ;

ClassImp(AliHelperPID)

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){
  
  for(Int_t ipart=0;ipart<kNSpecies;ipart++)
    for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++)
      fnsigmas[ipart][ipid]=999.;
  
  for(Int_t ipart=0;ipart<kNSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;
  
  fOutputList = new TList;
  fOutputList->SetOwner();
  fOutputList->SetName("OutputList");
  
  //nsigma plot
  for(Int_t ipart=0;ipart<kNSpecies;ipart++){
    for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
      Double_t miny=-30;
      Double_t maxy=30;
      if(ipid==kNSigmaTPCTOF){miny=0;maxy=50;}
      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);
      fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
      fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
      fOutputList->Add(fHistoNSigma);
    }
  }
  
  //nsigmaRec plot
  for(Int_t ipart=0;ipart<kNSpecies;ipart++){
    for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
      Double_t miny=-10;
      Double_t maxy=10;
      if(ipid==kNSigmaTPCTOF){miny=0;maxy=20;}
      TH2F *fHistoNSigma=new TH2F(Form("NSigmaRec_%d_%d",ipart,ipid),
				  Form("n#sigma for reconstructed %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,100,miny,maxy);
      fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
      fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
      fOutputList->Add(fHistoNSigma);
    }
  }
  
  //nsigmaDC plot
  for(Int_t ipart=0;ipart<kNSpecies;ipart++){
    for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
      Double_t miny=-10;
      Double_t maxy=10;
      if(ipid==kNSigmaTPCTOF){miny=0;maxy=20;}
      TH2F *fHistoNSigma=new TH2F(Form("NSigmaDC_%d_%d",ipart,ipid),
				  Form("n#sigma for double counting %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,100,miny,maxy);
      fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
      fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
      fOutputList->Add(fHistoNSigma);
    }
  }
  
  //nsigmaMC plot
  for(Int_t ipart=0;ipart<kNSpecies;ipart++){
    for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
      Double_t miny=-30;
      Double_t maxy=30;
      if(ipid==kNSigmaTPCTOF){miny=0;maxy=50;}
      TH2F *fHistoNSigma=new TH2F(Form("NSigmaMC_%d_%d",ipart,ipid),
				  Form("n#sigma for MC %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,100,miny,maxy);
      fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
      fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
      fOutputList->Add(fHistoNSigma);
    }
  }
  
  //PID signal plot
  for(Int_t idet=0;idet<kNDetectors;idet++){
    Double_t maxy=1000;
    if(idet==kTOF)maxy=200000;
    TH2F *fHistoPID=new TH2F(Form("PID_%d",idet),Form("%s signal",kDetectorName[idet]),200,0,10,200,-maxy,maxy);
    fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
    fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
    fOutputList->Add(fHistoPID);
  }
}

//////////////////////////////////////////////////////////////////////////////////////////////////

TH2F* AliHelperPID::GetHistogram2D(const char * name){
  // returns histo named name
  return (TH2F*) fOutputList->FindObject(name);
}

//////////////////////////////////////////////////////////////////////////////////////////////////

Int_t AliHelperPID::GetParticleSpecies(AliVTrack * trk, Bool_t FIllQAHistos){ 
  //return the specie according to the minimum nsigma value
  //no double counting, this has to be evaluated using CheckDoubleCounting()
  //Printf("fPtTOFPID %.1f, fRequestTOFPID %d, fNSigmaPID %.1f, fPIDType %d",fPtTOFPID,fRequestTOFPID,fNSigmaPID,fPIDType);
  
  CalculateNSigmas(trk,FIllQAHistos);//fill the data member fnsigmas with the nsigmas value [ipart][iPID]
  
  if(fUseExclusiveNSigma){
    Bool_t *HasDC;
    HasDC=GetDoubleCounting(trk,kFALSE);
    for(Int_t ipart=0;ipart<kNSpecies;ipart++){
      if(HasDC[ipart]==kTRUE)  return kSpUndefined;
    }
    return FindMinNSigma(trk,FIllQAHistos);//NSigmaRec distr filled here
  }
  else return FindMinNSigma(trk,FIllQAHistos);//NSigmaRec distr filled here
  
}

//////////////////////////////////////////////////////////////////////////////////////////////////

void AliHelperPID::CalculateNSigmas(AliVTrack * trk, Bool_t FIllQAHistos){ 
  //defines data member fnsigmas
  if(!fPIDResponse) {
    AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
    AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
    fPIDResponse = inputHandler->GetPIDResponse();
  }
  if(!fPIDResponse) {
    AliFatal("Cannot get pid response");
  }
  
  // Compute nsigma for each hypthesis
  AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(trk);
  // --- TPC
  Double_t nsigmaTPCkProton = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton);
  Double_t nsigmaTPCkKaon   = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon); 
  Double_t nsigmaTPCkPion   = fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion); 
  // --- TOF
  Double_t nsigmaTOFkProton=999.,nsigmaTOFkKaon=999.,nsigmaTOFkPion=999.;
  Double_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;

  CheckTOF(trk);
  
  if(fHasTOFPID && trk->Pt()>fPtTOFPID){//use TOF information
    nsigmaTOFkProton = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton);
    nsigmaTOFkKaon   = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon); 
    nsigmaTOFkPion   = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion); 
    // --- combined
    nsigmaTPCTOFkProton = TMath::Sqrt((nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton)/2.);
    nsigmaTPCTOFkKaon   = TMath::Sqrt((nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon)/2.);
    nsigmaTPCTOFkPion   = TMath::Sqrt((nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion)/2.);
  }else{
    // --- combined
    // if TOF is missing and below fPtTOFPID only the TPC information is used
    nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
    nsigmaTPCTOFkKaon   = TMath::Abs(nsigmaTPCkKaon);
    nsigmaTPCTOFkPion   = TMath::Abs(nsigmaTPCkPion);
  }
  
  //set data member fnsigmas
  fnsigmas[kSpPion][kNSigmaTPC]=nsigmaTPCkPion;
  fnsigmas[kSpKaon][kNSigmaTPC]=nsigmaTPCkKaon;
  fnsigmas[kSpProton][kNSigmaTPC]=nsigmaTPCkProton;
  fnsigmas[kSpPion][kNSigmaTOF]=nsigmaTOFkPion;
  fnsigmas[kSpKaon][kNSigmaTOF]=nsigmaTOFkKaon;
  fnsigmas[kSpProton][kNSigmaTOF]=nsigmaTOFkProton;
  fnsigmas[kSpPion][kNSigmaTPCTOF]=nsigmaTPCTOFkPion;
  fnsigmas[kSpKaon][kNSigmaTPCTOF]=nsigmaTPCTOFkKaon;
  fnsigmas[kSpProton][kNSigmaTPCTOF]=nsigmaTPCTOFkProton;
  
  if(FIllQAHistos){
    //Fill NSigma SeparationPlot
    for(Int_t ipart=0;ipart<kNSpecies;ipart++){
      for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
	if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && trk->Pt()<fPtTOFPID)continue;//not filling TOF and combined if no TOF PID
	TH2F *h=GetHistogram2D(Form("NSigma_%d_%d",ipart,ipid));
	h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
      }
    }
    //Fill PID signal plot
    for(Int_t idet=0;idet<kNDetectors;idet++){
      TH2F *h=GetHistogram2D(Form("PID_%d",idet));
      if(idet==kTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
      if(idet==kTOF && fHasTOFPID)h->Fill(trk->P(),trk->GetTOFsignal()*trk->Charge());
    }
  }
}

//////////////////////////////////////////////////////////////////////////////////////////////////

Int_t AliHelperPID::FindMinNSigma(AliVTrack * trk,Bool_t FillQAHistos){ 
  
  CheckTOF(trk);  
  if(fRequestTOFPID && (!fHasTOFPID) && trk->Pt()>fPtTOFPID)return kSpUndefined;
  
  //get the identity of the particle with the minimum Nsigma
  Double_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
  switch (fPIDType){
  case kNSigmaTPC:
    nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPC]);
    nsigmaKaon	  =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPC])  ;
    nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPC])  ;
    break;
  case kNSigmaTOF:
    nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTOF]);
    nsigmaKaon	  =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTOF])  ;
    nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTOF])  ;
    break;
  case kNSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
    nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
    nsigmaKaon	  =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF])  ;
    nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF])  ;
    break;
  }
  
  // guess the particle based on the smaller nsigma (within fNSigmaPID)
  if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return kSpUndefined;//if is the default value for the three
  
  if( ( nsigmaKaon   < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon   < fNSigmaPID)){
    if(FillQAHistos){
      for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
	if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
	TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",kSpKaon,ipid));
	h->Fill(trk->Pt(),fnsigmas[kSpKaon][ipid]);
      }
    }
    return kSpKaon;
  }
  if( ( nsigmaPion   < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion   < fNSigmaPID)){
    if(FillQAHistos){
      for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
	if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
	TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",kSpPion,ipid));
	h->Fill(trk->Pt(),fnsigmas[kSpPion][ipid]);
      }
    }
    return kSpPion;
  }
  if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)){
    if(FillQAHistos){
      for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
	if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
	TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",kSpProton,ipid));
	h->Fill(trk->Pt(),fnsigmas[kSpProton][ipid]);
      }
    }
    return kSpProton;
  }
  // else, return undefined
  return kSpUndefined;
}

//////////////////////////////////////////////////////////////////////////////////////////////////

Bool_t* AliHelperPID::GetDoubleCounting(AliVTrack * trk,Bool_t FIllQAHistos){ 
  //if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE
  //fill DC histos
  for(Int_t ipart=0;ipart<=kNSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track
  
  Int_t MinNSigma=FindMinNSigma(trk,kFALSE);//not filling the NSigmaRec histos
  
  CheckTOF(trk);
  
  if(MinNSigma==kSpUndefined)return fHasDoubleCounting;//in case of undefined no Double counting
  
  Double_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
  switch (fPIDType) {
  case kNSigmaTPC:
    nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPC]);
    nsigmaKaon	  =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPC])  ;
    nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPC])  ;
    break;
  case kNSigmaTOF:
    nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTOF]);
    nsigmaKaon	  =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTOF])  ;
    nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTOF])  ;
    break;
  case kNSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
    nsigmaProton  =  TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
    nsigmaKaon	  =  TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF])  ;
    nsigmaPion    =  TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF])  ;
    break;
  }
  if(nsigmaPion<fNSigmaPID && MinNSigma!=kSpPion)fHasDoubleCounting[kSpPion]=kTRUE;
  if(nsigmaKaon<fNSigmaPID && MinNSigma!=kSpKaon)fHasDoubleCounting[kSpKaon]=kTRUE;
  if(nsigmaProton<fNSigmaPID && MinNSigma!=kSpProton)fHasDoubleCounting[kSpProton]=kTRUE;
  
  if(FIllQAHistos){
    //fill NSigma distr for double counting
    for(Int_t ipart=0;ipart<kNSpecies;ipart++){
      if(fHasDoubleCounting[ipart]){
	for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
	  if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
	  TH2F *h=GetHistogram2D(Form("NSigmaDC_%d_%d",ipart,ipid));
	  h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
	}
      }
    }
  }
  return fHasDoubleCounting;
}

//////////////////////////////////////////////////////////////////////////////////////////////////

Int_t AliHelperPID::GetMCParticleSpecie(AliVEvent* event, AliVTrack * trk, Bool_t FillQAHistos){ 
  //return the specie according to the MC truth
  CheckTOF(trk);
  
  if(!fisMC)AliFatal("Error: AliHelperPID::GetMCParticleSpecie called on data\n");
  
  Int_t code=999;
  Bool_t isAOD=kFALSE;
  Bool_t isESD=kFALSE;
  TString nameoftrack(event->ClassName());  
  if(!nameoftrack.CompareTo("AliESDEvent"))isESD=kTRUE;
  else if(!nameoftrack.CompareTo("AliAODEvent"))isAOD=kTRUE;
  else AliFatal("Not processing AODs nor ESDs") ;
  
  if(isAOD){
    TClonesArray *arrayMC = 0;
    arrayMC = (TClonesArray*) event->GetList()->FindObject(AliAODMCParticle::StdBranchName());
    if (!arrayMC)AliFatal("Error: MC particles branch not found!\n");
    AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(trk->GetLabel()));
    if (!partMC){
      AliError("Cannot get MC particle");
      return kSpUndefined;
    }
    code=partMC->GetPdgCode();
  }
  if(isESD){
    AliStack* stack=0x0;
    AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
    if (!eventHandler) AliFatal("ERROR: Could not retrieve MC event handler");
    AliMCEvent* mcEvent = eventHandler->MCEvent();
    if (!mcEvent)AliFatal("ERROR: Could not retrieve MC event");
    stack = mcEvent->Stack();
    if (!stack) AliFatal("ERROR: stack not available\n");
    TParticle *part=0;
    part = (TParticle*)stack->Particle(TMath::Abs(trk->GetLabel()));
    if (!part){
      AliError("Cannot get MC particle");
      return kSpUndefined;
    }
    code = part->GetPdgCode();
  }
  switch(TMath::Abs(code)){
  case 2212:
    if(FillQAHistos){
      for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
	if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
	TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",kSpProton,ipid));
	h->Fill(trk->Pt(),fnsigmas[kSpProton][ipid]);
      }
    }
    return kSpProton;
    break;
  case 321:
    if(FillQAHistos){
      for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
	if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
	TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",kSpKaon,ipid));
	h->Fill(trk->Pt(),fnsigmas[kSpKaon][ipid]);
      }
    }
    return kSpKaon;
    break;
  case 211:
    if(FillQAHistos){
      for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
	if((ipid!=kNSigmaTPC) && (!fHasTOFPID) && (trk->Pt()<fPtTOFPID))continue;//not filling TOF and combined if no TOF PID
	TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",kSpPion,ipid));
	h->Fill(trk->Pt(),fnsigmas[kSpPion][ipid]);
      }
    }
    return kSpPion;
    break;
  default:
    return kSpUndefined;
  }

}

//////////////////////////////////////////////////////////////////////////////////////////////////

void AliHelperPID::CheckTOF(AliVTrack * trk)
{
  //check if the particle has TOF Matching
  UInt_t status;
  status=trk->GetStatus();
  if((status&AliVTrack::kTOFout)==0 || (status&AliVTrack::kTIME)==0)fHasTOFPID=kFALSE;
  else fHasTOFPID=kTRUE;
  //in addition to KTOFout and kTIME we look at the pt
  if(trk->Pt()<fPtTOFPID)fHasTOFPID=kFALSE;
  
  if(fRemoveTracksT0Fill)
    {
      //get the PIDResponse
      if(!fPIDResponse) {
	AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
	AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
	fPIDResponse = inputHandler->GetPIDResponse();
      }
      if(!fPIDResponse) {
	AliFatal("Cannot get pid response");
      }
      Int_t startTimeMask = fPIDResponse->GetTOFResponse().GetStartTimeMask(trk->P());
      if (startTimeMask < 0)fHasTOFPID=kFALSE; 
    }
}

//////////////////////////////////////////////////////////////////////////////////////////////////

Long64_t AliHelperPID::Merge(TCollection* list)
{
  // Merging interface.
  // Returns the number of merged objects (including this).
  Printf("Merging");
  if (!list)
    return 0;
  if (list->IsEmpty())
    return 1;
  TIterator* iter = list->MakeIterator();
  TObject* obj;
  // collections of TList with output histos
  TList collections;
  Int_t count = 0;
  while ((obj = iter->Next())) {
    AliHelperPID* entry = dynamic_cast<AliHelperPID*> (obj);
    if (entry == 0) 
      continue;
    TList * outputlist = entry->GetOutputList();      
    collections.Add(outputlist);
    count++;
  }
  fOutputList->Merge(&collections);
  delete iter;
  Printf("OK");
  return count+1;
}
Commit	Line	Data
32e49607	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
b46d749a	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){
32e49607	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
5838bc91	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	}
2b8b578d	146	return FindMinNSigma(trk,FIllQAHistos);//NSigmaRec distr filled here
5838bc91	147	}
5838bc91	148	else return FindMinNSigma(trk,FIllQAHistos);//NSigmaRec distr filled here
32e49607	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((nsigmaTPCkProtonnsigmaTPCkProton+nsigmaTOFkProtonnsigmaTOFkProton)/2.);
	183	nsigmaTPCTOFkKaon = TMath::Sqrt((nsigmaTPCkKaonnsigmaTPCkKaon+nsigmaTOFkKaonnsigmaTOFkKaon)/2.);
	184	nsigmaTPCTOFkPion = TMath::Sqrt((nsigmaTPCkPionnsigmaTPCkPion+nsigmaTOFkPionnsigmaTOFkPion)/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;
b46d749a	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	}
32e49607	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	}