[u/mrichter/AliRoot.git] / PWGPP / EvTrkSelection / AliSingleTrackEffCuts.cxx

#include "AliVEvent.h"
#include "AliMCEvent.h"
#include "TParticle.h"
#include "AliVParticle.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliVVertex.h"
#include "AliLog.h"
#include "AliGenEventHeader.h"
#include "AliInputEventHandler.h"
#include "AliAnalysisManager.h"
#include "AliAODMCHeader.h"
#include "AliAODMCParticle.h"
#include "AliMCParticle.h"
#include "AliVParticle.h"
#include "AliStack.h"
#include "AliMCEventHandler.h"
#include "AliPIDResponse.h"

#include "AliSingleTrackEffCuts.h"

using std::cout;
using std::endl;

ClassImp(AliSingleTrackEffCuts)

//_____________________________________________________
AliSingleTrackEffCuts::AliSingleTrackEffCuts():
AliAnalysisCuts(),
  fisAOD(kTRUE),
  fIsPdgCode(kFALSE),
  fPdgCode(0),
  fEtaMin(-12),
  fEtaMax(12),
  fYMin(-12),
  fYMax(12),
  fPtMin(-15),
  fPtMax(15),
  fIsCharged(kTRUE),
  fTriggerMask(AliVEvent::kAny),
  fMinVtxType(0),
  fMinVtxContr(1),
  fMaxVtxZ(10.),
  fCutOnZVertexSPD(0),
  fnClusITS(0),
  fnClusTPC(0),
  fnClusTOF(0),
  fnClusMUON(0),
  fusePid(kFALSE),
  fParticlePid(AliPID::kPion),
  fuseTPCPid(kTRUE),
  fnPTPCBins(0),
  fnPTPCBinLimits(0),
  fPTPCBinLimits(0),
  fnSigmaTPC(0),
  fPmaxTPC(9999.),
  fuseTOFPid(kTRUE),
  fnPTOFBins(0),
  fnPTOFBinLimits(0),
  fPTOFBinLimits(0),
  fnSigmaTOF(0),
  fPmaxTOF(9999.)
{
  //
  // Default constructor
  //
}

//________________________________________________________________________________
AliSingleTrackEffCuts::AliSingleTrackEffCuts(const char* name, const char* title):
AliAnalysisCuts(name,title),
  fisAOD(kTRUE),
  fIsPdgCode(kFALSE),
  fPdgCode(0),
  fEtaMin(-12),
  fEtaMax(12),
  fYMin(-12),
  fYMax(12),
  fPtMin(-15),
  fPtMax(15),
  fIsCharged(kTRUE),
  fTriggerMask(AliVEvent::kAny),
  fMinVtxType(0),
  fMinVtxContr(1),
  fMaxVtxZ(10.),
  fCutOnZVertexSPD(0),
  fnClusITS(0),
  fnClusTPC(0),
  fnClusTOF(0),
  fnClusMUON(0),
  fusePid(kFALSE),
  fParticlePid(AliPID::kPion),
  fuseTPCPid(kTRUE),
  fnPTPCBins(0),
  fnPTPCBinLimits(0),
  fPTPCBinLimits(0),
  fnSigmaTPC(0),
  fPmaxTPC(9999.),
  fuseTOFPid(kTRUE),
  fnPTOFBins(0),
  fnPTOFBinLimits(0),
  fPTOFBinLimits(0),
  fnSigmaTOF(0),
  fPmaxTOF(9999.)
{
  //
  // Default constructor
  //
}

//_________________________________________________________________________________
AliSingleTrackEffCuts::AliSingleTrackEffCuts(const AliSingleTrackEffCuts &source):
  AliAnalysisCuts(source),
  fisAOD(source.fisAOD),
  fIsPdgCode(source.fIsPdgCode),
  fPdgCode(source.fPdgCode),
  fEtaMin(source.fEtaMin),
  fEtaMax(source.fEtaMax),
  fYMin(source.fYMin),
  fYMax(source.fYMax),
  fPtMin(source.fPtMin),
  fPtMax(source.fPtMax),
  fIsCharged(source.fIsCharged),
  fTriggerMask(source.fTriggerMask),
  fMinVtxType(source.fMinVtxType),
  fMinVtxContr(source.fMinVtxContr),
  fMaxVtxZ(source.fMaxVtxZ),
  fCutOnZVertexSPD(source.fCutOnZVertexSPD),
  fnClusITS(source.fnClusITS),
  fnClusTPC(source.fnClusTPC),
  fnClusTOF(source.fnClusTOF),
  fnClusMUON(source.fnClusMUON),
  fusePid(source.fusePid),
  fParticlePid(source.fParticlePid),
  fuseTPCPid(source.fuseTPCPid),
  fnPTPCBins(source.fnPTPCBins),
  fnPTPCBinLimits(source.fnPTPCBinLimits),
  fPTPCBinLimits(source.fPTPCBinLimits),
  fnSigmaTPC(source.fnSigmaTPC),
  fPmaxTPC(source.fPmaxTPC),
  fuseTOFPid(source.fuseTOFPid),
  fnPTOFBins(source.fnPTOFBins),
  fnPTOFBinLimits(source.fnPTOFBinLimits),
  fPTOFBinLimits(source.fPTOFBinLimits),
  fnSigmaTOF(source.fnSigmaTOF),
  fPmaxTOF(source.fPmaxTOF)
{
  //
  // Copy constructor
  //
}


