[u/mrichter/AliRoot.git] / PWG3 / 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 <AliLog.h>
#include <AliESDtrack.h>

#include "AliDielectronVarManager.h"

#include "AliDielectronPID.h"

ClassImp(AliDielectronPID)

TGraph *AliDielectronPID::fgFitCorr=0x0;
Double_t AliDielectronPID::fgCorr=0.0;

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

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

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

//______________________________________________
void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, Double_t nSigmaUp/*=-99999.*/,
                              Double_t pMin/*=0*/, Double_t pMax/*=0*/, Bool_t exclude/*=kFALSE*/,
                              UInt_t pidBitType/*=AliDielectronPID::kRequire*/)
{
  //
  // Add a pid nsigma cut
  // use response of detector 'det' in the momentum range ('pMin') to ['pMax']
  // use a sigma band betwee '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));
  }
  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;
  fPmin[fNcuts]=pMin;
  fPmax[fNcuts]=pMax;
  fExclude[fNcuts]=exclude;
  fRequirePIDbit[fNcuts]=pidBitType;
  ++fNcuts;
  
}

//______________________________________________
void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, TF1 * const funUp,
                              Double_t pMin/*=0*/, Double_t pMax/*=0*/, Bool_t exclude/*=kFALSE*/,
                              UInt_t pidBitType/*=AliDielectronPID::kRequire*/)
{
  //
  // 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.,pMin,pMax,exclude,pidBitType);
}

//______________________________________________
void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, TF1 * const funLow, Double_t nSigmaUp,
                              Double_t pMin/*=0*/, Double_t pMax/*=0*/, Bool_t exclude/*=kFALSE*/,
                              UInt_t pidBitType/*=AliDielectronPID::kRequire*/)
{
  //
  // 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,pMin,pMax,exclude,pidBitType);
}

//______________________________________________
void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, TF1 * const funLow, TF1 * const funUp,
                              Double_t pMin/*=0*/, Double_t pMax/*=0*/, Bool_t exclude/*=kFALSE*/,
                              UInt_t pidBitType/*=AliDielectronPID::kRequire*/)
{
  //
  // 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.,pMin,pMax,exclude,pidBitType);
}

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

  //loop over all cuts
  AliVTrack *part=static_cast<AliVTrack*>(track);
  //TODO: Which momentum to use?
  //      Different momenta for different detectors?
  Double_t mom=part->P();
  
  Bool_t selected=kFALSE;
  fESDpid=AliDielectronVarManager::GetESDpid();
  
  for (UChar_t icut=0; icut<fNcuts; ++icut){
    Double_t pMin=fPmin[icut];
    Double_t pMax=fPmax[icut];

    // test momentum range. In case pMin==pMax use all momenta
    if ( (TMath::Abs(pMin-pMax)>1e-20) && (mom<=pMin || mom>pMax) ) continue;

    // 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);

    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 kTOF:
      selected = IsSelectedTOF(part,icut);
      break;
    }
    if (!selected) return kFALSE;
  }

  return selected;
}

//______________________________________________
Bool_t AliDielectronPID::IsSelectedITS(AliVTrack * const part, Int_t icut) const
{
  //
  // ITS part of the PID check
  // Don't accept the track if there was no pid bit set
  //
  Float_t numberOfSigmas=-1000.;
  
  if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliESDtrack::kITSpid)) return kFALSE;
  if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliESDtrack::kITSpid)) return kTRUE;
  
  if (part->IsA()==AliESDtrack::Class()){
    // ESD case in case the PID bit is not set, don't use this track!
    AliESDtrack *track=static_cast<AliESDtrack*>(part);
    numberOfSigmas=fESDpid->NumberOfSigmasITS(track, fPartType[icut]);
  }else{
    // AOD case
    // FIXME: Is there a place to check whether the PID is was set in ESD???
    AliAODTrack *track=static_cast<AliAODTrack*>(part);
    numberOfSigmas=NumberOfSigmasITS(track, fPartType[icut]);
  }
  Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
  return selected;
}

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

  if (part->IsA()==AliESDtrack::Class()){
    // ESD case in case the PID bit is not set, don't use this track!
    AliESDtrack *track=static_cast<AliESDtrack*>(part);
    numberOfSigmas=fESDpid->NumberOfSigmasTPC(track, fPartType[icut]);
  }else{
    // AOD case
    // FIXME: Is there a place to check whether the PID is was set in ESD???
    AliAODTrack *track=static_cast<AliAODTrack*>(part);
    numberOfSigmas=NumberOfSigmasTPC(track, fPartType[icut]);
  }
