[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronPID.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.                  *
**************************************************************************/

///////////////////////////////////////////////////////////////////////////
//                Dielectron PID                                  //
//                                                                       //
//                                                                       //
/*
Detailed description


*/
//                                                                       //
///////////////////////////////////////////////////////////////////////////

#include <TMath.h>
#include <TF1.h>
#include <TGraph.h>

#include <AliVTrack.h>
#include <AliVCluster.h>
#include <AliLog.h>
#include <AliExternalTrackParam.h>
#include <AliPIDResponse.h>
#include <AliESDtrack.h> //!!!!! Remove once Eta correction is treated in the tender

#include "AliDielectronVarManager.h"

#include "AliDielectronPID.h"

ClassImp(AliDielectronPID)

TGraph *AliDielectronPID::fgFitCorr=0x0;
Double_t AliDielectronPID::fgCorr=0.0;
Double_t AliDielectronPID::fgCorrdEdx=1.0;
TF1 *AliDielectronPID::fgFunEtaCorr=0x0;
TGraph *AliDielectronPID::fgdEdxRunCorr=0x0;

AliDielectronPID::AliDielectronPID() :
  AliAnalysisCuts(),
  fNcuts(0),
  fPIDResponse(0x0)
{
  //
  // Default Constructor
  //
  for (Int_t icut=0; icut<kNmaxPID; ++icut){
    fDetType[icut]=kTPC;
    fPartType[icut]=AliPID::kPion;
    fNsigmaLow[icut]=0;
    fNsigmaUp[icut]=0;
    fmin[icut]=0;
    fmax[icut]=0;
    fExclude[icut]=kFALSE;
    fFunUpperCut[icut]=0x0;
    fFunLowerCut[icut]=0x0;
    fRequirePIDbit[icut]=0;
    fActiveCuts[icut]=-1;
  }
}

//______________________________________________
AliDielectronPID::AliDielectronPID(const char* name, const char* title) :
  AliAnalysisCuts(name, title),
  fNcuts(0),
  fPIDResponse(0x0)
{
  //
  // Named Constructor
  //
  for (Int_t icut=0; icut<kNmaxPID; ++icut){
    fDetType[icut]=kTPC;
    fPartType[icut]=AliPID::kPion;
    fNsigmaLow[icut]=0;
    fNsigmaUp[icut]=0;
    fmin[icut]=0;
    fmax[icut]=0;
    fExclude[icut]=kFALSE;
    fFunUpperCut[icut]=0x0;
    fFunLowerCut[icut]=0x0;
    fRequirePIDbit[icut]=0;
    fActiveCuts[icut]=0;
  }
}

//______________________________________________
AliDielectronPID::~AliDielectronPID()
{
  //
  // Default Destructor
  //
}

//______________________________________________
void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, Double_t nSigmaUp/*=-99999.*/,
                              Double_t min/*=0*/, Double_t max/*=0*/, Bool_t exclude/*=kFALSE*/,
                              UInt_t pidBitType/*=AliDielectronPID::kRequire*/, Int_t var/*=-1*/)
{
  //
  // Add a pid nsigma cut
  // use response of detector 'det' in the range ['min'] to ['max'] for var
  // use a sigma band between 'nSigmaLow' and 'nSigmaUp'
  // if nSigmaUp==-99999. then nSigmaLow will be uesd as a symmetric band +- nSigmaLow
  // specify whether to 'exclude' the given band
  //

  if (fNcuts>=kNmaxPID){
    AliError(Form("only %d pid cut ranges allowed",kNmaxPID));
    return;
  }
  if (TMath::Abs(nSigmaUp+99999.)<1e-20){
    nSigmaUp=TMath::Abs(nSigmaLow);
    nSigmaLow=-1.*nSigmaUp;
  }
  fDetType[fNcuts]=det;
  fPartType[fNcuts]=type;
  fNsigmaLow[fNcuts]=nSigmaLow;
  fNsigmaUp[fNcuts]=nSigmaUp;
  fmin[fNcuts]=min;
  fmax[fNcuts]=max;
  fExclude[fNcuts]=exclude;
  fRequirePIDbit[fNcuts]=pidBitType;
  fActiveCuts[fNcuts]=(var==-1 ? AliDielectronVarManager::kP : var);

  AliInfo(Form("Add PID cut %d: sigma [% .1f,% .1f] \t cut [% .1f,% .f] \t var %d->%s \n",
	       fNcuts,nSigmaLow,nSigmaUp,min,max,fActiveCuts[fNcuts],AliDielectronVarManager::GetValueName(fActiveCuts[fNcuts])));  
  
  ++fNcuts;

}

//______________________________________________
void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, TF1 * const funUp,
                              Double_t min/*=0*/, Double_t max/*=0*/, Bool_t exclude/*=kFALSE*/,
                              UInt_t pidBitType/*=AliDielectronPID::kRequire*/, Int_t var/*=-1*/)
{
  //
  // cut using a TF1 as upper cut
  //
  if (funUp==0x0){
    AliError("A valid function is required for the upper cut. Not adding the cut!");
    return;
  }
  fFunUpperCut[fNcuts]=funUp;
  AddCut(det,type,nSigmaLow,0.,min,max,exclude,pidBitType,var);
}