//_________________________________________________________________________________________
AliSingleTrackEffCuts &AliSingleTrackEffCuts::operator=(const AliSingleTrackEffCuts &source)
{
  //
  // assignment operator
  //
  if(&source == this) return *this;

  fisAOD = source.fisAOD;

  fIsPdgCode = source.fIsPdgCode;
  fPdgCode = source.fPdgCode;
  fEtaMin = source.fEtaMin;
  fEtaMax = source.fEtaMax;
  fYMin = source.fYMin;
  fYMax = source.fYMax;
  fPtMin = source.fPtMin;
  fPtMax = source.fPtMax;
  fIsCharged = source.fIsCharged;

  fTriggerMask = source.fTriggerMask;
  fMinVtxType = source.fMinVtxType;
  fMinVtxContr = source.fMinVtxContr;
  fMaxVtxZ = source.fMaxVtxZ;
  fCutOnZVertexSPD = source.fCutOnZVertexSPD;

  fnClusITS = source.fnClusITS;
  fnClusTPC = source.fnClusTPC;
  fnClusTOF = source.fnClusTOF;
  fnClusMUON = source.fnClusMUON;

  fusePid = source.fusePid;
  fParticlePid = source.fParticlePid;
  fuseTPCPid = source.fuseTPCPid;
  fnPTPCBins = source.fnPTPCBins;
  fnPTPCBinLimits = source.fnPTPCBinLimits;
  if(source.fPTPCBinLimits && source.fnSigmaTPC) 
    SetTPCSigmaPtBins(source.fnPTPCBins,source.fPTPCBinLimits,source.fnSigmaTPC);
  fPmaxTPC = source.fPmaxTPC;
  fuseTOFPid = source.fuseTOFPid;
  fnPTOFBins = source.fnPTOFBins;
  fnPTOFBinLimits = source.fnPTOFBinLimits;
  if(source.fPTOFBinLimits && source.fnSigmaTOF) 
    SetTOFSigmaPtBins(source.fnPTOFBins,source.fPTOFBinLimits,source.fnSigmaTOF);
  fPmaxTOF = source.fPmaxTOF;

  return *this;
}

//______________________________________________
AliSingleTrackEffCuts::~AliSingleTrackEffCuts()
{
  //
  // Destructor
  //

  if(fPTPCBinLimits){
    delete [] fPTPCBinLimits;
    fPTPCBinLimits=NULL;
  }
  if(fnSigmaTPC){
    delete [] fnSigmaTPC;
    fnSigmaTPC=NULL;
  }
  if(fPTOFBinLimits){
    delete [] fPTOFBinLimits;
    fPTOFBinLimits=NULL;
  }
  if(fnSigmaTOF){
    delete [] fnSigmaTOF;
    fnSigmaTOF=NULL;
  }

}

//______________________________________________
Bool_t AliSingleTrackEffCuts::IsMCEventSelected(TObject* obj)
{
  //
  // Event selection at MC level 
  //  |zvtx|<=fMaxVtxZ
  //
  if(!obj) return  kFALSE;
  Bool_t isSelected = kTRUE;
  
  AliMCEvent *event=0;
  AliGenEventHeader *genHeader=0; //ESDs
  AliAODMCHeader *mcHeader=0; // AODs
  Bool_t isAOD = obj->IsA()->InheritsFrom("AliAODEvent");

  if(!isAOD) {
    event = dynamic_cast<AliMCEvent*>(obj);
    if (!event) return  kFALSE;
    genHeader = event->GenEventHeader();
    if (!genHeader) return  kFALSE;
  } else {
    AliAODEvent *aodEvent = dynamic_cast<AliAODEvent*> (obj);
    if (!aodEvent) return  kFALSE;
    mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
    if (!mcHeader) {
      AliError("Could not find MC Header in AOD");
      return kFALSE;
    }

    // Check for not null trigger mask to evict pPb MC buggy-events
    Int_t runnumber = aodEvent->GetRunNumber();
    if(aodEvent->GetTriggerMask()==0 && 
       (runnumber>=195344 && runnumber<=195677)){
      AliDebug(3,"Event rejected because of null trigger mask");
      return kFALSE;
    }
  }


  // Cut on Z Vertex ( |z-vtx| <= fMaxVtxZ )
  TArrayF vtxPos(3);
  Double_t zMCVertex=-1000;
  if(!isAOD) {
    genHeader->PrimaryVertex(vtxPos);
    zMCVertex = vtxPos[2];
  } else {
    zMCVertex = mcHeader->GetVtxZ();
  }  
  if( TMath::Abs(zMCVertex)>fMaxVtxZ ) isSelected = kFALSE;


  return isSelected;
}

//__________________________________________________________
Bool_t AliSingleTrackEffCuts::IsMCParticleGenerated(TObject* obj)
{
  //
  // Check generated particles (primary, charged, pdgcode) 
  //
  
  if(!obj) return kFALSE;
  if(!obj->InheritsFrom("AliVParticle")) {
    AliError("object must derived from AliVParticle !");
    return kFALSE;
  }
  AliVParticle* particle = dynamic_cast<AliVParticle *>(obj);
  if(!particle) return kFALSE;

  Bool_t isSelected = kTRUE;

  // Check particle Pdg Code
  if(fIsPdgCode && TMath::Abs( particle->PdgCode() )!= fPdgCode) isSelected = kFALSE;

  // Charge selection
  if(fIsCharged && (particle->Charge()==0)) isSelected = kFALSE;

  // Selection of Physical Primary particles
  if(!fisAOD) { // check on ESDs
    AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());  	     
    AliMCEvent* mcevent = mcinfo->MCEvent();
    AliStack* stack = ((AliMCEvent*)mcevent)->Stack();
    AliMCParticle* mcPart = dynamic_cast<AliMCParticle *>(obj);	
    if (!mcPart) return kFALSE;
    if(!stack->IsPhysicalPrimary(mcPart->GetLabel())) {
      isSelected = kFALSE;
    }
  } else { // Check on AODs
    AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle *>(obj);	    
    if (!mcPart) return kFALSE;
    if(!mcPart->IsPhysicalPrimary()) isSelected = kFALSE;
  }

  return isSelected;
}


//_____________________________________________________________________
Bool_t AliSingleTrackEffCuts::IsMCParticleInKineAcceptance(TObject *obj)
{
  //
  // Check generated particles (eta, y, pt)
  //

  if(!obj) return  kFALSE;
  if(!obj->InheritsFrom("AliVParticle")) AliError("object must derived from AliVParticle !");
  AliVParticle* particle = dynamic_cast<AliVParticle *>(obj);
  if(!particle) return kFALSE;

  Bool_t isSelected = kTRUE;

  // Cut on eta
  if(particle->Eta()<fEtaMin || particle->Eta()>fEtaMax) isSelected = kFALSE;

  // Cut on y
  Double_t energy = particle->E();
  Double_t pz = particle->Pz();
  Double_t particleY = (TMath::Abs(energy-pz)>1e-10) ?  0.5*TMath::Log( (energy+pz)/(energy-pz) ) : 1e6;
  if(particleY<fYMin || particleY>fYMax) isSelected = kFALSE;

  // Cut on pt
  if(particle->Pt()<fPtMin || particle->Pt()>fPtMax) isSelected = kFALSE;

  return isSelected;
}