//   if (fPartType[icut]==AliPID::kElectron){
//     numberOfSigmas-=fgCorr;
//   }
  Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
  return selected;
}

//______________________________________________
Bool_t AliDielectronPID::IsSelectedTRD(AliVTrack * const /*part*/, Int_t /*icut*/) const
{
  //   
  // TRD part of the pid check
  //
  return kFALSE;
}

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

  if (part->IsA()==AliESDtrack::Class()){
    // ESD case in case the PID bit is not set, don't use this track!
    AliESDtrack *track=static_cast<AliESDtrack*>(part);
    numberOfSigmas=fESDpid->NumberOfSigmasTOF(track, fPartType[icut], fESDpid->GetTOFResponse().GetTimeZero());
  }else{
    // AOD case
    // FIXME: Is there a place to check whether the PID is was set in ESD???
    AliAODTrack *track=static_cast<AliAODTrack*>(part);
    numberOfSigmas=NumberOfSigmasTOF(track, fPartType[icut]);
  }
  
  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) {
    // 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==4) {
    // 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);
    
  }
}

//______________________________________________
void AliDielectronPID::SetCorrVal(Double_t run)
{
  //
  // set correction value for run
  //
  fgCorr=0;
  if (!fgFitCorr) return;
  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];
//   printf("Corr: %f\n",fgCorr);
}
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>
8df8e382	33	#include <AliLog.h>
	34	#include <AliESDtrack.h>
	35
	36	#include "AliDielectronVarManager.h"
	37
	38	#include "AliDielectronPID.h"
	39
	40	ClassImp(AliDielectronPID)
	41
48609e3d	42	TGraph *AliDielectronPID::fgFitCorr=0x0;
	43	Double_t AliDielectronPID::fgCorr=0.0;
	44
8df8e382	45	AliDielectronPID::AliDielectronPID() :
	46	AliAnalysisCuts(),
	47	fNcuts(0),
8df8e382	48	fESDpid(0x0)
	49	{
	50	//
	51	// Default Constructor
	52	//
	53	for (Int_t icut=0; icut<kNmaxPID; ++icut){
	54	fDetType[icut]=kTPC;
	55	fPartType[icut]=AliPID::kPion;
	56	fNsigmaLow[icut]=0;
	57	fNsigmaUp[icut]=0;
	58	fPmin[icut]=0;
	59	fPmax[icut]=0;
	60	fExclude[icut]=kFALSE;
	61	fFunUpperCut[icut]=0x0;
	62	fFunLowerCut[icut]=0x0;
61d106d3	63	fRequirePIDbit[icut]=0;
8df8e382	64	}
	65	}
	66
	67	//______________________________________________
	68	AliDielectronPID::AliDielectronPID(const char* name, const char* title) :
	69	AliAnalysisCuts(name, title),
	70	fNcuts(0),
8df8e382	71	fESDpid(0x0)
	72	{
	73	//
	74	// Named Constructor
	75	//
	76	for (Int_t icut=0; icut<kNmaxPID; ++icut){
	77	fDetType[icut]=kTPC;
	78	fPartType[icut]=AliPID::kPion;
	79	fNsigmaLow[icut]=0;
	80	fNsigmaUp[icut]=0;
	81	fPmin[icut]=0;
	82	fPmax[icut]=0;
	83	fExclude[icut]=kFALSE;
	84	fFunUpperCut[icut]=0x0;
	85	fFunLowerCut[icut]=0x0;
61d106d3	86	fRequirePIDbit[icut]=0;
8df8e382	87	}
	88	}
	89
	90	//______________________________________________
	91	AliDielectronPID::~AliDielectronPID()
	92	{
	93	//
	94	// Default Destructor
	95	//
	96	}
	97
	98	//______________________________________________
	99	void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, Double_t nSigmaUp/=-99999./,