//______________________________________________
void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, TF1 * const funLow, Double_t nSigmaUp,
                              Double_t min/*=0*/, Double_t max/*=0*/, Bool_t exclude/*=kFALSE*/,
                              UInt_t pidBitType/*=AliDielectronPID::kRequire*/, Int_t var/*=-1*/)
{
  //
  // cut using a TF1 as lower cut
  //
  if (funLow==0x0){
    AliError("A valid function is required for the lower cut. Not adding the cut!");
    return;
  }
  fFunLowerCut[fNcuts]=funLow;
  AddCut(det,type,0.,nSigmaUp,min,max,exclude,pidBitType,var);
}

//______________________________________________
void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, TF1 * const funLow, TF1 * const funUp,
                              Double_t min/*=0*/, Double_t max/*=0*/, Bool_t exclude/*=kFALSE*/,
                              UInt_t pidBitType/*=AliDielectronPID::kRequire*/, Int_t var/*=-1*/)
{
  //
  // cut using a TF1 as lower and upper cut
  //
  if ( (funUp==0x0) || (funLow==0x0) ){
    AliError("A valid function is required for upper and lower cut. Not adding the cut!");
    return;
  }
  fFunUpperCut[fNcuts]=funUp;
  fFunLowerCut[fNcuts]=funLow;
  AddCut(det,type,0.,0.,min,max,exclude,pidBitType,var);
}

//______________________________________________
Bool_t AliDielectronPID::IsSelected(TObject* track)
{
  //
  // perform PID cuts
  //

  //loop over all cuts
  AliVTrack *part=static_cast<AliVTrack*>(track);
  AliESDtrack *esdTrack=0x0;
  Double_t origdEdx=-1;
  
  // apply ETa correction, remove once this is in the tender
  if( (part->IsA() == AliESDtrack::Class()) && fgFunEtaCorr ){
    esdTrack=static_cast<AliESDtrack*>(part);
    origdEdx=esdTrack->GetTPCsignal();
    esdTrack->SetTPCsignal(origdEdx/GetEtaCorr(esdTrack)/fgCorrdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN());
  }

  
  //Fill values
  Double_t values[AliDielectronVarManager::kNMaxValues];
  AliDielectronVarManager::Fill(track,values);

  Bool_t selected=kFALSE;
  fPIDResponse=AliDielectronVarManager::GetPIDResponse();
  for (UChar_t icut=0; icut<fNcuts; ++icut){
    Double_t min=fmin[icut];
    Double_t max=fmax[icut];
    Double_t val=values[fActiveCuts[icut]];
    
    // test var range. In case min==max do not cut
    if ( (TMath::Abs(min-max)>1e-20) && (val<min || val>=max) ) continue;
    
    switch (fDetType[icut]){
    case kITS:
      selected = IsSelectedITS(part,icut);
      break;
    case kTPC:
      selected = IsSelectedTPC(part,icut);
      break;
    case kTRD:
      selected = IsSelectedTRD(part,icut);
      break;
    case kTRDeleEff:
      selected = IsSelectedTRDeleEff(part,icut);
      break;
    case kTOF:
      selected = IsSelectedTOF(part,icut);
      break;
    case kEMCAL:
      selected = IsSelectedEMCAL(part,icut);
      break;
    }
    if (!selected) {
      if (esdTrack) esdTrack->SetTPCsignal(origdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN());

      return kFALSE;
    }
  }

  if (esdTrack) esdTrack->SetTPCsignal(origdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN());
  return selected;
}

//______________________________________________
Bool_t AliDielectronPID::IsSelectedITS(AliVTrack * const part, Int_t icut)
{
  //
  // ITS part of the PID check
  // Don't accept the track if there was no pid bit set
  //
  
  if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliVTrack::kITSpid)) return kFALSE;
  if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliVTrack::kITSpid)) return kTRUE;

  Double_t mom=part->P();
  
  Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasITS(part, fPartType[icut]);
  
  // test if we are supposed to use a function for the cut
  if (fFunUpperCut[icut]) fNsigmaUp[icut] =fFunUpperCut[icut]->Eval(mom);
  if (fFunLowerCut[icut]) fNsigmaLow[icut]=fFunLowerCut[icut]->Eval(mom);
  
  Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
  return selected;
}

//______________________________________________
Bool_t AliDielectronPID::IsSelectedTPC(AliVTrack * const part, Int_t icut)
{
  //
  // TPC part of the PID check
  // Don't accept the track if there was no pid bit set
  //
  if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliVTrack::kTPCpid)) return kFALSE;
  if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliVTrack::kTPCpid)) return kTRUE;

  Double_t mom=part->GetTPCmomentum();
  
  Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasTPC(part, fPartType[icut]);

  if (fPartType[icut]==AliPID::kElectron){
    numberOfSigmas-=fgCorr;
  }
  
  // test if we are supposed to use a function for the cut
  if (fFunUpperCut[icut]) fNsigmaUp[icut] =fFunUpperCut[icut]->Eval(mom);
  if (fFunLowerCut[icut]) fNsigmaLow[icut]=fFunLowerCut[icut]->Eval(mom);
  
  Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
  return selected;
}