//_______________________________________________________________________
Bool_t AliSingleTrackEffCuts::IsMCParticleInReconstructable(TObject *obj)
{
  //
  // Check if particle has left enough hits in the detectors (only at ESD level)
  //

  if(!obj) return kFALSE;
  TString className(obj->ClassName());
  if (className.CompareTo("AliMCParticle") != 0) {
    AliError("obj must point to an AliMCParticle !");
    return kTRUE; // <===================================== FIX ME !!
  }

  AliMCParticle * part = dynamic_cast<AliMCParticle*>(obj);
  if(!part) return kFALSE;

  Bool_t isSelected = kTRUE;

  Int_t nHitsITS=0, nHitsTPC=0, nHitsTRD=0, nHitsTOF=0, nHitsMUON=0;
  for(Int_t iTrackRef=0; iTrackRef<part->GetNumberOfTrackReferences(); iTrackRef++) {
    AliTrackReference * trackRef = part->GetTrackReference(iTrackRef);
    if(trackRef){
      Int_t detectorId = trackRef->DetectorId();
      switch(detectorId) {
      case AliTrackReference::kITS  : 
	nHitsITS++; 
	break;
      case AliTrackReference::kTPC  : 
	nHitsTPC++; 
	break;
      case AliTrackReference::kTRD  : 
	nHitsTRD++; 
	break;
      case AliTrackReference::kTOF  : 
	nHitsTOF++; 
	break;
      case AliTrackReference::kMUON : 
	nHitsMUON++; 
	break;
      default : break;
      }
    }
  }

  if(nHitsITS<fnClusITS) isSelected = kFALSE;
  if(nHitsTPC<fnClusTPC) isSelected = kFALSE;
  if(nHitsTOF<fnClusTOF) isSelected = kFALSE;
  if(nHitsMUON<fnClusMUON) isSelected = kFALSE;

  return isSelected;
}


//_____________________________________________________________
Bool_t AliSingleTrackEffCuts::IsRecoEventSelected(TObject* obj)
{
  //
  // Event selection at reconstructed level (trigger, zvtx)
  //
 
  AliVEvent *event = dynamic_cast<AliVEvent*>(obj);
  if(!event) return kFALSE;

  Bool_t isSelected = kTRUE;

  // Check if event is accepted by the Physics selection
  UInt_t trigFired = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
  Bool_t isEvtSelected = (trigFired & fTriggerMask);
  if(!isEvtSelected) isSelected = kFALSE;

  // Vertex selection
  Bool_t isVtxSelected = IsVertexSelected(event);
  if(!isVtxSelected) isSelected = kFALSE;

  return isSelected;
}


//______________________________________________________________________
Bool_t AliSingleTrackEffCuts::IsRecoParticleKineAcceptance(TObject *obj)
{
  //
  // Check if reconstructed particle is in the acceptance (eta, y, pt)
  //

  Bool_t isSelected = kTRUE;

  AliVParticle *track = dynamic_cast<AliVParticle*>(obj);
  if (!track) return kFALSE;
 
  // Cut on eta
  if(track->Eta()<fEtaMin || track->Eta()>fEtaMax) isSelected = kFALSE;

  // Cut on y
  if(track->Y()<fYMin || track->Y()>fYMax) isSelected = kFALSE;

  // Cut on pt
  if(track->Pt()<fPtMin || track->Pt() >fPtMax) isSelected = kFALSE;

  return isSelected;
}

//_______________________________________________________________
Bool_t AliSingleTrackEffCuts::IsVertexSelected(AliVEvent *event)
{
  //
  // Check if the reconstructed vertex is selected
  //

  Bool_t accept = kTRUE;
  Bool_t isAOD = event->IsA()->InheritsFrom("AliAODEvent");

  const AliVVertex *vertex = event->GetPrimaryVertex();
  if(!vertex){
    accept = kFALSE;
    AliInfo("no vtx");
    return accept;
  }

  // Cut on vertex type  
  TString title=vertex->GetTitle();
  if(title.Contains("Z") && fMinVtxType>1){
    accept=kFALSE;
  } else if(title.Contains("3D") && fMinVtxType>2){
    accept=kFALSE;
  }
  
  // cut on minimum number of contributors
  if(vertex->GetNContributors()<fMinVtxContr){
    AliInfo(Form("too few contributors %d",vertex->GetNContributors()));
    accept=kFALSE;
  }

  // cut on absolute |z| of the vertex
  if(TMath::Abs(vertex->GetZ())>fMaxVtxZ) {
    AliInfo("outside the Vtx range");
    accept=kFALSE;
  }

  // cut on distance of SPD and TRK vertexes
  if(fCutOnZVertexSPD==1) {
    const AliVVertex *vSPD = NULL;
    if(isAOD) {
      vSPD = ((AliAODEvent*)event)->GetPrimaryVertexSPD(); 
    }else {    
      vSPD = ((AliESDEvent*)event)->GetPrimaryVertexSPD();  
    }

    if(vSPD && vSPD->GetNContributors()>=fMinVtxContr) {      
      if(TMath::Abs(vSPD->GetZ()-vertex->GetZ())>0.5) accept = kFALSE;
    }
  }

  return accept;
}

//_________________________________________________________________________________________________
void AliSingleTrackEffCuts::SetTPCSigmaPtBins(Int_t nPtBins, Float_t *pBinLimits, Float_t *sigmaBin)
{
  //
  // Set TPC Pid number-of-Sigma in P bins
  //

  // Set the pt bins
  if(fPTPCBinLimits) {
    delete [] fPTPCBinLimits;
    fPTPCBinLimits = NULL;
    printf("Changing the TPC cut P bins\n");
  }
  if(fnSigmaTPC) {
    delete [] fnSigmaTPC;
    fnSigmaTPC = NULL;
    printf("Changing the TPC Sigma cut per p bin\n");
  }

  fnPTPCBins = nPtBins;
  fnPTPCBinLimits = nPtBins+1;
  fPTPCBinLimits = new Float_t[fnPTPCBinLimits];
  for(Int_t ib=0; ib<nPtBins+1; ib++) fPTPCBinLimits[ib]=pBinLimits[ib];
  fnSigmaTPC = new Float_t[fnPTPCBins];
  for(Int_t ib=0; ib<nPtBins; ib++) fnSigmaTPC[ib]=sigmaBin[ib];

  return;
}
//_____________________________________________________________________________________________________
void AliSingleTrackEffCuts::SetTOFSigmaPtBins(Int_t nPtBins, Float_t *pBinLimits, Float_t *sigmaBin)
{
  //
  // Set TOF Pid number-of-Sigma in P bins
  //

  // Set the pt bins
  if(fPTOFBinLimits) {
    delete [] fPTOFBinLimits;
    fPTOFBinLimits = NULL;
    printf("Changing the TOF cut P bins\n");
  }
  if(fnSigmaTOF) {
    delete [] fnSigmaTOF;
    fnSigmaTOF = NULL;
    printf("Changing the TOF Sigma cut per p bin\n");
  }

  fnPTOFBins = nPtBins;
  fnPTOFBinLimits = nPtBins+1;
  fPTOFBinLimits = new Float_t[fnPTOFBinLimits];
  for(Int_t ib=0; ib<nPtBins+1; ib++) fPTOFBinLimits[ib]=pBinLimits[ib];
  fnSigmaTOF = new Float_t[fnPTOFBins];
  for(Int_t ib=0; ib<nPtBins; ib++) fnSigmaTOF[ib]=sigmaBin[ib];

  return;
}