61d106d3	100	Double_t pMin/=0/, Double_t pMax/=0/, Bool_t exclude/=kFALSE/,
61d106d3	101	UInt_t pidBitType/=AliDielectronPID::kRequire/)
8df8e382	102	{
	103	//
	104	// Add a pid nsigma cut
	105	// use response of detector 'det' in the momentum range ('pMin') to ['pMax']
	106	// use a sigma band betwee 'nSigmaLow' and 'nSigmaUp'
	107	// if nSigmaUp==-99999. then nSigmaLow will be uesd as a symmetric band +- nSigmaLow
	108	// specify whether to 'exclude' the given band
	109	//
	110
	111	if (fNcuts==kNmaxPID){
	112	AliError(Form("only %d pid cut ranges allowed",kNmaxPID));
	113	}
	114	if (TMath::Abs(nSigmaUp+99999.)<1e-20){
	115	nSigmaUp=TMath::Abs(nSigmaLow);
	116	nSigmaLow=-1.*nSigmaUp;
	117	}
	118	fDetType[fNcuts]=det;
	119	fPartType[fNcuts]=type;
	120	fNsigmaLow[fNcuts]=nSigmaLow;
	121	fNsigmaUp[fNcuts]=nSigmaUp;
	122	fPmin[fNcuts]=pMin;
	123	fPmax[fNcuts]=pMax;
	124	fExclude[fNcuts]=exclude;
61d106d3	125	fRequirePIDbit[fNcuts]=pidBitType;
8df8e382	126	++fNcuts;
	127
	128	}
	129
	130	//______________________________________________
	131	void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, TF1 * const funUp,
61d106d3	132	Double_t pMin/=0/, Double_t pMax/=0/, Bool_t exclude/=kFALSE/,
61d106d3	133	UInt_t pidBitType/=AliDielectronPID::kRequire/)
8df8e382	134	{
	135	//
	136	// cut using a TF1 as upper cut
	137	//
	138	if (funUp==0x0){
	139	AliError("A valid function is required for the upper cut. Not adding the cut!");
	140	return;
	141	}
	142	fFunUpperCut[fNcuts]=funUp;
61d106d3	143	AddCut(det,type,nSigmaLow,0.,pMin,pMax,exclude,pidBitType);
8df8e382	144	}
	145
	146	//______________________________________________
	147	void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, TF1 * const funLow, Double_t nSigmaUp,
61d106d3	148	Double_t pMin/=0/, Double_t pMax/=0/, Bool_t exclude/=kFALSE/,
61d106d3	149	UInt_t pidBitType/=AliDielectronPID::kRequire/)
8df8e382	150	{
	151	//
	152	// cut using a TF1 as lower cut
	153	//
	154	if (funLow==0x0){
	155	AliError("A valid function is required for the lower cut. Not adding the cut!");
	156	return;
	157	}
	158	fFunLowerCut[fNcuts]=funLow;
61d106d3	159	AddCut(det,type,0.,nSigmaUp,pMin,pMax,exclude,pidBitType);
8df8e382	160	}
	161
	162	//______________________________________________
	163	void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, TF1 * const funLow, TF1 * const funUp,