//______________________________________________
Bool_t AliDielectronPID::IsSelectedTRD(AliVTrack * const part, Int_t icut)
{
  //   
  // TRD part of the pid check
  // the TRD checks on the probabilities.
  //

  if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliVTrack::kTRDpid)) return kFALSE;
  if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliVTrack::kTRDpid)) return kTRUE;
  if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable && (part->GetTRDntrackletsPID()<4)) return kTRUE;

  Double_t p[AliPID::kSPECIES]={0.};
  fPIDResponse->ComputeTRDProbability(part,AliPID::kSPECIES,p);
  Float_t particleProb=p[fPartType[icut]];

  Bool_t selected=((particleProb>=fNsigmaLow[icut])&&(particleProb<=fNsigmaUp[icut]))^fExclude[icut];
  return selected;
}

//______________________________________________
Bool_t AliDielectronPID::IsSelectedTRDeleEff(AliVTrack * const part, Int_t icut)
{
  //
  // TRD part of the pid check using electron efficiency requirement
  // in this case the upper limit as well as the particle specie is ignored 
  //   and the lower limit regarded as the requested electron efficiency
  //

  if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliVTrack::kTRDpid)) return kFALSE;
  if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliVTrack::kTRDpid)) return kTRUE;

  Bool_t selected=fPIDResponse->IdentifiedAsElectronTRD(part,fNsigmaLow[icut])^fExclude[icut];
  return selected;
}

//______________________________________________
Bool_t AliDielectronPID::IsSelectedTOF(AliVTrack * const part, Int_t icut)
{
  //
  // TOF part of the PID check
  // Don't accept the track if there was no pid bit set
  //
  
  if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliVTrack::kTOFpid)) return kFALSE;
  if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliVTrack::kTOFpid)) return kTRUE;

  Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasTOF(part, fPartType[icut]);
  
  Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
  return selected;
}

//______________________________________________
Bool_t AliDielectronPID::IsSelectedEMCAL(AliVTrack * const part, Int_t icut)
{
  //
  // emcal pid selecttion
  //

  //TODO: correct way to check for emcal pid?
  Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasEMCAL(part, fPartType[icut]);

  Bool_t hasPID=numberOfSigmas>-998.;

  if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!hasPID) return kFALSE;
  if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!hasPID) return kTRUE;


  Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
  return selected;
}

//______________________________________________
void AliDielectronPID::SetDefaults(Int_t def){
  //
  // initialise default pid strategies
  //

  if (def==0){
    // 2sigma bands TPC:
    // - include e
    // - exclude mu,K,pi,p
    // -complete p range
    AddCut(kTPC,AliPID::kElectron,2);
    AddCut(kTPC,AliPID::kMuon,-2.,2.,0.,0.,kTRUE);
    AddCut(kTPC,AliPID::kPion,-2.,2.,0.,0.,kTRUE);
    AddCut(kTPC,AliPID::kKaon,-2.,2.,0.,0.,kTRUE);
    AddCut(kTPC,AliPID::kProton,-2.,2.,0.,0.,kTRUE);
  } else if (def==1) {
    // 2sigma bands TPC:
    // - include e 0<p<inf
    // - exclude mu,K,pi,p 0<p<2
    AddCut(kTPC,AliPID::kElectron,2.);
    AddCut(kTPC,AliPID::kMuon,-2.,2.,0.,2.,kTRUE);
    AddCut(kTPC,AliPID::kPion,-2.,2.,0.,2.,kTRUE);
    AddCut(kTPC,AliPID::kKaon,-2.,2.,0.,2.,kTRUE);
    AddCut(kTPC,AliPID::kProton,-2.,2.,0.,2.,kTRUE);
    
  } else if (def==2) {
    // include 2sigma e TPC
    // 3sigma bands TOF
    // - exclude K,P
    AddCut(kTPC,AliPID::kElectron,-3.,3.);
    AddCut(kTPC,AliPID::kPion,-3.,3.,0.,0.,kTRUE);
    AddCut(kTPC,AliPID::kProton,-3.,3.,0.,0.,kTRUE);

  } else if (def==3 || def==4) { // def3 and def4 are the same??
    // include 2sigma e TPC
    // 3sigma bands TOF
    // - exclude K 0<p<1
    // - exclude P 0<p<2
    AddCut(kTPC,AliPID::kElectron,2);
    AddCut(kTOF,AliPID::kKaon,-3.,3.,0.,1.,kTRUE);
    AddCut(kTOF,AliPID::kProton,-6.,6.,0.,1.,kTRUE);
    AddCut(kTOF,AliPID::kProton,-3.,3.,1.,2.,kTRUE);
    
  } else if (def==5) {
    AddCut(kTPC,AliPID::kElectron,-0.5,3);
    AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
  } else if (def==6) {
    // lower cut TPC: parametrisation by HFE
    // upper cut TPC: 3 sigma
    // TOF ele band 3sigma 0<p<1.5GeV
    TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)", 0, 100);
    lowerCut->SetParameters(-2.7,-0.4357);
    AddCut(kTPC,AliPID::kElectron,lowerCut,3.);
    AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
  } else if (def==7) {
    // wide TPC cut
    // TOF ele band 3sigma 0<p<1.5GeV
    AddCut(kTPC,AliPID::kElectron,10.);
    AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
  } else if (def==8) {
    // TOF 5 sigma inclusion if TOFpid available
    // this should reduce K,p,Pi to a large extent
    AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
  } else if (def==9) {
    // lower cut TPC: parametrisation by HFE
    // upper cut TPC: 3 sigma
    // TOF 5 sigma inclusion if TOFpid available
    // this should reduce K,p,Pi to a large extent
    TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)", 0, 100);
    lowerCut->SetParameters(-2.65,-0.6757);
    AddCut(kTPC,AliPID::kElectron,lowerCut,4.);
    AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);

  } else if (def==10) {
    AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
    AddCut(kTPC,AliPID::kElectron,3.);
    AddCut(kTPC,AliPID::kPion,-3.,3.,0.,0.,kTRUE);
    AddCut(kTPC,AliPID::kProton,-3.,3.,0.,0.,kTRUE);
    
  } else if (def==11) {
    // lower cut TPC: parametrisation by HFE
    // only use from period d on !!
    // upper cut TPC: 3 sigma
    // TOF ele band 3sigma 0<p<2.0GeV
    TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)+[2]", 0, 100);
    lowerCut->SetParameters(-3.718,-0.4,0.27);
    AddCut(kTPC,AliPID::kElectron,lowerCut,3.);
    AddCut(kTOF,AliPID::kElectron,-3,3,0,2.);
  } else if (def==12) {
    // lower cut TPC: parametrisation by HFE
    // only use from period d on !!
    // upper cut TPC: 3 sigma
    // TOF 5 sigma inclusion if TOFpid available
    // this should reduce K,p,Pi to a large extent
    TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)+[2]", 0, 100);
    lowerCut->SetParameters(-3.718,-0.4,0.27);
    AddCut(kTPC,AliPID::kElectron,lowerCut,4.);
    AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
  }

}