//___________________________________________________________________________
Bool_t AliSingleTrackEffCuts::CheckTPCPIDStatus(AliAODTrack *track) const{
  //
  // Check TPC PID status
  //

  if ((track->GetStatus() & AliESDtrack::kTPCin)==0) return kFALSE;
  UShort_t nTPCClus=track->GetTPCClusterMap().CountBits();
  if (nTPCClus<70) return kFALSE;
  return kTRUE;
}

//___________________________________________________________________________
Bool_t AliSingleTrackEffCuts::CheckTOFPIDStatus(AliAODTrack *track) const{
  //
  // Check TOC PID status
  //

  if ((track->GetStatus()&AliESDtrack::kTOFout)==0)   return kFALSE;
  if ((track->GetStatus()&AliESDtrack::kTIME)==0)     return kFALSE;
  if ((track->GetStatus()&AliESDtrack::kTOFpid)==0)   return kFALSE;
  if ((track->GetStatus()&AliESDtrack::kTOFmismatch)!=0) return kFALSE;
  return kTRUE;
}

//___________________________________________________________________________
Bool_t AliSingleTrackEffCuts::IsRecoParticlePID(TObject *obj)
{
  //
  // Check Particle PID (AOD only for now!)
  //

  Bool_t isSelected = kFALSE;
  Bool_t isAOD = obj->IsA()->InheritsFrom("AliAODTrack");

  if(!isAOD || !GetUsePid()) return isSelected;
  if(!fuseTPCPid && !fuseTOFPid) return isSelected;

  AliAODTrack *track = dynamic_cast<AliAODTrack*>(obj);
  if(!track) { cout<<"No track found"<<endl; return isSelected; }

  // AliAODPid *pid = track->GetDetPid();
  // if(!pid) { cout<<"No AliAODPid found"<<endl; return isSelected; }

  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
  AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
  AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
  if(!pidResp) { cout<<"No PidResponse found"<<endl; return isSelected;}

  Double_t pPart = track->P();
  
  // Check detector status
  Bool_t okTPC = CheckTPCPIDStatus(track);
  Bool_t okTOF = CheckTOFPIDStatus(track);

  //  Check Number of Sigmas
  Double_t nsigmaTPC=pidResp->NumberOfSigmasTPC((AliVParticle*)track,(AliPID::EParticleType)fParticlePid);
  Double_t nsigmaTOF=pidResp->NumberOfSigmasTOF((AliVParticle*)track,(AliPID::EParticleType)fParticlePid);
  Bool_t isTPCPid=false, isTOFPid=false;

  // If use TPC and TPC infos are ok, check whether the sigma is ok in the given p range
  if(fuseTPCPid && okTPC) {
    for(Int_t j=0; j<fnPTPCBins; j++) {
      //      cout<<" checking bin: ("<<fPTPCBinLimits[j]<<","<<fPTPCBinLimits[j+1]<<") should be nsigma < "<<fnSigmaTPC[j]<<endl;
      if ((pPart>fPTPCBinLimits[j]) && (pPart<fPTPCBinLimits[j+1]) && nsigmaTPC<fnSigmaTPC[j]) isTPCPid=true;
    }
    if(pPart>fPmaxTPC) isTPCPid=true;
  }

  // If use TPC and TPC infos are ok, check whether the sigma is ok in the given p range
  if(fuseTOFPid && okTOF) {
    for(Int_t j=0; j<fnPTOFBins; j++) {
      if ((pPart>fPTOFBinLimits[j]) && (pPart<fPTOFBinLimits[j+1]) && nsigmaTOF<fnSigmaTOF[j]) isTPCPid=true;
    }
    if(pPart>fPmaxTOF) isTOFPid=true;
  }

  isSelected = (isTPCPid || isTOFPid) ? true : false;

  return isSelected;
}
Commit	Line	Data
c0757458	1	#include "AliVEvent.h"
	2	#include "AliMCEvent.h"
	3	#include "TParticle.h"
	4	#include "AliVParticle.h"
	5	#include "AliESDEvent.h"
	6	#include "AliAODEvent.h"
	7	#include "AliVVertex.h"
	8	#include "AliLog.h"
	9	#include "AliGenEventHeader.h"
	10	#include "AliInputEventHandler.h"
	11	#include "AliAnalysisManager.h"
	12	#include "AliAODMCHeader.h"
	13	#include "AliAODMCParticle.h"
	14	#include "AliMCParticle.h"
	15	#include "AliVParticle.h"
	16	#include "AliStack.h"
	17	#include "AliMCEventHandler.h"
	18	#include "AliPIDResponse.h"
	19
	20	#include "AliSingleTrackEffCuts.h"
	21