61d106d3	164	Double_t pMin/=0/, Double_t pMax/=0/, Bool_t exclude/=kFALSE/,
61d106d3	165	UInt_t pidBitType/=AliDielectronPID::kRequire/)
8df8e382	166	{
	167	//
	168	// cut using a TF1 as lower and upper cut
	169	//
	170	if ( (funUp==0x0) \|\| (funLow==0x0) ){
	171	AliError("A valid function is required for upper and lower cut. Not adding the cut!");
	172	return;
	173	}
	174	fFunUpperCut[fNcuts]=funUp;
	175	fFunLowerCut[fNcuts]=funLow;
61d106d3	176	AddCut(det,type,0.,0.,pMin,pMax,exclude,pidBitType);
8df8e382	177	}
	178
	179	//______________________________________________
	180	Bool_t AliDielectronPID::IsSelected(TObject* track)
	181	{
	182	//
	183	// perform PID cuts
	184	//
	185
	186	//loop over all cuts
164bfb53	187	AliVTrack part=static_cast<AliVTrack>(track);
8df8e382	188	//TODO: Which momentum to use?
	189	// Different momenta for different detectors?
	190	Double_t mom=part->P();
	191
	192	Bool_t selected=kFALSE;
	193	fESDpid=AliDielectronVarManager::GetESDpid();
	194
	195	for (UChar_t icut=0; icut<fNcuts; ++icut){
	196	Double_t pMin=fPmin[icut];
	197	Double_t pMax=fPmax[icut];
	198
	199	// test momentum range. In case pMin==pMax use all momenta
	200	if ( (TMath::Abs(pMin-pMax)>1e-20) && (mom<=pMin \|\| mom>pMax) ) continue;
	201
	202	// test if we are supposed to use a function for the cut
	203	if (fFunUpperCut[icut]) fNsigmaUp[icut] =fFunUpperCut[icut]->Eval(mom);
	204	if (fFunLowerCut[icut]) fNsigmaLow[icut]=fFunLowerCut[icut]->Eval(mom);
	205
	206	switch (fDetType[icut]){
	207	case kITS:
	208	selected = IsSelectedITS(part,icut);
	209	break;
	210	case kTPC:
	211	selected = IsSelectedTPC(part,icut);
	212	break;
	213	case kTRD:
	214	selected = IsSelectedTRD(part,icut);
	215	break;
	216	case kTOF:
	217	selected = IsSelectedTOF(part,icut);
	218	break;
	219	}
	220	if (!selected) return kFALSE;
	221	}
	222
	223	return selected;
	224	}
	225
	226	//______________________________________________
164bfb53	227	Bool_t AliDielectronPID::IsSelectedITS(AliVTrack * const part, Int_t icut) const
8df8e382	228	{
	229	//
	230	// ITS part of the PID check
	231	// Don't accept the track if there was no pid bit set
	232	//
	233	Float_t numberOfSigmas=-1000.;
	234
164bfb53	235	if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliESDtrack::kITSpid)) return kFALSE;
	236	if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliESDtrack::kITSpid)) return kTRUE;
	237
8df8e382	238	if (part->IsA()==AliESDtrack::Class()){
	239	// ESD case in case the PID bit is not set, don't use this track!
	240	AliESDtrack track=static_cast<AliESDtrack>(part);
8df8e382	241	numberOfSigmas=fESDpid->NumberOfSigmasITS(track, fPartType[icut]);
	242	}else{
	243	// AOD case
	244	// FIXME: Is there a place to check whether the PID is was set in ESD???
	245	AliAODTrack track=static_cast<AliAODTrack>(part);
	246	numberOfSigmas=NumberOfSigmasITS(track, fPartType[icut]);
	247	}
	248	Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
	249	return selected;
	250	}
	251
	252	//______________________________________________
164bfb53	253	Bool_t AliDielectronPID::IsSelectedTPC(AliVTrack * const part, Int_t icut) const
8df8e382	254	{
	255	//
	256	// TPC part of the PID check
	257	// Don't accept the track if there was no pid bit set
	258	//
	259	Float_t numberOfSigmas=-1000.;
	260
48609e3d	261	if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliESDtrack::kTPCpid)) return kFALSE;
	262	if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliESDtrack::kTPCpid)) return kTRUE;
	263