//______________________________________________
void AliDielectronPID::SetCorrVal(Double_t run)
{
  //
  // set correction value for run
  //
  fgCorr=0.;
  fgCorrdEdx=1.;
  
  if (fgFitCorr){
    fgCorr=fgFitCorr->Eval(run);
    if (run<fgFitCorr->GetX()[0]) fgCorr=fgFitCorr->GetY()[0];
    if (run>fgFitCorr->GetX()[fgFitCorr->GetN()-1]) fgCorr=fgFitCorr->GetY()[fgFitCorr->GetN()-1];
  }

  if (fgdEdxRunCorr){
    fgCorrdEdx=fgdEdxRunCorr->Eval(run);
    if (run<fgdEdxRunCorr->GetX()[0]) fgCorrdEdx=fgdEdxRunCorr->GetY()[0];
    if (run>fgdEdxRunCorr->GetX()[fgFitCorr->GetN()-1]) fgCorrdEdx=fgdEdxRunCorr->GetY()[fgdEdxRunCorr->GetN()-1];
  }
}

//______________________________________________
Double_t AliDielectronPID::GetEtaCorr(AliVTrack *track)
{
  //
  // return eta correction
  //
  if (!fgFunEtaCorr) return 1;
  return fgFunEtaCorr->Eval(track->Eta());
}
Commit	Line	Data
8df8e382	1	/*************************************************************************
	2	* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	///////////////////////////////////////////////////////////////////////////
	17	// Dielectron PID //
	18	// //
	19	// //
	20	/*
	21	Detailed description
	22
	23
	24	*/
	25	// //
	26	///////////////////////////////////////////////////////////////////////////
	27
	28	#include <TMath.h>
	29	#include <TF1.h>
48609e3d	30	#include <TGraph.h>
8df8e382	31
164bfb53	32	#include <AliVTrack.h>
5720c765	33	#include <AliVCluster.h>
8df8e382	34	#include <AliLog.h>
ffbede40	35	#include <AliExternalTrackParam.h>
5720c765	36	#include <AliPIDResponse.h>
5720c765	37	#include <AliESDtrack.h> //!!!!! Remove once Eta correction is treated in the tender
8df8e382	38
	39	#include "AliDielectronVarManager.h"
	40
	41	#include "AliDielectronPID.h"
	42
	43	ClassImp(AliDielectronPID)
	44
48609e3d	45	TGraph *AliDielectronPID::fgFitCorr=0x0;
48609e3d	46	Double_t AliDielectronPID::fgCorr=0.0;
d327d9cd	47	Double_t AliDielectronPID::fgCorrdEdx=1.0;
5720c765	48	TF1 *AliDielectronPID::fgFunEtaCorr=0x0;
d327d9cd	49	TGraph *AliDielectronPID::fgdEdxRunCorr=0x0;
48609e3d	50
8df8e382	51	AliDielectronPID::AliDielectronPID() :
	52	AliAnalysisCuts(),
	53	fNcuts(0),
5720c765	54	fPIDResponse(0x0)
8df8e382	55	{
	56	//
	57	// Default Constructor
	58	//
	59	for (Int_t icut=0; icut<kNmaxPID; ++icut){
	60	fDetType[icut]=kTPC;
	61	fPartType[icut]=AliPID::kPion;
	62	fNsigmaLow[icut]=0;
	63	fNsigmaUp[icut]=0;
5720c765	64	fmin[icut]=0;
5720c765	65	fmax[icut]=0;
8df8e382	66	fExclude[icut]=kFALSE;
	67	fFunUpperCut[icut]=0x0;
	68	fFunLowerCut[icut]=0x0;
61d106d3	69	fRequirePIDbit[icut]=0;
5720c765	70	fActiveCuts[icut]=-1;
8df8e382	71	}
	72	}
	73
	74	//______________________________________________
	75	AliDielectronPID::AliDielectronPID(const char* name, const char* title) :
	76	AliAnalysisCuts(name, title),
	77	fNcuts(0),
