[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronTRDpidCut.cxx

/**************************************************************************
* Copyright(c) 1998-1999, 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.                  *
**************************************************************************/
//
// Class for TRD PID
// Class further contains TRD specific cuts and QA histograms 
//  
// Authors:
//   Jens Wiechula <Jens.Wiechula@cern.ch>
//   Markus Fasel <M.Fasel@gsi.de>  
// 
#include <TList.h>
#include <TString.h>

#include <AliVParticle.h>
#include <AliESDtrack.h>
#include <AliLog.h>
#include <AliPID.h>

#include "AliDielectronTRDpidCut.h"

ClassImp(AliDielectronTRDpidCut)

const Double_t AliDielectronTRDpidCut::fgkVerySmall = 1e-12;

//___________________________________________________________________
AliDielectronTRDpidCut::AliDielectronTRDpidCut() :
    AliAnalysisCuts()
  , fMinP(1.)
  , fElectronEfficiency(0.91)
  , fPIDMethod(kNN)
  , fRequirePIDbit(0)
{
  //
  // default  constructor
  // 
  memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
  SetUseDefaultParameters();
}

//___________________________________________________________________
AliDielectronTRDpidCut::AliDielectronTRDpidCut(const char* name) :
    AliAnalysisCuts(name,name)
  , fMinP(1.)
  , fElectronEfficiency(0.91)
  , fPIDMethod(kNN)
  , fRequirePIDbit(0)
{
  //
  // default  constructor
  // 
  memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
  SetUseDefaultParameters();
}

//___________________________________________________________________
AliDielectronTRDpidCut::AliDielectronTRDpidCut(const AliDielectronTRDpidCut &ref):
    AliAnalysisCuts(ref.GetName(),ref.GetTitle())
  , fMinP(ref.fMinP)
  , fElectronEfficiency(ref.fElectronEfficiency)
  , fPIDMethod(ref.fPIDMethod)
  , fRequirePIDbit(ref.fRequirePIDbit)
{
  //
  // Copy constructor
  //
  memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
  ref.Copy(*this);
}

//___________________________________________________________________
AliDielectronTRDpidCut &AliDielectronTRDpidCut::operator=(const AliDielectronTRDpidCut &ref){
  //
  // Assignment operator
  //
  if(this != &ref){
    ref.Copy(*this);
  }
  return *this;
}

//___________________________________________________________________
void AliDielectronTRDpidCut::Copy(TObject &ref) const {
  //
  // Performs the copying of the object
  //
  AliDielectronTRDpidCut &target = dynamic_cast<AliDielectronTRDpidCut &>(ref);
  
  Bool_t defaultParameters = UseDefaultParameters();
  target.SetUseDefaultParameters(defaultParameters);
  target.fMinP = fMinP;
  target.fPIDMethod = fPIDMethod;
  target.fRequirePIDbit = fRequirePIDbit;
  target.fElectronEfficiency = fElectronEfficiency;
  memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams);
  AliAnalysisCuts::Copy(ref);
}

//___________________________________________________________________
AliDielectronTRDpidCut::~AliDielectronTRDpidCut(){
  //
  // Destructor
  //
}

//______________________________________________________
Bool_t AliDielectronTRDpidCut::InitializePID(){
  //
  // InitializePID: Load TRD thresholds and create the electron efficiency axis
  // to navigate 
  //
  if(UseDefaultParameters()){
    if(fPIDMethod == kLQ)
      InitParameters1DLQ();
    else
      InitParameters();
  }
  return kTRUE;
}

//______________________________________________________
Bool_t AliDielectronTRDpidCut::IsSelected(TObject *track) {
  //
  // Does PID for TRD alone:
  // PID thresholds based on 90% Electron Efficiency level approximated by a linear 
  // step function
  //
  const AliESDtrack *part = dynamic_cast<const AliESDtrack *>(track);
  if (!part) return kFALSE;
  
  if (fRequirePIDbit==AliDielectronTRDpidCut::kRequire&&!(part->GetStatus()&AliESDtrack::kTRDpid)) return kFALSE;
  if (fRequirePIDbit==AliDielectronTRDpidCut::kIfAvailable&&!(part->GetStatus()&AliESDtrack::kTRDpid)) return kTRUE;
  
  Double_t p = GetP((AliVParticle*)track);
  if(p < fMinP){ 
    AliDebug(1, Form("Track momentum below %f", fMinP));
    if(fRequirePIDbit==AliDielectronTRDpidCut::kRequire) return kFALSE;
	if(fRequirePIDbit==AliDielectronTRDpidCut::kIfAvailable) return kTRUE;
  }

  Double_t electronLike = GetElectronLikelihood((AliVParticle*)track);
  Double_t threshold = GetTRDthresholds(fElectronEfficiency, p);
  AliDebug(1, Form("Threshold: %f\n", threshold));
  if(electronLike > threshold){
    return kTRUE;
  }
  return kFALSE;

}