8df8e382	264	if (part->IsA()==AliESDtrack::Class()){
	265	// ESD case in case the PID bit is not set, don't use this track!
	266	AliESDtrack track=static_cast<AliESDtrack>(part);
8df8e382	267	numberOfSigmas=fESDpid->NumberOfSigmasTPC(track, fPartType[icut]);
	268	}else{
	269	// AOD case
	270	// FIXME: Is there a place to check whether the PID is was set in ESD???
	271	AliAODTrack track=static_cast<AliAODTrack>(part);
	272	numberOfSigmas=NumberOfSigmasTPC(track, fPartType[icut]);
	273	}
48609e3d	274	// if (fPartType[icut]==AliPID::kElectron){
	275	// numberOfSigmas-=fgCorr;
	276	// }
8df8e382	277	Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
	278	return selected;
	279	}
	280
	281	//______________________________________________
164bfb53	282	Bool_t AliDielectronPID::IsSelectedTRD(AliVTrack * const /part/, Int_t /icut/) const
8df8e382	283	{
	284	//
	285	// TRD part of the pid check
	286	//
	287	return kFALSE;
	288	}
	289
	290	//______________________________________________
164bfb53	291	Bool_t AliDielectronPID::IsSelectedTOF(AliVTrack * const part, Int_t icut) const
8df8e382	292	{
	293	//
	294	// TOF part of the PID check
	295	// Don't accept the track if there was no pid bit set
	296	//
	297	Float_t numberOfSigmas=-1000.;
	298
48609e3d	299	if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliESDtrack::kTOFpid)) return kFALSE;
	300	if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliESDtrack::kTOFpid)) return kTRUE;
	301