5720c765	78	fPIDResponse(0x0)
8df8e382	79	{
	80	//
	81	// Named Constructor
	82	//
	83	for (Int_t icut=0; icut<kNmaxPID; ++icut){
	84	fDetType[icut]=kTPC;
	85	fPartType[icut]=AliPID::kPion;
	86	fNsigmaLow[icut]=0;
	87	fNsigmaUp[icut]=0;
5720c765	88	fmin[icut]=0;
5720c765	89	fmax[icut]=0;
8df8e382	90	fExclude[icut]=kFALSE;
	91	fFunUpperCut[icut]=0x0;
	92	fFunLowerCut[icut]=0x0;
61d106d3	93	fRequirePIDbit[icut]=0;
5720c765	94	fActiveCuts[icut]=0;
8df8e382	95	}
	96	}
	97
	98	//______________________________________________
	99	AliDielectronPID::~AliDielectronPID()
	100	{
	101	//
	102	// Default Destructor
	103	//
	104	}
	105
	106	//______________________________________________
	107	void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, Double_t nSigmaUp/=-99999./,
5720c765	108	Double_t min/=0/, Double_t max/=0/, Bool_t exclude/=kFALSE/,
5720c765	109	UInt_t pidBitType/=AliDielectronPID::kRequire/, Int_t var/=-1/)
8df8e382	110	{
	111	//
	112	// Add a pid nsigma cut
5720c765	113	// use response of detector 'det' in the range ['min'] to ['max'] for var
5720c765	114	// use a sigma band between 'nSigmaLow' and 'nSigmaUp'
8df8e382	115	// if nSigmaUp==-99999. then nSigmaLow will be uesd as a symmetric band +- nSigmaLow
	116	// specify whether to 'exclude' the given band
	117	//
	118
45b2b1b8	119	if (fNcuts>=kNmaxPID){
8df8e382	120	AliError(Form("only %d pid cut ranges allowed",kNmaxPID));
45b2b1b8	121	return;
8df8e382	122	}
	123	if (TMath::Abs(nSigmaUp+99999.)<1e-20){
	124	nSigmaUp=TMath::Abs(nSigmaLow);
	125	nSigmaLow=-1.*nSigmaUp;
	126	}
	127	fDetType[fNcuts]=det;
	128	fPartType[fNcuts]=type;
	129	fNsigmaLow[fNcuts]=nSigmaLow;
	130	fNsigmaUp[fNcuts]=nSigmaUp;
5720c765	131	fmin[fNcuts]=min;
5720c765	132	fmax[fNcuts]=max;
8df8e382	133	fExclude[fNcuts]=exclude;
61d106d3	134	fRequirePIDbit[fNcuts]=pidBitType;
5720c765	135	fActiveCuts[fNcuts]=(var==-1 ? AliDielectronVarManager::kP : var);
	136
	137	AliInfo(Form("Add PID cut %d: sigma [% .1f,% .1f] \t cut [% .1f,% .f] \t var %d->%s \n",
	138	fNcuts,nSigmaLow,nSigmaUp,min,max,fActiveCuts[fNcuts],AliDielectronVarManager::GetValueName(fActiveCuts[fNcuts])));
8df8e382	139
5720c765	140	++fNcuts;
5720c765	141
8df8e382	142	}
	143
	144	//______________________________________________
	145	void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, TF1 * const funUp,
5720c765	146	Double_t min/=0/, Double_t max/=0/, Bool_t exclude/=kFALSE/,
5720c765	147	UInt_t pidBitType/=AliDielectronPID::kRequire/, Int_t var/=-1/)
8df8e382	148	{
	149	//
	150	// cut using a TF1 as upper cut
	151	//
	152	if (funUp==0x0){
	153	AliError("A valid function is required for the upper cut. Not adding the cut!");
	154	return;
	155	}
	156	fFunUpperCut[fNcuts]=funUp;
5720c765	157	AddCut(det,type,nSigmaLow,0.,min,max,exclude,pidBitType,var);
8df8e382	158	}
	159
	160	//______________________________________________
	161	void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, TF1 * const funLow, Double_t nSigmaUp,
5720c765	162	Double_t min/=0/, Double_t max/=0/, Bool_t exclude/=kFALSE/,
5720c765	163	UInt_t pidBitType/=AliDielectronPID::kRequire/, Int_t var/=-1/)
8df8e382	164	{
	165	//
	166	// cut using a TF1 as lower cut
	167	//
	168	if (funLow==0x0){
	169	AliError("A valid function is required for the lower cut. Not adding the cut!");
	170	return;
	171	}
	172	fFunLowerCut[fNcuts]=funLow;
5720c765	173	AddCut(det,type,0.,nSigmaUp,min,max,exclude,pidBitType,var);
8df8e382	174	}
	175
	176	//______________________________________________
	177	void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, TF1 * const funLow, TF1 * const funUp,