a0ec15ff	22	using std::cout;
a0ec15ff	23	using std::endl;
c0757458	24
	25	ClassImp(AliSingleTrackEffCuts)
	26
	27	//_____________________________________________________
	28	AliSingleTrackEffCuts::AliSingleTrackEffCuts():
	29	AliAnalysisCuts(),
	30	fisAOD(kTRUE),
	31	fIsPdgCode(kFALSE),
	32	fPdgCode(0),
	33	fEtaMin(-12),
	34	fEtaMax(12),
	35	fYMin(-12),
	36	fYMax(12),
	37	fPtMin(-15),
	38	fPtMax(15),
	39	fIsCharged(kTRUE),
	40	fTriggerMask(AliVEvent::kAny),
	41	fMinVtxType(0),
	42	fMinVtxContr(1),
	43	fMaxVtxZ(10.),
	44	fCutOnZVertexSPD(0),
	45	fnClusITS(0),
	46	fnClusTPC(0),
	47	fnClusTOF(0),
	48	fnClusMUON(0),
	49	fusePid(kFALSE),
	50	fParticlePid(AliPID::kPion),
	51	fuseTPCPid(kTRUE),
	52	fnPTPCBins(0),
	53	fnPTPCBinLimits(0),
	54	fPTPCBinLimits(0),
	55	fnSigmaTPC(0),
	56	fPmaxTPC(9999.),
	57	fuseTOFPid(kTRUE),
	58	fnPTOFBins(0),
	59	fnPTOFBinLimits(0),
	60	fPTOFBinLimits(0),
	61	fnSigmaTOF(0),
	62	fPmaxTOF(9999.)
	63	{
	64	//
	65	// Default constructor
	66	//
	67	}
	68
	69	//________________________________________________________________________________
	70	AliSingleTrackEffCuts::AliSingleTrackEffCuts(const char* name, const char* title):
	71	AliAnalysisCuts(name,title),
	72	fisAOD(kTRUE),
	73	fIsPdgCode(kFALSE),
	74	fPdgCode(0),
	75	fEtaMin(-12),
	76	fEtaMax(12),
	77	fYMin(-12),
	78	fYMax(12),
	79	fPtMin(-15),
	80	fPtMax(15),
	81	fIsCharged(kTRUE),
	82	fTriggerMask(AliVEvent::kAny),
	83	fMinVtxType(0),
	84	fMinVtxContr(1),
	85	fMaxVtxZ(10.),
	86	fCutOnZVertexSPD(0),
	87	fnClusITS(0),