//___________________________________________________________________
Double_t AliDielectronTRDpidCut::GetTRDthresholds(Double_t electronEff, Double_t p) const {
  //
  // Return momentum dependent and electron efficiency dependent TRD thresholds
  // 
  Double_t params[4];
  GetParameters(electronEff, params);
  Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
  return TMath::Max(TMath::Min(threshold, 0.99), 0.2); // truncate the threshold upperwards to 0.999 and lowerwards to 0.2 and exclude unphysical values
}

//___________________________________________________________________
void AliDielectronTRDpidCut::SetThresholdParameters(Double_t electronEff, Double_t *params){
  //
  // Set threshold parameters for the given bin
  //
  if(electronEff >= 1. || electronEff < 0.7) return;
  Int_t effbin = static_cast<Int_t>((electronEff - 0.7)/0.05); 
  memcpy(&fThreshParams[effbin * 4], params, sizeof(Double_t) * 4); 
  SetUseDefaultParameters(kFALSE);
}

//___________________________________________________________________
void AliDielectronTRDpidCut::InitParameters(){
  //
  // Fill the Parameters into an array
  //

  AliDebug(2, "Loading threshold Parameter");
  // Parameters for 6 Layers
  fThreshParams[0] = -0.001839; // 0.7 electron eff
  fThreshParams[1] = 0.000276;
  fThreshParams[2] = 0.044902; 
  fThreshParams[3] = 1.726751;
  fThreshParams[4] = -0.002405; // 0.75 electron eff
  fThreshParams[5] = 0.000372;
  fThreshParams[6] = 0.061775;
  fThreshParams[7] = 1.739371;
  fThreshParams[8] = -0.003178; // 0.8 electron eff
  fThreshParams[9] = 0.000521;
  fThreshParams[10] = 0.087585;
  fThreshParams[11] = 1.749154;
  fThreshParams[12] = -0.004058; // 0.85 electron eff
  fThreshParams[13] = 0.000748;
  fThreshParams[14] = 0.129583;
  fThreshParams[15] = 1.782323;
  fThreshParams[16] = -0.004967; // 0.9 electron eff
  fThreshParams[17] = 0.001216;
  fThreshParams[18] = 0.210128;
  fThreshParams[19] = 1.807665;
  fThreshParams[20] = -0.000996; // 0.95 electron eff
  fThreshParams[21] = 0.002627;
  fThreshParams[22] = 0.409099;
  fThreshParams[23] = 1.787076;
}

//___________________________________________________________________
void AliDielectronTRDpidCut::InitParameters1DLQ(){
  //
  // Init Parameters for 1DLQ PID (M. Fasel, Sept. 6th, 2010)
  //

  // Parameters for 6 Layers
  AliDebug(2, Form("Loading threshold parameter for Method 1DLQ"));
  fThreshParams[0] = -0.02241; // 0.7 electron eff
  fThreshParams[1] = 0.05043;
  fThreshParams[2] = 0.7925; 
  fThreshParams[3] = 2.625;
  fThreshParams[4] = 0.07438; // 0.75 electron eff
  fThreshParams[5] = 0.05158;
  fThreshParams[6] = 2.864;
  fThreshParams[7] = 4.356;
  fThreshParams[8] = 0.1977; // 0.8 electron eff
  fThreshParams[9] = 0.05956;
  fThreshParams[10] = 2.853;
  fThreshParams[11] = 3.713;
  fThreshParams[12] = 0.5206; // 0.85 electron eff
  fThreshParams[13] = 0.03077;
  fThreshParams[14] = 2.966;
  fThreshParams[15] = 4.07;
  fThreshParams[16] = 0.8808; // 0.9 electron eff
  fThreshParams[17] = 0.002092;
  fThreshParams[18] = 1.17;
  fThreshParams[19] = 4.506;
  fThreshParams[20] = 1.; // 0.95 electron eff
  fThreshParams[21] = 0.;
  fThreshParams[22] = 0.;
  fThreshParams[23] = 0.;

}