5720c765	178	Double_t min/=0/, Double_t max/=0/, Bool_t exclude/=kFALSE/,
5720c765	179	UInt_t pidBitType/=AliDielectronPID::kRequire/, Int_t var/=-1/)
8df8e382	180	{
	181	//
	182	// cut using a TF1 as lower and upper cut
	183	//
	184	if ( (funUp==0x0) \|\| (funLow==0x0) ){
	185	AliError("A valid function is required for upper and lower cut. Not adding the cut!");
	186	return;
	187	}
	188	fFunUpperCut[fNcuts]=funUp;
	189	fFunLowerCut[fNcuts]=funLow;
5720c765	190	AddCut(det,type,0.,0.,min,max,exclude,pidBitType,var);
8df8e382	191	}
	192
	193	//______________________________________________
	194	Bool_t AliDielectronPID::IsSelected(TObject* track)
	195	{
	196	//
	197	// perform PID cuts
	198	//
	199
	200	//loop over all cuts
164bfb53	201	AliVTrack part=static_cast<AliVTrack>(track);
5720c765	202	AliESDtrack *esdTrack=0x0;
5720c765	203	Double_t origdEdx=-1;
8df8e382	204
5720c765	205	// apply ETa correction, remove once this is in the tender
	206	if( (part->IsA() == AliESDtrack::Class()) && fgFunEtaCorr ){
	207	esdTrack=static_cast<AliESDtrack*>(part);
	208	origdEdx=esdTrack->GetTPCsignal();
d327d9cd	209	esdTrack->SetTPCsignal(origdEdx/GetEtaCorr(esdTrack)/fgCorrdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN());
5720c765	210	}
8df8e382	211
5720c765	212
	213	//Fill values
	214	Double_t values[AliDielectronVarManager::kNMaxValues];
	215	AliDielectronVarManager::Fill(track,values);
8df8e382	216
5720c765	217	Bool_t selected=kFALSE;
	218	fPIDResponse=AliDielectronVarManager::GetPIDResponse();
	219	for (UChar_t icut=0; icut<fNcuts; ++icut){
	220	Double_t min=fmin[icut];
	221	Double_t max=fmax[icut];
	222	Double_t val=values[fActiveCuts[icut]];
	223
	224	// test var range. In case min==max do not cut
	225	if ( (TMath::Abs(min-max)>1e-20) && (val<min \|\| val>=max) ) continue;
	226
8df8e382	227	switch (fDetType[icut]){
	228	case kITS:
	229	selected = IsSelectedITS(part,icut);
	230	break;
	231	case kTPC:
	232	selected = IsSelectedTPC(part,icut);
	233	break;
	234	case kTRD:
	235	selected = IsSelectedTRD(part,icut);
	236	break;
5720c765	237	case kTRDeleEff:
	238	selected = IsSelectedTRDeleEff(part,icut);
	239	break;
8df8e382	240	case kTOF:
	241	selected = IsSelectedTOF(part,icut);
	242	break;
5720c765	243	case kEMCAL:
	244	selected = IsSelectedEMCAL(part,icut);
	245	break;
	246	}
	247	if (!selected) {
	248	if (esdTrack) esdTrack->SetTPCsignal(origdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN());
	249
	250	return kFALSE;
8df8e382	251	}
8df8e382	252	}
8df8e382	253
5720c765	254	if (esdTrack) esdTrack->SetTPCsignal(origdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN());
8df8e382	255	return selected;
	256	}
	257
	258	//______________________________________________
ffbede40	259	Bool_t AliDielectronPID::IsSelectedITS(AliVTrack * const part, Int_t icut)
8df8e382	260	{
	261	//
	262	// ITS part of the PID check
	263	// Don't accept the track if there was no pid bit set
	264	//
8df8e382	265
5720c765	266	if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliVTrack::kITSpid)) return kFALSE;
5720c765	267	if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliVTrack::kITSpid)) return kTRUE;
ffbede40	268
ffbede40	269	Double_t mom=part->P();
164bfb53	270
5720c765	271	Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasITS(part, fPartType[icut]);
ffbede40	272
	273	// test if we are supposed to use a function for the cut
	274	if (fFunUpperCut[icut]) fNsigmaUp[icut] =fFunUpperCut[icut]->Eval(mom);
	275	if (fFunLowerCut[icut]) fNsigmaLow[icut]=fFunLowerCut[icut]->Eval(mom);
	276
8df8e382	277	Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
	278	return selected;
	279	}
	280
	281	//______________________________________________
ffbede40	282	Bool_t AliDielectronPID::IsSelectedTPC(AliVTrack * const part, Int_t icut)
8df8e382	283	{
	284	//
	285	// TPC part of the PID check
	286	// Don't accept the track if there was no pid bit set
	287	//
5720c765	288	if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliVTrack::kTPCpid)) return kFALSE;
5720c765	289	if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliVTrack::kTPCpid)) return kTRUE;
48609e3d	290
5720c765	291	Double_t mom=part->GetTPCmomentum();
ffbede40	292
5720c765	293	Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasTPC(part, fPartType[icut]);
5720c765	294
2a14a7b1	295	if (fPartType[icut]==AliPID::kElectron){
2a14a7b1	296	numberOfSigmas-=fgCorr;
8df8e382	297	}
ffbede40	298
	299	// test if we are supposed to use a function for the cut
	300	if (fFunUpperCut[icut]) fNsigmaUp[icut] =fFunUpperCut[icut]->Eval(mom);
	301	if (fFunLowerCut[icut]) fNsigmaLow[icut]=fFunLowerCut[icut]->Eval(mom);
	302