8df8e382	302	if (part->IsA()==AliESDtrack::Class()){
	303	// ESD case in case the PID bit is not set, don't use this track!
	304	AliESDtrack track=static_cast<AliESDtrack>(part);
8df8e382	305	numberOfSigmas=fESDpid->NumberOfSigmasTOF(track, fPartType[icut], fESDpid->GetTOFResponse().GetTimeZero());
	306	}else{
	307	// AOD case
	308	// FIXME: Is there a place to check whether the PID is was set in ESD???
	309	AliAODTrack track=static_cast<AliAODTrack>(part);
	310	numberOfSigmas=NumberOfSigmasTOF(track, fPartType[icut]);
	311	}
	312
	313	Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
	314	return selected;
	315	}
	316
	317	//______________________________________________
	318	void AliDielectronPID::SetDefaults(Int_t def){
	319	//
	320	// initialise default pid strategies
	321	//
	322
	323	if (def==0){
	324	// 2sigma bands TPC:
	325	// - include e
	326	// - exclude mu,K,pi,p
	327	// -complete p range
	328	AddCut(kTPC,AliPID::kElectron,2);
	329	AddCut(kTPC,AliPID::kMuon,-2.,2.,0.,0.,kTRUE);
	330	AddCut(kTPC,AliPID::kPion,-2.,2.,0.,0.,kTRUE);
	331	AddCut(kTPC,AliPID::kKaon,-2.,2.,0.,0.,kTRUE);
	332	AddCut(kTPC,AliPID::kProton,-2.,2.,0.,0.,kTRUE);
	333	} else if (def==1) {
	334	// 2sigma bands TPC:
	335	// - include e 0<p<inf
	336	// - exclude mu,K,pi,p 0<p<2
	337
	338	AddCut(kTPC,AliPID::kElectron,2.);
	339	AddCut(kTPC,AliPID::kMuon,-2.,2.,0.,2.,kTRUE);
	340	AddCut(kTPC,AliPID::kPion,-2.,2.,0.,2.,kTRUE);
	341	AddCut(kTPC,AliPID::kKaon,-2.,2.,0.,2.,kTRUE);
	342	AddCut(kTPC,AliPID::kProton,-2.,2.,0.,2.,kTRUE);
	343
	344	} else if (def==2) {
	345	// include 2sigma e TPC
	346	// 3sigma bands TOF
	347	// - exclude K,P
48609e3d	348	AddCut(kTPC,AliPID::kElectron,-3.,3.);
3505bfad	349	AddCut(kTPC,AliPID::kPion,-3.,3.,0.,0.,kTRUE);
3505bfad	350	AddCut(kTPC,AliPID::kProton,-3.,3.,0.,0.,kTRUE);
8df8e382	351
	352	} else if (def==3) {
	353	// include 2sigma e TPC
	354	// 3sigma bands TOF
	355	// - exclude K 0<p<1
	356	// - exclude P 0<p<2
	357	AddCut(kTPC,AliPID::kElectron,2);
	358	AddCut(kTOF,AliPID::kKaon,-3.,3.,0.,1.,kTRUE);
	359	AddCut(kTOF,AliPID::kProton,-6.,6.,0.,1.,kTRUE);
	360	AddCut(kTOF,AliPID::kProton,-3.,3.,1.,2.,kTRUE);
	361
	362	} else if (def==4) {
	363	// include 2sigma e TPC
	364	// 3sigma bands TOF
	365	// - exclude K 0<p<1
	366	// - exclude P 0<p<2
	367	AddCut(kTPC,AliPID::kElectron,2.);
	368	AddCut(kTOF,AliPID::kKaon,-3.,3.,0.,1.,kTRUE);
	369	AddCut(kTOF,AliPID::kProton,-6.,6.,0.,1.,kTRUE);
	370	AddCut(kTOF,AliPID::kProton,-3.,3.,1.,2.,kTRUE);
8df8e382	371	} else if (def==5) {
	372	AddCut(kTPC,AliPID::kElectron,-0.5,3);
	373	AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
	374	} else if (def==6) {
	375	// lower cut TPC: parametrisation by HFE
	376	// upper cut TPC: 3 sigma
	377	// TOF ele band 3sigma 0<p<1.5GeV
	378	TF1 lowerCut=new TF1("lowerCut", "[0] TMath::Exp([1]*x)", 0, 100);
	379	lowerCut->SetParameters(-2.7,-0.4357);
	380	AddCut(kTPC,AliPID::kElectron,lowerCut,3.);
	381	AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
	382	} else if (def==7) {
61d106d3	383	// wide TPC cut
8df8e382	384	// TOF ele band 3sigma 0<p<1.5GeV
	385	AddCut(kTPC,AliPID::kElectron,10.);
	386	AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
61d106d3	387	} else if (def==8) {
	388	// TOF 5 sigma inclusion if TOFpid available
	389	// this should reduce K,p,Pi to a large extent
	390	AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
	391	} else if (def==9) {
	392	// lower cut TPC: parametrisation by HFE
	393	// upper cut TPC: 3 sigma
	394	// TOF 5 sigma inclusion if TOFpid available
	395	// this should reduce K,p,Pi to a large extent
	396	TF1 lowerCut=new TF1("lowerCut", "[0] TMath::Exp([1]*x)", 0, 100);
164bfb53	397	lowerCut->SetParameters(-2.65,-0.6757);
164bfb53	398	AddCut(kTPC,AliPID::kElectron,lowerCut,4.);
61d106d3	399	AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
	400	} else if (def==10) {
	401	AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
	402	AddCut(kTPC,AliPID::kElectron,3.);
	403	AddCut(kTPC,AliPID::kPion,-3.,3.,0.,0.,kTRUE);
	404	AddCut(kTPC,AliPID::kProton,-3.,3.,0.,0.,kTRUE);
	405
8df8e382	406	}
	407	}
	408
48609e3d	409	//______________________________________________
	410	void AliDielectronPID::SetCorrVal(Double_t run)
	411	{
	412	//
	413	// set correction value for run
	414	//
	415	fgCorr=0;
	416	if (!fgFitCorr) return;
	417	fgCorr=fgFitCorr->Eval(run);
	418	if (run<fgFitCorr->GetX()[0]) fgCorr=fgFitCorr->GetY()[0];
	419	if (run>fgFitCorr->GetX()[fgFitCorr->GetN()-1]) fgCorr=fgFitCorr->GetY()[fgFitCorr->GetN()-1];
	420	// printf("Corr: %f\n",fgCorr);
	421	}
	422