//___________________________________________________________________
void AliDielectronTRDpidCut::RenormalizeElPi(const Double_t *likein, Double_t *likeout) const {
  //
  // Renormalize likelihoods for electrons and pions neglecting the 
  // likelihoods for protons, kaons and muons
  //
  memset(likeout, 0, sizeof(Double_t) * AliPID::kSPECIES);
  Double_t norm = likein[AliPID::kElectron] + likein[AliPID::kPion];
  if(norm == 0.) norm = 1.;   // Safety
  likeout[AliPID::kElectron] = likein[AliPID::kElectron] / norm;
  likeout[AliPID::kPion] = likein[AliPID::kPion] / norm;
}

//___________________________________________________________________
void AliDielectronTRDpidCut::GetParameters(Double_t electronEff, Double_t *parameters) const {
  //
  // return parameter set for the given efficiency bin
  //
  Int_t effbin = static_cast<Int_t>((electronEff - 0.7)/0.05);
  memcpy(parameters, fThreshParams + effbin * 4, sizeof(Double_t) * 4);
}

//___________________________________________________________________
Double_t AliDielectronTRDpidCut::GetElectronLikelihood(const AliVParticle *track) const {
  //
  // Get TRD likelihoods for ESD respectively AOD tracks
  //
  Double_t pidProbs[AliPID::kSPECIES]; memset(pidProbs, 0, sizeof(Double_t) * AliPID::kSPECIES);
  const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
  if(esdtrack) esdtrack->GetTRDpid(pidProbs);
  if(!IsRenormalizeElPi()) return pidProbs[AliPID::kElectron];
  Double_t probsNew[AliPID::kSPECIES];
  RenormalizeElPi(pidProbs, probsNew);
  return probsNew[AliPID::kElectron];
}

//___________________________________________________________________
Double_t AliDielectronTRDpidCut::GetP(const AliVParticle *track) const {
  //
  // Get the Momentum in the TRD
  //
  Double_t p = 0.;
  const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
  if(esdtrack) p = esdtrack->GetOuterParam() ? esdtrack->GetOuterParam()->P() : esdtrack->P();
  return p;
}

//___________________________________________________________________
Double_t AliDielectronTRDpidCut::GetChargeLayer(const AliVParticle *track, UInt_t layer) const {
  //
  // Get the Charge in a single TRD layer
  //
  if(layer >= 6) return 0.;
  Double_t charge = 0.;
  const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
  if(esdtrack)
    for(Int_t islice = 0; islice < esdtrack->GetNumberOfTRDslices(); islice++)
      charge += esdtrack->GetTRDslice(static_cast<UInt_t>(layer), islice);
  return charge;
}

//___________________________________________________________________
void AliDielectronTRDpidCut::GetTRDmomenta(const AliVParticle *track, Double_t *mom) const {
  //
  // Fill Array with momentum information at the TRD tracklet
  //
  const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
  if(esdtrack)
    for(Int_t itl = 0; itl < 6; itl++)
      mom[itl] = esdtrack->GetTRDmomentum(itl);
}