8df8e382	303	Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
	304	return selected;
	305	}
	306
	307	//______________________________________________
5720c765	308	Bool_t AliDielectronPID::IsSelectedTRD(AliVTrack * const part, Int_t icut)
8df8e382	309	{
	310	//
	311	// TRD part of the pid check
5720c765	312	// the TRD checks on the probabilities.
	313	//
	314
	315	if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliVTrack::kTRDpid)) return kFALSE;
	316	if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliVTrack::kTRDpid)) return kTRUE;
	317	if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable && (part->GetTRDntrackletsPID()<4)) return kTRUE;
	318
	319	Double_t p[AliPID::kSPECIES]={0.};
	320	fPIDResponse->ComputeTRDProbability(part,AliPID::kSPECIES,p);
	321	Float_t particleProb=p[fPartType[icut]];
	322
	323	Bool_t selected=((particleProb>=fNsigmaLow[icut])&&(particleProb<=fNsigmaUp[icut]))^fExclude[icut];
	324	return selected;
	325	}
	326
	327	//______________________________________________
	328	Bool_t AliDielectronPID::IsSelectedTRDeleEff(AliVTrack * const part, Int_t icut)
	329	{
8df8e382	330	//
5720c765	331	// TRD part of the pid check using electron efficiency requirement
	332	// in this case the upper limit as well as the particle specie is ignored
	333	// and the lower limit regarded as the requested electron efficiency
	334	//
	335
	336	if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliVTrack::kTRDpid)) return kFALSE;
	337	if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliVTrack::kTRDpid)) return kTRUE;
	338
	339	Bool_t selected=fPIDResponse->IdentifiedAsElectronTRD(part,fNsigmaLow[icut])^fExclude[icut];
	340	return selected;
8df8e382	341	}
	342
	343	//______________________________________________
ffbede40	344	Bool_t AliDielectronPID::IsSelectedTOF(AliVTrack * const part, Int_t icut)
8df8e382	345	{
	346	//
	347	// TOF part of the PID check
	348	// Don't accept the track if there was no pid bit set
	349	//
8df8e382	350
5720c765	351	if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliVTrack::kTOFpid)) return kFALSE;
	352	if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliVTrack::kTOFpid)) return kTRUE;
	353
	354	Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasTOF(part, fPartType[icut]);
8df8e382	355
	356	Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
	357	return selected;
	358	}
	359
5720c765	360	//______________________________________________
	361	Bool_t AliDielectronPID::IsSelectedEMCAL(AliVTrack * const part, Int_t icut)
	362	{
	363	//
	364	// emcal pid selecttion
	365	//
	366
	367	//TODO: correct way to check for emcal pid?
	368	Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasEMCAL(part, fPartType[icut]);
	369
	370	Bool_t hasPID=numberOfSigmas>-998.;
	371
	372	if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!hasPID) return kFALSE;
	373	if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!hasPID) return kTRUE;
	374
	375
	376	Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
	377	return selected;
	378	}
	379