88	fnClusTPC(0),
89	fnClusTOF(0),
90	fnClusMUON(0),
91	fusePid(kFALSE),
92	fParticlePid(AliPID::kPion),
93	fuseTPCPid(kTRUE),
94	fnPTPCBins(0),
95	fnPTPCBinLimits(0),
96	fPTPCBinLimits(0),
97	fnSigmaTPC(0),
98	fPmaxTPC(9999.),
99	fuseTOFPid(kTRUE),
100	fnPTOFBins(0),
101	fnPTOFBinLimits(0),
102	fPTOFBinLimits(0),
103	fnSigmaTOF(0),
104	fPmaxTOF(9999.)
105	{
106	//
107	// Default constructor
108	//
109	}
110
111	//_________________________________________________________________________________
112	AliSingleTrackEffCuts::AliSingleTrackEffCuts(const AliSingleTrackEffCuts &source):
113	AliAnalysisCuts(source),
114	fisAOD(source.fisAOD),
115	fIsPdgCode(source.fIsPdgCode),
116	fPdgCode(source.fPdgCode),
117	fEtaMin(source.fEtaMin),
118	fEtaMax(source.fEtaMax),
119	fYMin(source.fYMin),
120	fYMax(source.fYMax),
121	fPtMin(source.fPtMin),
122	fPtMax(source.fPtMax),
123	fIsCharged(source.fIsCharged),
124	fTriggerMask(source.fTriggerMask),
125	fMinVtxType(source.fMinVtxType),
126	fMinVtxContr(source.fMinVtxContr),
127	fMaxVtxZ(source.fMaxVtxZ),
128	fCutOnZVertexSPD(source.fCutOnZVertexSPD),
129	fnClusITS(source.fnClusITS),
130	fnClusTPC(source.fnClusTPC),
131	fnClusTOF(source.fnClusTOF),
132	fnClusMUON(source.fnClusMUON),
133	fusePid(source.fusePid),
134	fParticlePid(source.fParticlePid),
135	fuseTPCPid(source.fuseTPCPid),
136	fnPTPCBins(source.fnPTPCBins),
137	fnPTPCBinLimits(source.fnPTPCBinLimits),
138	fPTPCBinLimits(source.fPTPCBinLimits),
139	fnSigmaTPC(source.fnSigmaTPC),
140	fPmaxTPC(source.fPmaxTPC),
141	fuseTOFPid(source.fuseTOFPid),
142	fnPTOFBins(source.fnPTOFBins),
143	fnPTOFBinLimits(source.fnPTOFBinLimits),
144	fPTOFBinLimits(source.fPTOFBinLimits),
145	fnSigmaTOF(source.fnSigmaTOF),
146	fPmaxTOF(source.fPmaxTOF)
147	{
148	//
149	// Copy constructor
150	//
151	}
152
153
154	//_________________________________________________________________________________________
155	AliSingleTrackEffCuts &AliSingleTrackEffCuts::operator=(const AliSingleTrackEffCuts &source)
156	{
157	//
158	// assignment operator
159	//
160	if(&source == this) return *this;
161
162	fisAOD = source.fisAOD;
163
164	fIsPdgCode = source.fIsPdgCode;
165	fPdgCode = source.fPdgCode;
166	fEtaMin = source.fEtaMin;
167	fEtaMax = source.fEtaMax;
168	fYMin = source.fYMin;
169	fYMax = source.fYMax;
170	fPtMin = source.fPtMin;
171	fPtMax = source.fPtMax;
172	fIsCharged = source.fIsCharged;
173
174	fTriggerMask = source.fTriggerMask;
175	fMinVtxType = source.fMinVtxType;
176	fMinVtxContr = source.fMinVtxContr;
177	fMaxVtxZ = source.fMaxVtxZ;
178	fCutOnZVertexSPD = source.fCutOnZVertexSPD;
179
180	fnClusITS = source.fnClusITS;
181	fnClusTPC = source.fnClusTPC;
182	fnClusTOF = source.fnClusTOF;
183	fnClusMUON = source.fnClusMUON;
184
185	fusePid = source.fusePid;
186	fParticlePid = source.fParticlePid;
187	fuseTPCPid = source.fuseTPCPid;
188	fnPTPCBins = source.fnPTPCBins;
189	fnPTPCBinLimits = source.fnPTPCBinLimits;
190	if(source.fPTPCBinLimits && source.fnSigmaTPC)
191	SetTPCSigmaPtBins(source.fnPTPCBins,source.fPTPCBinLimits,source.fnSigmaTPC);
192	fPmaxTPC = source.fPmaxTPC;
193	fuseTOFPid = source.fuseTOFPid;
194	fnPTOFBins = source.fnPTOFBins;
195	fnPTOFBinLimits = source.fnPTOFBinLimits;
196	if(source.fPTOFBinLimits && source.fnSigmaTOF)
197	SetTOFSigmaPtBins(source.fnPTOFBins,source.fPTOFBinLimits,source.fnSigmaTOF);
198	fPmaxTOF = source.fPmaxTOF;
199
200	return *this;
201	}
202
203	//______________________________________________
204	AliSingleTrackEffCuts::~AliSingleTrackEffCuts()
205	{
206	//
207	// Destructor
208	//
209
210	if(fPTPCBinLimits){
211	delete [] fPTPCBinLimits;
212	fPTPCBinLimits=NULL;
213	}
214	if(fnSigmaTPC){
215	delete [] fnSigmaTPC;
216	fnSigmaTPC=NULL;
217	}
218	if(fPTOFBinLimits){
219	delete [] fPTOFBinLimits;
220	fPTOFBinLimits=NULL;
221	}
222	if(fnSigmaTOF){
223	delete [] fnSigmaTOF;
224	fnSigmaTOF=NULL;
225	}
226
227	}
228
229	//______________________________________________
230	Bool_t AliSingleTrackEffCuts::IsMCEventSelected(TObject* obj)
231	{
232	//
233	// Event selection at MC level
234	// \|zvtx\|<=fMaxVtxZ
235	//
236	if(!obj) return kFALSE;
237	Bool_t isSelected = kTRUE;
238
239	AliMCEvent *event=0;
240	AliGenEventHeader *genHeader=0; //ESDs
241	AliAODMCHeader *mcHeader=0; // AODs
242	Bool_t isAOD = obj->IsA()->InheritsFrom("AliAODEvent");
243
244	if(!isAOD) {
245	event = dynamic_cast<AliMCEvent*>(obj);
246	if (!event) return kFALSE;
247	genHeader = event->GenEventHeader();
248	if (!genHeader) return kFALSE;
249	} else {
250	AliAODEvent aodEvent = dynamic_cast<AliAODEvent> (obj);
251	if (!aodEvent) return kFALSE;
252	mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
253	if (!mcHeader) {
254	AliError("Could not find MC Header in AOD");
255	return kFALSE;
256	}
257
258	// Check for not null trigger mask to evict pPb MC buggy-events
259	Int_t runnumber = aodEvent->GetRunNumber();
260	if(aodEvent->GetTriggerMask()==0 &&
261	(runnumber>=195344 && runnumber<=195677)){
262	AliDebug(3,"Event rejected because of null trigger mask");
263	return kFALSE;
264	}
265	}
266
267
268	// Cut on Z Vertex ( \|z-vtx\| <= fMaxVtxZ )
269	TArrayF vtxPos(3);
270	Double_t zMCVertex=-1000;
271	if(!isAOD) {
272	genHeader->PrimaryVertex(vtxPos);
273	zMCVertex = vtxPos[2];
274	} else {
275	zMCVertex = mcHeader->GetVtxZ();
276	}
277	if( TMath::Abs(zMCVertex)>fMaxVtxZ ) isSelected = kFALSE;
278
279
280	return isSelected;
281	}
282
283	//__________________________________________________________
284	Bool_t AliSingleTrackEffCuts::IsMCParticleGenerated(TObject* obj)
285	{
286	//
287	// Check generated particles (primary, charged, pdgcode)
288	//
289
290	if(!obj) return kFALSE;
291	if(!obj->InheritsFrom("AliVParticle")) {
292	AliError("object must derived from AliVParticle !");
293	return kFALSE;
294	}
295	AliVParticle* particle = dynamic_cast<AliVParticle *>(obj);
296	if(!particle) return kFALSE;
297
298	Bool_t isSelected = kTRUE;
299
300	// Check particle Pdg Code
301	if(fIsPdgCode && TMath::Abs( particle->PdgCode() )!= fPdgCode) isSelected = kFALSE;
302
303	// Charge selection
64aa1114	304	if(fIsCharged && (particle->Charge()==0)) isSelected = kFALSE;
c0757458	305
	306	// Selection of Physical Primary particles
	307	if(!fisAOD) { // check on ESDs
	308	AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
	309	AliMCEvent* mcevent = mcinfo->MCEvent();
	310	AliStack* stack = ((AliMCEvent*)mcevent)->Stack();
c37908f3	311	AliMCParticle* mcPart = dynamic_cast<AliMCParticle *>(obj);
c37908f3	312	if (!mcPart) return kFALSE;
c0757458	313	if(!stack->IsPhysicalPrimary(mcPart->GetLabel())) {
	314	isSelected = kFALSE;
	315	}
	316	} else { // Check on AODs
	317	AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle *>(obj);
77a76537	318	if (!mcPart) return kFALSE;
c0757458	319	if(!mcPart->IsPhysicalPrimary()) isSelected = kFALSE;
	320	}
	321
	322	return isSelected;
	323	}
	324
	325
	326	//_____________________________________________________________________
	327	Bool_t AliSingleTrackEffCuts::IsMCParticleInKineAcceptance(TObject *obj)
	328	{
	329	//
	330	// Check generated particles (eta, y, pt)
	331	//
	332
	333	if(!obj) return kFALSE;
	334	if(!obj->InheritsFrom("AliVParticle")) AliError("object must derived from AliVParticle !");
	335	AliVParticle* particle = dynamic_cast<AliVParticle *>(obj);
	336	if(!particle) return kFALSE;
	337
	338	Bool_t isSelected = kTRUE;
	339
	340	// Cut on eta
	341	if(particle->Eta()<fEtaMin \|\| particle->Eta()>fEtaMax) isSelected = kFALSE;
	342
	343	// Cut on y
	344	Double_t energy = particle->E();
	345	Double_t pz = particle->Pz();
	346	Double_t particleY = (TMath::Abs(energy-pz)>1e-10) ? 0.5*TMath::Log( (energy+pz)/(energy-pz) ) : 1e6;
	347	if(particleY<fYMin \|\| particleY>fYMax) isSelected = kFALSE;
	348
	349	// Cut on pt
	350	if(particle->Pt()<fPtMin \|\| particle->Pt()>fPtMax) isSelected = kFALSE;
	351
	352	return isSelected;
	353	}
	354
	355
	356	//_______________________________________________________________________
	357	Bool_t AliSingleTrackEffCuts::IsMCParticleInReconstructable(TObject *obj)
	358	{
	359	//
	360	// Check if particle has left enough hits in the detectors (only at ESD level)
	361	//
	362
	363	if(!obj) return kFALSE;
	364	TString className(obj->ClassName());
	365	if (className.CompareTo("AliMCParticle") != 0) {
	366	AliError("obj must point to an AliMCParticle !");
	367	return kTRUE; // <===================================== FIX ME !!
	368	}
	369
	370	AliMCParticle * part = dynamic_cast<AliMCParticle*>(obj);
	371	if(!part) return kFALSE;
	372
	373	Bool_t isSelected = kTRUE;
	374
	375	Int_t nHitsITS=0, nHitsTPC=0, nHitsTRD=0, nHitsTOF=0, nHitsMUON=0;
	376	for(Int_t iTrackRef=0; iTrackRef<part->GetNumberOfTrackReferences(); iTrackRef++) {
	377	AliTrackReference * trackRef = part->GetTrackReference(iTrackRef);
	378	if(trackRef){
	379	Int_t detectorId = trackRef->DetectorId();
	380	switch(detectorId) {
5cc7b66d	381	case AliTrackReference::kITS :
	382	nHitsITS++;
	383	break;
	384	case AliTrackReference::kTPC :
	385	nHitsTPC++;
	386	break;
	387	case AliTrackReference::kTRD :
	388	nHitsTRD++;
	389	break;
	390	case AliTrackReference::kTOF :
	391	nHitsTOF++;
	392	break;
	393	case AliTrackReference::kMUON :
	394	nHitsMUON++;
	395	break;
	396	default : break;
c0757458	397	}
	398	}
	399	}
	400
	401	if(nHitsITS<fnClusITS) isSelected = kFALSE;
	402	if(nHitsTPC<fnClusTPC) isSelected = kFALSE;
	403	if(nHitsTOF<fnClusTOF) isSelected = kFALSE;
	404	if(nHitsMUON<fnClusMUON) isSelected = kFALSE;
	405
	406	return isSelected;
	407	}
	408
	409
	410	//_____________________________________________________________
	411	Bool_t AliSingleTrackEffCuts::IsRecoEventSelected(TObject* obj)
	412	{
	413	//
	414	// Event selection at reconstructed level (trigger, zvtx)
	415	//
	416
	417	AliVEvent event = dynamic_cast<AliVEvent>(obj);
	418	if(!event) return kFALSE;
	419
	420	Bool_t isSelected = kTRUE;
	421
	422	// Check if event is accepted by the Physics selection
	423	UInt_t trigFired = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
	424	Bool_t isEvtSelected = (trigFired & fTriggerMask);
	425	if(!isEvtSelected) isSelected = kFALSE;
	426
	427	// Vertex selection
	428	Bool_t isVtxSelected = IsVertexSelected(event);
	429	if(!isVtxSelected) isSelected = kFALSE;
	430
	431	return isSelected;
	432	}
	433
	434
	435	//______________________________________________________________________
	436	Bool_t AliSingleTrackEffCuts::IsRecoParticleKineAcceptance(TObject *obj)
	437	{
	438	//
	439	// Check if reconstructed particle is in the acceptance (eta, y, pt)
	440	//
	441
	442	Bool_t isSelected = kTRUE;
	443
	444	AliVParticle track = dynamic_cast<AliVParticle>(obj);
c37908f3	445	if (!track) return kFALSE;
c0757458	446
	447	// Cut on eta
	448	if(track->Eta()<fEtaMin \|\| track->Eta()>fEtaMax) isSelected = kFALSE;
	449
	450	// Cut on y
	451	if(track->Y()<fYMin \|\| track->Y()>fYMax) isSelected = kFALSE;
	452
	453	// Cut on pt
	454	if(track->Pt()<fPtMin \|\| track->Pt() >fPtMax) isSelected = kFALSE;
	455
	456	return isSelected;
	457	}
	458
	459	//_______________________________________________________________
	460	Bool_t AliSingleTrackEffCuts::IsVertexSelected(AliVEvent *event)
	461	{
	462	//
	463	// Check if the reconstructed vertex is selected
	464	//
	465
	466	Bool_t accept = kTRUE;
	467	Bool_t isAOD = event->IsA()->InheritsFrom("AliAODEvent");
	468
	469	const AliVVertex *vertex = event->GetPrimaryVertex();
	470	if(!vertex){
	471	accept = kFALSE;
	472	AliInfo("no vtx");
	473	return accept;
	474	}
	475
	476	// Cut on vertex type
	477	TString title=vertex->GetTitle();
	478	if(title.Contains("Z") && fMinVtxType>1){
	479	accept=kFALSE;
	480	} else if(title.Contains("3D") && fMinVtxType>2){
	481	accept=kFALSE;
	482	}
	483
	484	// cut on minimum number of contributors
	485	if(vertex->GetNContributors()<fMinVtxContr){
	486	AliInfo(Form("too few contributors %d",vertex->GetNContributors()));
	487	accept=kFALSE;
	488	}
	489
	490	// cut on absolute \|z\| of the vertex
	491	if(TMath::Abs(vertex->GetZ())>fMaxVtxZ) {
	492	AliInfo("outside the Vtx range");
	493	accept=kFALSE;
	494	}
	495
	496	// cut on distance of SPD and TRK vertexes
	497	if(fCutOnZVertexSPD==1) {
	498	const AliVVertex *vSPD = NULL;
	499	if(isAOD) {
	500	vSPD = ((AliAODEvent*)event)->GetPrimaryVertexSPD();
	501	}else {
	502	vSPD = ((AliESDEvent*)event)->GetPrimaryVertexSPD();
	503	}
	504
	505	if(vSPD && vSPD->GetNContributors()>=fMinVtxContr) {
	506	if(TMath::Abs(vSPD->GetZ()-vertex->GetZ())>0.5) accept = kFALSE;
	507	}
	508	}
	509
510	return accept;
511	}
512
513	//_________________________________________________________________________________________________
514	void AliSingleTrackEffCuts::SetTPCSigmaPtBins(Int_t nPtBins, Float_t pBinLimits, Float_t sigmaBin)
515	{
516	//
517	// Set TPC Pid number-of-Sigma in P bins
518	//
519
520	// Set the pt bins
521	if(fPTPCBinLimits) {
522	delete [] fPTPCBinLimits;
523	fPTPCBinLimits = NULL;
524	printf("Changing the TPC cut P bins\n");
525	}
526	if(fnSigmaTPC) {
527	delete [] fnSigmaTPC;
528	fnSigmaTPC = NULL;
529	printf("Changing the TPC Sigma cut per p bin\n");
530	}
531
532	fnPTPCBins = nPtBins;
533	fnPTPCBinLimits = nPtBins+1;
534	fPTPCBinLimits = new Float_t[fnPTPCBinLimits];
535	for(Int_t ib=0; ib<nPtBins+1; ib++) fPTPCBinLimits[ib]=pBinLimits[ib];
536	fnSigmaTPC = new Float_t[fnPTPCBins];
537	for(Int_t ib=0; ib<nPtBins; ib++) fnSigmaTPC[ib]=sigmaBin[ib];
538
539	return;
540	}
541	//_____________________________________________________________________________________________________
542	void AliSingleTrackEffCuts::SetTOFSigmaPtBins(Int_t nPtBins, Float_t pBinLimits, Float_t sigmaBin)
543	{
544	//
545	// Set TOF Pid number-of-Sigma in P bins
546	//
547
548	// Set the pt bins
549	if(fPTOFBinLimits) {
550	delete [] fPTOFBinLimits;
551	fPTOFBinLimits = NULL;
552	printf("Changing the TOF cut P bins\n");
553	}
554	if(fnSigmaTOF) {
555	delete [] fnSigmaTOF;
556	fnSigmaTOF = NULL;
557	printf("Changing the TOF Sigma cut per p bin\n");
558	}
559
560	fnPTOFBins = nPtBins;
561	fnPTOFBinLimits = nPtBins+1;
562	fPTOFBinLimits = new Float_t[fnPTOFBinLimits];
563	for(Int_t ib=0; ib<nPtBins+1; ib++) fPTOFBinLimits[ib]=pBinLimits[ib];
564	fnSigmaTOF = new Float_t[fnPTOFBins];
565	for(Int_t ib=0; ib<nPtBins; ib++) fnSigmaTOF[ib]=sigmaBin[ib];
566
567	return;
568	}
569
570	//___________________________________________________________________________
571	Bool_t AliSingleTrackEffCuts::CheckTPCPIDStatus(AliAODTrack *track) const{
572	//
573	// Check TPC PID status
574	//
575
576	if ((track->GetStatus() & AliESDtrack::kTPCin)==0) return kFALSE;
577	UShort_t nTPCClus=track->GetTPCClusterMap().CountBits();
578	if (nTPCClus<70) return kFALSE;
579	return kTRUE;
580	}
581
582	//___________________________________________________________________________
583	Bool_t AliSingleTrackEffCuts::CheckTOFPIDStatus(AliAODTrack *track) const{
584	//
585	// Check TOC PID status
586	//
587
588	if ((track->GetStatus()&AliESDtrack::kTOFout)==0) return kFALSE;
589	if ((track->GetStatus()&AliESDtrack::kTIME)==0) return kFALSE;
590	if ((track->GetStatus()&AliESDtrack::kTOFpid)==0) return kFALSE;
64aa1114	591	if ((track->GetStatus()&AliESDtrack::kTOFmismatch)!=0) return kFALSE;
c0757458	592	return kTRUE;
	593	}
	594
	595	//___________________________________________________________________________
	596	Bool_t AliSingleTrackEffCuts::IsRecoParticlePID(TObject *obj)
	597	{
	598	//
	599	// Check Particle PID (AOD only for now!)
	600	//
	601
	602	Bool_t isSelected = kFALSE;
	603	Bool_t isAOD = obj->IsA()->InheritsFrom("AliAODTrack");
	604
	605	if(!isAOD \|\| !GetUsePid()) return isSelected;
	606	if(!fuseTPCPid && !fuseTOFPid) return isSelected;
	607
	608	AliAODTrack track = dynamic_cast<AliAODTrack>(obj);
	609	if(!track) { cout<<"No track found"<<endl; return isSelected; }
	610
	611	// AliAODPid *pid = track->GetDetPid();
	612	// if(!pid) { cout<<"No AliAODPid found"<<endl; return isSelected; }
	613
	614	AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
	615	AliInputEventHandler inputHandler=(AliInputEventHandler)mgr->GetInputEventHandler();
	616	AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
	617	if(!pidResp) { cout<<"No PidResponse found"<<endl; return isSelected;}
	618
	619	Double_t pPart = track->P();
	620
	621	// Check detector status
	622	Bool_t okTPC = CheckTPCPIDStatus(track);
	623	Bool_t okTOF = CheckTOFPIDStatus(track);
	624
	625	// Check Number of Sigmas
	626	Double_t nsigmaTPC=pidResp->NumberOfSigmasTPC((AliVParticle*)track,(AliPID::EParticleType)fParticlePid);
	627	Double_t nsigmaTOF=pidResp->NumberOfSigmasTOF((AliVParticle*)track,(AliPID::EParticleType)fParticlePid);
	628	Bool_t isTPCPid=false, isTOFPid=false;
	629
	630	// If use TPC and TPC infos are ok, check whether the sigma is ok in the given p range
	631	if(fuseTPCPid && okTPC) {
	632	for(Int_t j=0; j<fnPTPCBins; j++) {
	633	// cout<<" checking bin: ("<<fPTPCBinLimits[j]<<","<<fPTPCBinLimits[j+1]<<") should be nsigma < "<<fnSigmaTPC[j]<<endl;
	634	if ((pPart>fPTPCBinLimits[j]) && (pPart<fPTPCBinLimits[j+1]) && nsigmaTPC<fnSigmaTPC[j]) isTPCPid=true;
	635	}
	636	if(pPart>fPmaxTPC) isTPCPid=true;
	637	}
	638
	639	// If use TPC and TPC infos are ok, check whether the sigma is ok in the given p range
	640	if(fuseTOFPid && okTOF) {
	641	for(Int_t j=0; j<fnPTOFBins; j++) {
	642	if ((pPart>fPTOFBinLimits[j]) && (pPart<fPTOFBinLimits[j+1]) && nsigmaTOF<fnSigmaTOF[j]) isTPCPid=true;
	643	}
	644	if(pPart>fPmaxTOF) isTOFPid=true;
	645	}
	646
	647	isSelected = (isTPCPid \|\| isTOFPid) ? true : false;
	648
	649	return isSelected;
	650	}