//___________________________________________________________________
Double_t AliDielectronTRDpidCut::GetTRDSignalV1(const AliESDtrack *track, Float_t truncation) const {
  //
  // Calculation of the TRD Signal via truncated mean
  // Method 1: Take all Slices available
  // cut out 0s
  // Order them in increasing order
  // Cut out the upper third
  // Calculate mean over the last 2/3 slices
  //
  const Int_t kNSlices = 48;
  const Int_t kSlicePerLayer = 7;
  // Weight the slice to equalize the MPV of the dQ/dl-distribution per slice to the one in the first slice
  // Pions are used as reference for the equalization
  const Double_t kWeightSlice[8] = {1., 2.122, 1.8, 1.635, 1.595, 1.614, 1.16, 7.0};
  AliDebug(3, Form("Number of Tracklets: %d\n", track->GetTRDntrackletsPID()));
  Double_t trdSlices[kNSlices], tmp[kNSlices];
  Int_t indices[48];
  Int_t icnt = 0;
  for(Int_t idet = 0; idet < 6; idet++)
    for(Int_t islice = 0; islice < kSlicePerLayer; islice++){
      AliDebug(2, Form("Chamber[%d], Slice[%d]: TRDSlice = %f", idet, islice, track->GetTRDslice(idet, islice)));
      if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
      trdSlices[icnt++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
    }
  AliDebug(1, Form("Number of Slices: %d\n", icnt));
  if(icnt < 6) return 0.;   // We need at least 6 Slices for the truncated mean
  TMath::Sort(icnt, trdSlices, indices, kFALSE);
  memcpy(tmp, trdSlices, sizeof(Double_t) * icnt);
  for(Int_t ien = 0; ien < icnt; ien++)
    trdSlices[ien] = tmp[indices[ien]];
  Double_t trdSignal = TMath::Mean(static_cast<Int_t>(static_cast<Float_t>(icnt) * truncation), trdSlices);
  Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
  AliDebug(3, Form("PID Meth. 1: p[%f], TRDSignal[%f]", mom, trdSignal));
  return trdSignal;
}

//___________________________________________________________________
Double_t AliDielectronTRDpidCut::GetTRDSignalV2(const AliESDtrack *track, Float_t truncation) const {
  //
  // Calculation of the TRD Signal via truncated mean
  // Method 2: Take only first 5 slices per chamber
  // Order them in increasing order
  // Cut out upper half 
  // Now do mean with the reamining 3 slices per chamber
  //
  const Double_t kWeightSlice[8] = {1., 2.122, 1.8, 1.635, 1.595, 1.614, 1.16, 7.0};
  const Int_t kLayers = 6;
  const Int_t kSlicesLow = 6;
  const Int_t kSlicesHigh = 1;
  Double_t trdSlicesLowTime[kLayers*kSlicesLow], trdSlicesRemaining[kLayers*(kSlicesHigh + kSlicesLow)];
  Int_t indices[kLayers*kSlicesLow];
  Int_t cntLowTime=0, cntRemaining = 0;
  for(Int_t idet = 0; idet < 6; idet++)
    for(Int_t islice = 0; islice < kSlicesLow+kSlicesHigh; islice++){
      if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
      if(islice < kSlicesLow){
        AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
        trdSlicesLowTime[cntLowTime++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
      } else{
        AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
        trdSlicesRemaining[cntRemaining++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
      }
    }
  if(cntLowTime < 4 || cntRemaining < 2) return 0.; // Min. Number of Slices at high time is 2 (matches with 1 layer), for the truncated mean we need at least 4 Slices
  TMath::Sort(cntLowTime, trdSlicesLowTime, indices, kFALSE);
  // Fill the second array with the lower half of the first time bins
  for(Int_t ien = 0; ien < static_cast<Int_t>(static_cast<Float_t>(cntLowTime) * truncation); ien++)
    trdSlicesRemaining[cntRemaining++] = trdSlicesLowTime[indices[ien]];
  Double_t trdSignal = TMath::Mean(cntRemaining, trdSlicesRemaining);
  Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
  AliDebug(3, Form("PID Meth. 2: p[%f], TRDSignal[%f]", mom, trdSignal));
  return trdSignal;
}
Commit	Line	Data
ba15fdfb	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, 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	// Class for TRD PID
	17	// Class further contains TRD specific cuts and QA histograms
	18	//
	19	// Authors:
	20	// Jens Wiechula <Jens.Wiechula@cern.ch>
	21	// Markus Fasel <M.Fasel@gsi.de>
	22	//
	23	#include <TList.h>
	24	#include <TString.h>
	25
	26	#include <AliVParticle.h>
	27	#include <AliESDtrack.h>
	28	#include <AliLog.h>
	29	#include <AliPID.h>
	30
	31	#include "AliDielectronTRDpidCut.h"
	32
	33	ClassImp(AliDielectronTRDpidCut)
	34
	35	const Double_t AliDielectronTRDpidCut::fgkVerySmall = 1e-12;
	36
	37	//___________________________________________________________________
	38	AliDielectronTRDpidCut::AliDielectronTRDpidCut() :
	39	AliAnalysisCuts()
	40	, fMinP(1.)
	41	, fElectronEfficiency(0.91)
	42	, fPIDMethod(kNN)
	43	, fRequirePIDbit(0)
	44	{
	45	//
	46	// default constructor
	47	//
	48	memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
	49	SetUseDefaultParameters();
	50	}
	51
	52	//___________________________________________________________________
	53	AliDielectronTRDpidCut::AliDielectronTRDpidCut(const char* name) :
	54	AliAnalysisCuts(name,name)
	55	, fMinP(1.)
	56	, fElectronEfficiency(0.91)
	57	, fPIDMethod(kNN)
	58	, fRequirePIDbit(0)
	59	{
	60	//
	61	// default constructor
	62	//
	63	memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
	64	SetUseDefaultParameters();
65	}
66
67	//___________________________________________________________________
68	AliDielectronTRDpidCut::AliDielectronTRDpidCut(const AliDielectronTRDpidCut &ref):
69	AliAnalysisCuts(ref.GetName(),ref.GetTitle())
70	, fMinP(ref.fMinP)
71	, fElectronEfficiency(ref.fElectronEfficiency)
72	, fPIDMethod(ref.fPIDMethod)
73	, fRequirePIDbit(ref.fRequirePIDbit)
74	{
75	//
76	// Copy constructor
77	//
78	memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
79	ref.Copy(*this);
80	}
81
82	//___________________________________________________________________
83	AliDielectronTRDpidCut &AliDielectronTRDpidCut::operator=(const AliDielectronTRDpidCut &ref){
84	//
85	// Assignment operator
86	//
87	if(this != &ref){
88	ref.Copy(*this);
89	}
90	return *this;
91	}
92
93	//___________________________________________________________________
94	void AliDielectronTRDpidCut::Copy(TObject &ref) const {
95	//
96	// Performs the copying of the object
97	//
98	AliDielectronTRDpidCut &target = dynamic_cast<AliDielectronTRDpidCut &>(ref);
99
100	Bool_t defaultParameters = UseDefaultParameters();
101	target.SetUseDefaultParameters(defaultParameters);
102	target.fMinP = fMinP;
103	target.fPIDMethod = fPIDMethod;
104	target.fRequirePIDbit = fRequirePIDbit;
105	target.fElectronEfficiency = fElectronEfficiency;
106	memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams);
107	AliAnalysisCuts::Copy(ref);
108	}
109
110	//___________________________________________________________________
111	AliDielectronTRDpidCut::~AliDielectronTRDpidCut(){
112	//
113	// Destructor
114	//
115	}
116
117	//______________________________________________________
118	Bool_t AliDielectronTRDpidCut::InitializePID(){
119	//
120	// InitializePID: Load TRD thresholds and create the electron efficiency axis
121	// to navigate
122	//
123	if(UseDefaultParameters()){
124	if(fPIDMethod == kLQ)
125	InitParameters1DLQ();
126	else
127	InitParameters();
128	}
129	return kTRUE;
130	}
131
132	//______________________________________________________
133	Bool_t AliDielectronTRDpidCut::IsSelected(TObject *track) {
134	//
135	// Does PID for TRD alone:
136	// PID thresholds based on 90% Electron Efficiency level approximated by a linear
137	// step function
138	//
139	const AliESDtrack part = dynamic_cast<const AliESDtrack >(track);
5cc8c825	140	if (!part) return kFALSE;
ba15fdfb	141
	142	if (fRequirePIDbit==AliDielectronTRDpidCut::kRequire&&!(part->GetStatus()&AliESDtrack::kTRDpid)) return kFALSE;
	143	if (fRequirePIDbit==AliDielectronTRDpidCut::kIfAvailable&&!(part->GetStatus()&AliESDtrack::kTRDpid)) return kTRUE;
	144
	145	Double_t p = GetP((AliVParticle*)track);
	146	if(p < fMinP){
	147	AliDebug(1, Form("Track momentum below %f", fMinP));
	148	if(fRequirePIDbit==AliDielectronTRDpidCut::kRequire) return kFALSE;
	149	if(fRequirePIDbit==AliDielectronTRDpidCut::kIfAvailable) return kTRUE;
	150	}
	151
	152	Double_t electronLike = GetElectronLikelihood((AliVParticle*)track);
	153	Double_t threshold = GetTRDthresholds(fElectronEfficiency, p);
	154	AliDebug(1, Form("Threshold: %f\n", threshold));
	155	if(electronLike > threshold){
	156	return kTRUE;
	157	}
	158	return kFALSE;
	159
	160	}
	161
	162	//___________________________________________________________________
	163	Double_t AliDielectronTRDpidCut::GetTRDthresholds(Double_t electronEff, Double_t p) const {
	164	//
	165	// Return momentum dependent and electron efficiency dependent TRD thresholds
	166	//
	167	Double_t params[4];
	168	GetParameters(electronEff, params);
	169	Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
	170	return TMath::Max(TMath::Min(threshold, 0.99), 0.2); // truncate the threshold upperwards to 0.999 and lowerwards to 0.2 and exclude unphysical values
	171	}
	172
	173	//___________________________________________________________________
	174	void AliDielectronTRDpidCut::SetThresholdParameters(Double_t electronEff, Double_t *params){
	175	//
	176	// Set threshold parameters for the given bin
	177	//
	178	if(electronEff >= 1. \|\| electronEff < 0.7) return;
	179	Int_t effbin = static_cast<Int_t>((electronEff - 0.7)/0.05);
	180	memcpy(&fThreshParams[effbin * 4], params, sizeof(Double_t) * 4);
	181	SetUseDefaultParameters(kFALSE);
	182	}
	183
	184	//___________________________________________________________________
	185	void AliDielectronTRDpidCut::InitParameters(){
	186	//
	187	// Fill the Parameters into an array
	188	//
	189
	190	AliDebug(2, "Loading threshold Parameter");
	191	// Parameters for 6 Layers
	192	fThreshParams[0] = -0.001839; // 0.7 electron eff
	193	fThreshParams[1] = 0.000276;
	194	fThreshParams[2] = 0.044902;
	195	fThreshParams[3] = 1.726751;
	196	fThreshParams[4] = -0.002405; // 0.75 electron eff
	197	fThreshParams[5] = 0.000372;
	198	fThreshParams[6] = 0.061775;
	199	fThreshParams[7] = 1.739371;
	200	fThreshParams[8] = -0.003178; // 0.8 electron eff
	201	fThreshParams[9] = 0.000521;
	202	fThreshParams[10] = 0.087585;
	203	fThreshParams[11] = 1.749154;
	204	fThreshParams[12] = -0.004058; // 0.85 electron eff
205	fThreshParams[13] = 0.000748;
206	fThreshParams[14] = 0.129583;
207	fThreshParams[15] = 1.782323;
208	fThreshParams[16] = -0.004967; // 0.9 electron eff
209	fThreshParams[17] = 0.001216;
210	fThreshParams[18] = 0.210128;
211	fThreshParams[19] = 1.807665;
212	fThreshParams[20] = -0.000996; // 0.95 electron eff
213	fThreshParams[21] = 0.002627;
214	fThreshParams[22] = 0.409099;
215	fThreshParams[23] = 1.787076;
216	}
217
218	//___________________________________________________________________
219	void AliDielectronTRDpidCut::InitParameters1DLQ(){
220	//
221	// Init Parameters for 1DLQ PID (M. Fasel, Sept. 6th, 2010)
222	//
223
224	// Parameters for 6 Layers
225	AliDebug(2, Form("Loading threshold parameter for Method 1DLQ"));
226	fThreshParams[0] = -0.02241; // 0.7 electron eff
227	fThreshParams[1] = 0.05043;
228	fThreshParams[2] = 0.7925;
229	fThreshParams[3] = 2.625;
230	fThreshParams[4] = 0.07438; // 0.75 electron eff
231	fThreshParams[5] = 0.05158;
232	fThreshParams[6] = 2.864;
233	fThreshParams[7] = 4.356;
234	fThreshParams[8] = 0.1977; // 0.8 electron eff
235	fThreshParams[9] = 0.05956;
236	fThreshParams[10] = 2.853;
237	fThreshParams[11] = 3.713;
238	fThreshParams[12] = 0.5206; // 0.85 electron eff
239	fThreshParams[13] = 0.03077;
240	fThreshParams[14] = 2.966;
241	fThreshParams[15] = 4.07;
242	fThreshParams[16] = 0.8808; // 0.9 electron eff
243	fThreshParams[17] = 0.002092;
244	fThreshParams[18] = 1.17;
245	fThreshParams[19] = 4.506;
246	fThreshParams[20] = 1.; // 0.95 electron eff
247	fThreshParams[21] = 0.;
248	fThreshParams[22] = 0.;
249	fThreshParams[23] = 0.;
250
251	}
252
253	//___________________________________________________________________
254	void AliDielectronTRDpidCut::RenormalizeElPi(const Double_t likein, Double_t likeout) const {
255	//
256	// Renormalize likelihoods for electrons and pions neglecting the
257	// likelihoods for protons, kaons and muons
258	//
259	memset(likeout, 0, sizeof(Double_t) * AliPID::kSPECIES);
260	Double_t norm = likein[AliPID::kElectron] + likein[AliPID::kPion];
261	if(norm == 0.) norm = 1.; // Safety
262	likeout[AliPID::kElectron] = likein[AliPID::kElectron] / norm;
263	likeout[AliPID::kPion] = likein[AliPID::kPion] / norm;
264	}
265
266	//___________________________________________________________________
267	void AliDielectronTRDpidCut::GetParameters(Double_t electronEff, Double_t *parameters) const {
268	//
269	// return parameter set for the given efficiency bin
270	//
271	Int_t effbin = static_cast<Int_t>((electronEff - 0.7)/0.05);
272	memcpy(parameters, fThreshParams + effbin * 4, sizeof(Double_t) * 4);
273	}
274
275	//___________________________________________________________________
276	Double_t AliDielectronTRDpidCut::GetElectronLikelihood(const AliVParticle *track) const {
277	//
278	// Get TRD likelihoods for ESD respectively AOD tracks
279	//
280	Double_t pidProbs[AliPID::kSPECIES]; memset(pidProbs, 0, sizeof(Double_t) * AliPID::kSPECIES);
281	const AliESDtrack esdtrack = dynamic_cast<const AliESDtrack >(track);
282	if(esdtrack) esdtrack->GetTRDpid(pidProbs);
283	if(!IsRenormalizeElPi()) return pidProbs[AliPID::kElectron];
284	Double_t probsNew[AliPID::kSPECIES];
285	RenormalizeElPi(pidProbs, probsNew);
286	return probsNew[AliPID::kElectron];
287	}
288
289	//___________________________________________________________________
290	Double_t AliDielectronTRDpidCut::GetP(const AliVParticle *track) const {
291	//
292	// Get the Momentum in the TRD
293	//
294	Double_t p = 0.;
295	const AliESDtrack esdtrack = dynamic_cast<const AliESDtrack >(track);
296	if(esdtrack) p = esdtrack->GetOuterParam() ? esdtrack->GetOuterParam()->P() : esdtrack->P();
297	return p;
298	}
299
300	//___________________________________________________________________
301	Double_t AliDielectronTRDpidCut::GetChargeLayer(const AliVParticle *track, UInt_t layer) const {
302	//
303	// Get the Charge in a single TRD layer
304	//
305	if(layer >= 6) return 0.;
306	Double_t charge = 0.;
307	const AliESDtrack esdtrack = dynamic_cast<const AliESDtrack >(track);
308	if(esdtrack)
309	for(Int_t islice = 0; islice < esdtrack->GetNumberOfTRDslices(); islice++)
310	charge += esdtrack->GetTRDslice(static_cast<UInt_t>(layer), islice);
311	return charge;
312	}
313
314	//___________________________________________________________________
315	void AliDielectronTRDpidCut::GetTRDmomenta(const AliVParticle track, Double_t mom) const {
316	//
317	// Fill Array with momentum information at the TRD tracklet
318	//
319	const AliESDtrack esdtrack = dynamic_cast<const AliESDtrack >(track);
320	if(esdtrack)
321	for(Int_t itl = 0; itl < 6; itl++)
322	mom[itl] = esdtrack->GetTRDmomentum(itl);
323	}
324
325	//___________________________________________________________________
326	Double_t AliDielectronTRDpidCut::GetTRDSignalV1(const AliESDtrack *track, Float_t truncation) const {
327	//
328	// Calculation of the TRD Signal via truncated mean
329	// Method 1: Take all Slices available
330	// cut out 0s
331	// Order them in increasing order
332	// Cut out the upper third
333	// Calculate mean over the last 2/3 slices
334	//
335	const Int_t kNSlices = 48;
336	const Int_t kSlicePerLayer = 7;
337	// Weight the slice to equalize the MPV of the dQ/dl-distribution per slice to the one in the first slice
338	// Pions are used as reference for the equalization
339	const Double_t kWeightSlice[8] = {1., 2.122, 1.8, 1.635, 1.595, 1.614, 1.16, 7.0};
340	AliDebug(3, Form("Number of Tracklets: %d\n", track->GetTRDntrackletsPID()));
341	Double_t trdSlices[kNSlices], tmp[kNSlices];
342	Int_t indices[48];
343	Int_t icnt = 0;
344	for(Int_t idet = 0; idet < 6; idet++)
345	for(Int_t islice = 0; islice < kSlicePerLayer; islice++){
346	AliDebug(2, Form("Chamber[%d], Slice[%d]: TRDSlice = %f", idet, islice, track->GetTRDslice(idet, islice)));
347	if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
348	trdSlices[icnt++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
349	}
350	AliDebug(1, Form("Number of Slices: %d\n", icnt));
351	if(icnt < 6) return 0.; // We need at least 6 Slices for the truncated mean
352	TMath::Sort(icnt, trdSlices, indices, kFALSE);
353	memcpy(tmp, trdSlices, sizeof(Double_t) * icnt);
354	for(Int_t ien = 0; ien < icnt; ien++)
355	trdSlices[ien] = tmp[indices[ien]];
356	Double_t trdSignal = TMath::Mean(static_cast<Int_t>(static_cast<Float_t>(icnt) * truncation), trdSlices);
357	Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
358	AliDebug(3, Form("PID Meth. 1: p[%f], TRDSignal[%f]", mom, trdSignal));
359	return trdSignal;
360	}
361
362	//___________________________________________________________________
363	Double_t AliDielectronTRDpidCut::GetTRDSignalV2(const AliESDtrack *track, Float_t truncation) const {
364	//
365	// Calculation of the TRD Signal via truncated mean
366	// Method 2: Take only first 5 slices per chamber
367	// Order them in increasing order
368	// Cut out upper half
369	// Now do mean with the reamining 3 slices per chamber
370	//
371	const Double_t kWeightSlice[8] = {1., 2.122, 1.8, 1.635, 1.595, 1.614, 1.16, 7.0};
372	const Int_t kLayers = 6;
373	const Int_t kSlicesLow = 6;
374	const Int_t kSlicesHigh = 1;
375	Double_t trdSlicesLowTime[kLayerskSlicesLow], trdSlicesRemaining[kLayers(kSlicesHigh + kSlicesLow)];
376	Int_t indices[kLayers*kSlicesLow];
377	Int_t cntLowTime=0, cntRemaining = 0;
378	for(Int_t idet = 0; idet < 6; idet++)
379	for(Int_t islice = 0; islice < kSlicesLow+kSlicesHigh; islice++){
380	if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
381	if(islice < kSlicesLow){
382	AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
383	trdSlicesLowTime[cntLowTime++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
384	} else{
385	AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
386	trdSlicesRemaining[cntRemaining++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
387	}
388	}
389	if(cntLowTime < 4 \|\| cntRemaining < 2) return 0.; // Min. Number of Slices at high time is 2 (matches with 1 layer), for the truncated mean we need at least 4 Slices
390	TMath::Sort(cntLowTime, trdSlicesLowTime, indices, kFALSE);
391	// Fill the second array with the lower half of the first time bins
392	for(Int_t ien = 0; ien < static_cast<Int_t>(static_cast<Float_t>(cntLowTime) * truncation); ien++)
393	trdSlicesRemaining[cntRemaining++] = trdSlicesLowTime[indices[ien]];
394	Double_t trdSignal = TMath::Mean(cntRemaining, trdSlicesRemaining);
395	Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
396	AliDebug(3, Form("PID Meth. 2: p[%f], TRDSignal[%f]", mom, trdSignal));
397	return trdSignal;
398	}