8df8e382	380	//______________________________________________
	381	void AliDielectronPID::SetDefaults(Int_t def){
	382	//
	383	// initialise default pid strategies
	384	//
	385
	386	if (def==0){
	387	// 2sigma bands TPC:
	388	// - include e
	389	// - exclude mu,K,pi,p
	390	// -complete p range
	391	AddCut(kTPC,AliPID::kElectron,2);
	392	AddCut(kTPC,AliPID::kMuon,-2.,2.,0.,0.,kTRUE);
	393	AddCut(kTPC,AliPID::kPion,-2.,2.,0.,0.,kTRUE);
	394	AddCut(kTPC,AliPID::kKaon,-2.,2.,0.,0.,kTRUE);
	395	AddCut(kTPC,AliPID::kProton,-2.,2.,0.,0.,kTRUE);
	396	} else if (def==1) {
	397	// 2sigma bands TPC:
	398	// - include e 0<p<inf
	399	// - exclude mu,K,pi,p 0<p<2
8df8e382	400	AddCut(kTPC,AliPID::kElectron,2.);
	401	AddCut(kTPC,AliPID::kMuon,-2.,2.,0.,2.,kTRUE);
	402	AddCut(kTPC,AliPID::kPion,-2.,2.,0.,2.,kTRUE);
	403	AddCut(kTPC,AliPID::kKaon,-2.,2.,0.,2.,kTRUE);
	404	AddCut(kTPC,AliPID::kProton,-2.,2.,0.,2.,kTRUE);
	405
	406	} else if (def==2) {
	407	// include 2sigma e TPC
	408	// 3sigma bands TOF
	409	// - exclude K,P
48609e3d	410	AddCut(kTPC,AliPID::kElectron,-3.,3.);
3505bfad	411	AddCut(kTPC,AliPID::kPion,-3.,3.,0.,0.,kTRUE);
3505bfad	412	AddCut(kTPC,AliPID::kProton,-3.,3.,0.,0.,kTRUE);
8df8e382	413
5720c765	414	} else if (def==3 \|\| def==4) { // def3 and def4 are the same??
8df8e382	415	// include 2sigma e TPC
	416	// 3sigma bands TOF
	417	// - exclude K 0<p<1
	418	// - exclude P 0<p<2
	419	AddCut(kTPC,AliPID::kElectron,2);
	420	AddCut(kTOF,AliPID::kKaon,-3.,3.,0.,1.,kTRUE);
	421	AddCut(kTOF,AliPID::kProton,-6.,6.,0.,1.,kTRUE);
	422	AddCut(kTOF,AliPID::kProton,-3.,3.,1.,2.,kTRUE);
	423
8df8e382	424	} else if (def==5) {
	425	AddCut(kTPC,AliPID::kElectron,-0.5,3);
	426	AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
	427	} else if (def==6) {
	428	// lower cut TPC: parametrisation by HFE
	429	// upper cut TPC: 3 sigma
	430	// TOF ele band 3sigma 0<p<1.5GeV
	431	TF1 lowerCut=new TF1("lowerCut", "[0] TMath::Exp([1]*x)", 0, 100);
	432	lowerCut->SetParameters(-2.7,-0.4357);
	433	AddCut(kTPC,AliPID::kElectron,lowerCut,3.);
	434	AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
	435	} else if (def==7) {
61d106d3	436	// wide TPC cut
8df8e382	437	// TOF ele band 3sigma 0<p<1.5GeV
	438	AddCut(kTPC,AliPID::kElectron,10.);
	439	AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
61d106d3	440	} else if (def==8) {
	441	// TOF 5 sigma inclusion if TOFpid available
	442	// this should reduce K,p,Pi to a large extent
	443	AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
	444	} else if (def==9) {
	445	// lower cut TPC: parametrisation by HFE
	446	// upper cut TPC: 3 sigma
	447	// TOF 5 sigma inclusion if TOFpid available
	448	// this should reduce K,p,Pi to a large extent
	449	TF1 lowerCut=new TF1("lowerCut", "[0] TMath::Exp([1]*x)", 0, 100);
164bfb53	450	lowerCut->SetParameters(-2.65,-0.6757);
164bfb53	451	AddCut(kTPC,AliPID::kElectron,lowerCut,4.);
61d106d3	452	AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
ba15fdfb	453
61d106d3	454	} else if (def==10) {
	455	AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
	456	AddCut(kTPC,AliPID::kElectron,3.);
	457	AddCut(kTPC,AliPID::kPion,-3.,3.,0.,0.,kTRUE);
	458	AddCut(kTPC,AliPID::kProton,-3.,3.,0.,0.,kTRUE);
	459
ba15fdfb	460	} else if (def==11) {
	461	// lower cut TPC: parametrisation by HFE
	462	// only use from period d on !!
	463	// upper cut TPC: 3 sigma
	464	// TOF ele band 3sigma 0<p<2.0GeV
	465	TF1 lowerCut=new TF1("lowerCut", "[0] TMath::Exp([1]*x)+[2]", 0, 100);
	466	lowerCut->SetParameters(-3.718,-0.4,0.27);
	467	AddCut(kTPC,AliPID::kElectron,lowerCut,3.);
	468	AddCut(kTOF,AliPID::kElectron,-3,3,0,2.);
	469	} else if (def==12) {
	470	// lower cut TPC: parametrisation by HFE
	471	// only use from period d on !!
	472	// upper cut TPC: 3 sigma
	473	// TOF 5 sigma inclusion if TOFpid available
	474	// this should reduce K,p,Pi to a large extent
	475	TF1 lowerCut=new TF1("lowerCut", "[0] TMath::Exp([1]*x)+[2]", 0, 100);
	476	lowerCut->SetParameters(-3.718,-0.4,0.27);
	477	AddCut(kTPC,AliPID::kElectron,lowerCut,4.);
	478	AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
8df8e382	479	}
ba15fdfb	480
8df8e382	481	}
8df8e382	482
48609e3d	483	//______________________________________________
	484	void AliDielectronPID::SetCorrVal(Double_t run)
	485	{
	486	//
	487	// set correction value for run
	488	//
d327d9cd	489	fgCorr=0.;
	490	fgCorrdEdx=1.;
	491
	492	if (fgFitCorr){
	493	fgCorr=fgFitCorr->Eval(run);
	494	if (run<fgFitCorr->GetX()[0]) fgCorr=fgFitCorr->GetY()[0];
	495	if (run>fgFitCorr->GetX()[fgFitCorr->GetN()-1]) fgCorr=fgFitCorr->GetY()[fgFitCorr->GetN()-1];
	496	}
	497
	498	if (fgdEdxRunCorr){
	499	fgCorrdEdx=fgdEdxRunCorr->Eval(run);
	500	if (run<fgdEdxRunCorr->GetX()[0]) fgCorrdEdx=fgdEdxRunCorr->GetY()[0];
	501	if (run>fgdEdxRunCorr->GetX()[fgFitCorr->GetN()-1]) fgCorrdEdx=fgdEdxRunCorr->GetY()[fgdEdxRunCorr->GetN()-1];
	502	}
48609e3d	503	}
48609e3d	504
5720c765	505	//______________________________________________
	506	Double_t AliDielectronPID::GetEtaCorr(AliVTrack *track)
	507	{
	508	//
	509	// return eta correction
	510	//
	511	if (!fgFunEtaCorr) return 1;
	512	return fgFunEtaCorr->Eval(track->Eta());
	513	}