[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEpidTRD.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
// Implements the abstract base class AliHFEpidbase 
// Make PID does the PID decision 
// Class further contains TRD specific cuts and QA histograms 
//  
// Authors: 
//   Markus Fasel <M.Fasel@gsi.de>  
// 
#include <TH1F.h>
#include <TH2F.h>
#include <THnSparse.h>
#include <TList.h>
#include <TString.h>

#include "AliAODpidUtil.h"
#include "AliAODTrack.h"
#include "AliAODMCParticle.h"
#include "AliESDtrack.h"
#include "AliESDpid.h"
#include "AliLog.h"
#include "AliMCParticle.h"
#include "AliPID.h"

#include "AliHFEpidQAmanager.h"
#include "AliHFEpidTRD.h"

ClassImp(AliHFEpidTRD)

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

//___________________________________________________________________
AliHFEpidTRD::AliHFEpidTRD() :
    AliHFEpidBase()
  , fMinP(1.)
  , fElectronEfficiency(0.91)
  , fPIDMethod(kNN)
{
  //
  // default  constructor
  // 
}

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

//___________________________________________________________________
AliHFEpidTRD::AliHFEpidTRD(const AliHFEpidTRD &ref):
    AliHFEpidBase("")
  , fMinP(ref.fMinP)
  , fElectronEfficiency(ref.fElectronEfficiency)
  , fPIDMethod(ref.fPIDMethod)
{
  //
  // Copy constructor
  //
  memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
  ref.Copy(*this);
}

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

//___________________________________________________________________
void AliHFEpidTRD::Copy(TObject &ref) const {
  //
  // Performs the copying of the object
  //
  AliHFEpidTRD &target = dynamic_cast<AliHFEpidTRD &>(ref);

  target.fMinP = fMinP;
  target.fPIDMethod = fPIDMethod;
  target.fElectronEfficiency = fElectronEfficiency;
  memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams);
  AliHFEpidBase::Copy(ref);
}

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

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

//______________________________________________________
Int_t AliHFEpidTRD::IsSelected(AliHFEpidObject *track, AliHFEpidQAmanager *pidqa){
  //
  // Does PID for TRD alone:
  // PID thresholds based on 90% Electron Efficiency level approximated by a linear 
  // step function
  //
  AliDebug(1, "Applying TRD PID");
  if((!fESDpid && track->IsESDanalysis()) || (!fAODpid && track->IsAODanalysis())){ 
    AliDebug(1, "Cannot process track");
    return 0;
  }

  AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis: AliHFEpidObject::kAODanalysis;
  Double_t p = GetP(track->GetRecTrack(), anatype);
  if(p < fMinP){ 
    AliDebug(1, Form("Track momentum below %f", fMinP));
    return 0;
  }

  if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kBeforePID); 
  Double_t electronLike = GetElectronLikelihood(track->GetRecTrack(), anatype);
  Double_t threshold = GetTRDthresholds(fElectronEfficiency, p);
  AliDebug(1, Form("Threshold: %f\n", threshold));
  if(electronLike > threshold){
    if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kAfterPID);
    return 11;
  }
  return 211;

}

//___________________________________________________________________
Double_t AliHFEpidTRD::GetTRDthresholds(Double_t electronEff, Double_t p){ 
  //
  // 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 AliHFEpidTRD::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 AliHFEpidTRD::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 AliHFEpidTRD::GetParameters(Double_t electronEff, Double_t *parameters){
  //
  // 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 AliHFEpidTRD::GetElectronLikelihood(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType){
  //
  // Get TRD likelihoods for ESD respectively AOD tracks
  //
  Double_t pidProbs[AliPID::kSPECIES]; memset(pidProbs, 0, sizeof(Double_t) * AliPID::kSPECIES);
  if(anaType == AliHFEpidObject::kESDanalysis){
    const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
    esdtrack->GetTRDpid(pidProbs);
  } else {
    const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
    fAODpid->MakeTRDPID(const_cast<AliAODTrack *>(aodtrack), pidProbs);
  }
  return pidProbs[AliPID::kElectron];
}

//___________________________________________________________________
Double_t AliHFEpidTRD::GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType){
  //
  // Get the Momentum in the TRD
  //
  Double_t p = 0.;
  if(anaType == AliHFEpidObject::kESDanalysis){
    const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
    p = esdtrack->GetOuterParam() ? esdtrack->GetOuterParam()->P() : esdtrack->P();
  } else {
    const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
    p = aodtrack->P();
  }
  return p;
}

//___________________________________________________________________
Double_t AliHFEpidTRD::GetTRDSignalV1(const AliESDtrack *track){
  //
  // 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;
  AliDebug(3, Form("Number of Tracklets: %d\n", track->GetTRDpidQuality()));
  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 < 8; 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);
    }
  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<Double_t>(icnt) * 2./3.), 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 AliHFEpidTRD::GetTRDSignalV2(const AliESDtrack *track){
  //
  // 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
  //
  Double_t trdSlicesLowTime[30], trdSlicesRemaining[32];
  Int_t indices[30];
  Int_t cntLowTime=0, cntRemaining = 0;
  for(Int_t idet = 0; idet < 6; idet++)
    for(Int_t islice = 0; islice < 8; islice++){
      if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
      if(islice < 5){
        AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
        trdSlicesLowTime[cntLowTime++] = track->GetTRDslice(idet, 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);
      }
    }
  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<Double_t>(cntLowTime) * 0.5); 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
809a4336	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	**************************************************************************/
50685501	15	//
	16	// Class for TRD PID
	17	// Implements the abstract base class AliHFEpidbase
	18	// Make PID does the PID decision
	19	// Class further contains TRD specific cuts and QA histograms
	20	//
	21	// Authors:
	22	// Markus Fasel <M.Fasel@gsi.de>
	23	//
809a4336	24	#include <TH1F.h>
75d81601	25	#include <TH2F.h>
e3fc062d	26	#include <THnSparse.h>
50685501	27	#include <TList.h>
809a4336	28	#include <TString.h>
809a4336	29
3a72645a	30	#include "AliAODpidUtil.h"
722347d8	31	#include "AliAODTrack.h"
722347d8	32	#include "AliAODMCParticle.h"
809a4336	33	#include "AliESDtrack.h"
e3fc062d	34	#include "AliESDpid.h"
75d81601	35	#include "AliLog.h"
722347d8	36	#include "AliMCParticle.h"
809a4336	37	#include "AliPID.h"
809a4336	38
3a72645a	39	#include "AliHFEpidQAmanager.h"
809a4336	40	#include "AliHFEpidTRD.h"
	41
	42	ClassImp(AliHFEpidTRD)
	43
75d81601	44	const Double_t AliHFEpidTRD::fgkVerySmall = 1e-12;
75d81601	45
faee3b18	46	//___________________________________________________________________
	47	AliHFEpidTRD::AliHFEpidTRD() :
	48	AliHFEpidBase()
e3fc062d	49	, fMinP(1.)
67fe7bd0	50	, fElectronEfficiency(0.91)
faee3b18	51	, fPIDMethod(kNN)
faee3b18	52	{
	53	//
	54	// default constructor
	55	//
faee3b18	56	}
faee3b18	57
809a4336	58	//___________________________________________________________________
	59	AliHFEpidTRD::AliHFEpidTRD(const char* name) :
	60	AliHFEpidBase(name)
e3fc062d	61	, fMinP(1.)
67fe7bd0	62	, fElectronEfficiency(0.91)
809a4336	63	, fPIDMethod(kNN)
809a4336	64	{
	65	//
	66	// default constructor
	67	//
dbe3abbe	68	memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
809a4336	69	}
	70
	71	//___________________________________________________________________
	72	AliHFEpidTRD::AliHFEpidTRD(const AliHFEpidTRD &ref):
	73	AliHFEpidBase("")
67fe7bd0	74	, fMinP(ref.fMinP)
	75	, fElectronEfficiency(ref.fElectronEfficiency)
	76	, fPIDMethod(ref.fPIDMethod)
809a4336	77	{
	78	//
	79	// Copy constructor
	80	//
dbe3abbe	81	memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
809a4336	82	ref.Copy(*this);
	83	}
	84
	85	//___________________________________________________________________
	86	AliHFEpidTRD &AliHFEpidTRD::operator=(const AliHFEpidTRD &ref){
	87	//
	88	// Assignment operator
	89	//
	90	if(this != &ref){
	91	ref.Copy(*this);
	92	}
	93	return *this;
	94	}
	95
	96	//___________________________________________________________________
	97	void AliHFEpidTRD::Copy(TObject &ref) const {
	98	//
	99	// Performs the copying of the object
	100	//
	101	AliHFEpidTRD &target = dynamic_cast<AliHFEpidTRD &>(ref);
	102
e3fc062d	103	target.fMinP = fMinP;
809a4336	104	target.fPIDMethod = fPIDMethod;
67fe7bd0	105	target.fElectronEfficiency = fElectronEfficiency;
dbe3abbe	106	memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams);
809a4336	107	AliHFEpidBase::Copy(ref);
	108	}
	109
	110	//___________________________________________________________________
	111	AliHFEpidTRD::~AliHFEpidTRD(){
	112	//
	113	// Destructor
	114	//
809a4336	115	}
	116
	117	//______________________________________________________
	118	Bool_t AliHFEpidTRD::InitializePID(){
	119	//
	120	// InitializePID: Load TRD thresholds and create the electron efficiency axis
	121	// to navigate
	122	//
e3fc062d	123	if(fPIDMethod == kLQ)
	124	InitParameters1DLQ();
	125	else
	126	InitParameters();
809a4336	127	return kTRUE;
	128	}
	129
	130	//______________________________________________________
3a72645a	131	Int_t AliHFEpidTRD::IsSelected(AliHFEpidObject track, AliHFEpidQAmanager pidqa){
809a4336	132	//
	133	// Does PID for TRD alone:
	134	// PID thresholds based on 90% Electron Efficiency level approximated by a linear
	135	// step function
	136	//
3a72645a	137	AliDebug(1, "Applying TRD PID");
	138	if((!fESDpid && track->IsESDanalysis()) \|\| (!fAODpid && track->IsAODanalysis())){
	139	AliDebug(1, "Cannot process track");
	140	return 0;
	141	}
722347d8	142
3a72645a	143	AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis: AliHFEpidObject::kAODanalysis;
	144	Double_t p = GetP(track->GetRecTrack(), anatype);
	145	if(p < fMinP){
	146	AliDebug(1, Form("Track momentum below %f", fMinP));
	147	return 0;
	148	}
809a4336	149
3a72645a	150	if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kBeforePID);
3a72645a	151	Double_t electronLike = GetElectronLikelihood(track->GetRecTrack(), anatype);
67fe7bd0	152	Double_t threshold = GetTRDthresholds(fElectronEfficiency, p);
0792aa82	153	AliDebug(1, Form("Threshold: %f\n", threshold));
3a72645a	154	if(electronLike > threshold){
3a72645a	155	if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kAfterPID);
0792aa82	156	return 11;
	157	}
	158	return 211;
3a72645a	159
809a4336	160	}
	161
	162	//___________________________________________________________________
dbe3abbe	163	Double_t AliHFEpidTRD::GetTRDthresholds(Double_t electronEff, Double_t p){
	164	//
	165	// Return momentum dependent and electron efficiency dependent TRD thresholds
	166	//
	167	Double_t params[4];
	168	GetParameters(electronEff, params);
0792aa82	169	Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
0792aa82	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
809a4336	171	}
	172
	173	//___________________________________________________________________
dbe3abbe	174	void AliHFEpidTRD::InitParameters(){
	175	//
	176	// Fill the Parameters into an array
	177	//
	178
e3fc062d	179	AliDebug(2, "Loading threshold Parameter");
dbe3abbe	180	// Parameters for 6 Layers
	181	fThreshParams[0] = -0.001839; // 0.7 electron eff
	182	fThreshParams[1] = 0.000276;
	183	fThreshParams[2] = 0.044902;
	184	fThreshParams[3] = 1.726751;
	185	fThreshParams[4] = -0.002405; // 0.75 electron eff
	186	fThreshParams[5] = 0.000372;
	187	fThreshParams[6] = 0.061775;
	188	fThreshParams[7] = 1.739371;
	189	fThreshParams[8] = -0.003178; // 0.8 electron eff
	190	fThreshParams[9] = 0.000521;
	191	fThreshParams[10] = 0.087585;
	192	fThreshParams[11] = 1.749154;
	193	fThreshParams[12] = -0.004058; // 0.85 electron eff
	194	fThreshParams[13] = 0.000748;
	195	fThreshParams[14] = 0.129583;
	196	fThreshParams[15] = 1.782323;
	197	fThreshParams[16] = -0.004967; // 0.9 electron eff
	198	fThreshParams[17] = 0.001216;
	199	fThreshParams[18] = 0.210128;
	200	fThreshParams[19] = 1.807665;
	201	fThreshParams[20] = -0.000996; // 0.95 electron eff
	202	fThreshParams[21] = 0.002627;
	203	fThreshParams[22] = 0.409099;
	204	fThreshParams[23] = 1.787076;
809a4336	205	}
809a4336	206
e3fc062d	207	//___________________________________________________________________
	208	void AliHFEpidTRD::InitParameters1DLQ(){
	209	//
	210	// Init Parameters for 1DLQ PID (M. Fasel, Sept. 6th, 2010)
	211	//
	212
	213	// Parameters for 6 Layers
	214	AliDebug(2, Form("Loading threshold parameter for Method 1DLQ"));
	215	fThreshParams[0] = -0.02241; // 0.7 electron eff
	216	fThreshParams[1] = 0.05043;
	217	fThreshParams[2] = 0.7925;
	218	fThreshParams[3] = 2.625;
	219	fThreshParams[4] = 0.07438; // 0.75 electron eff
	220	fThreshParams[5] = 0.05158;
	221	fThreshParams[6] = 2.864;
	222	fThreshParams[7] = 4.356;
	223	fThreshParams[8] = 0.1977; // 0.8 electron eff
	224	fThreshParams[9] = 0.05956;
	225	fThreshParams[10] = 2.853;
	226	fThreshParams[11] = 3.713;
	227	fThreshParams[12] = 0.5206; // 0.85 electron eff
	228	fThreshParams[13] = 0.03077;
	229	fThreshParams[14] = 2.966;
	230	fThreshParams[15] = 4.07;
	231	fThreshParams[16] = 0.8808; // 0.9 electron eff
	232	fThreshParams[17] = 0.002092;
	233	fThreshParams[18] = 1.17;
	234	fThreshParams[19] = 4.506;
	235	fThreshParams[20] = 1.; // 0.95 electron eff
	236	fThreshParams[21] = 0.;
	237	fThreshParams[22] = 0.;
	238	fThreshParams[23] = 0.;
	239
	240	}
	241
809a4336	242	//___________________________________________________________________
dbe3abbe	243	void AliHFEpidTRD::GetParameters(Double_t electronEff, Double_t *parameters){
	244	//
	245	// return parameter set for the given efficiency bin
	246	//
0792aa82	247	Int_t effbin = static_cast<Int_t>((electronEff - 0.7)/0.05);
dbe3abbe	248	memcpy(parameters, fThreshParams + effbin * 4, sizeof(Double_t) * 4);
809a4336	249	}
75d81601	250
75d81601	251	//___________________________________________________________________
3a72645a	252	Double_t AliHFEpidTRD::GetElectronLikelihood(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType){
	253	//
	254	// Get TRD likelihoods for ESD respectively AOD tracks
	255	//
	256	Double_t pidProbs[AliPID::kSPECIES]; memset(pidProbs, 0, sizeof(Double_t) * AliPID::kSPECIES);
	257	if(anaType == AliHFEpidObject::kESDanalysis){
	258	const AliESDtrack esdtrack = dynamic_cast<const AliESDtrack >(track);
	259	esdtrack->GetTRDpid(pidProbs);
	260	} else {
	261	const AliAODTrack aodtrack = dynamic_cast<const AliAODTrack >(track);
	262	fAODpid->MakeTRDPID(const_cast<AliAODTrack *>(aodtrack), pidProbs);
	263	}
	264	return pidProbs[AliPID::kElectron];
	265	}
	266
	267	//___________________________________________________________________
	268	Double_t AliHFEpidTRD::GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType){
	269	//
	270	// Get the Momentum in the TRD
	271	//
	272	Double_t p = 0.;
	273	if(anaType == AliHFEpidObject::kESDanalysis){
	274	const AliESDtrack esdtrack = dynamic_cast<const AliESDtrack >(track);
	275	p = esdtrack->GetOuterParam() ? esdtrack->GetOuterParam()->P() : esdtrack->P();
	276	} else {
	277	const AliAODTrack aodtrack = dynamic_cast<const AliAODTrack >(track);
	278	p = aodtrack->P();
	279	}
	280	return p;
	281	}
	282
	283	//___________________________________________________________________
	284	Double_t AliHFEpidTRD::GetTRDSignalV1(const AliESDtrack *track){
75d81601	285	//
	286	// Calculation of the TRD Signal via truncated mean
	287	// Method 1: Take all Slices available
	288	// cut out 0s
	289	// Order them in increasing order
	290	// Cut out the upper third
	291	// Calculate mean over the last 2/3 slices
	292	//
	293	const Int_t kNSlices = 48;
e3fc062d	294	AliDebug(3, Form("Number of Tracklets: %d\n", track->GetTRDpidQuality()));
75d81601	295	Double_t trdSlices[kNSlices], tmp[kNSlices];
	296	Int_t indices[48];
	297	Int_t icnt = 0;
	298	for(Int_t idet = 0; idet < 6; idet++)
	299	for(Int_t islice = 0; islice < 8; islice++){
	300	AliDebug(2, Form("Chamber[%d], Slice[%d]: TRDSlice = %f", idet, islice, track->GetTRDslice(idet, islice)));
	301	if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
	302	trdSlices[icnt++] = track->GetTRDslice(idet, islice);
	303	}
	304	AliDebug(1, Form("Number of Slices: %d\n", icnt));
	305	if(icnt < 6) return 0.; // We need at least 6 Slices for the truncated mean
	306	TMath::Sort(icnt, trdSlices, indices, kFALSE);
	307	memcpy(tmp, trdSlices, sizeof(Double_t) * icnt);
	308	for(Int_t ien = 0; ien < icnt; ien++)
	309	trdSlices[ien] = tmp[indices[ien]];
	310	Double_t trdSignal = TMath::Mean(static_cast<Int_t>(static_cast<Double_t>(icnt) * 2./3.), trdSlices);
	311	Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
e3fc062d	312	AliDebug(3, Form("PID Meth. 1: p[%f], TRDSignal[%f]", mom, trdSignal));
75d81601	313	return trdSignal;
	314	}
	315
	316	//___________________________________________________________________
3a72645a	317	Double_t AliHFEpidTRD::GetTRDSignalV2(const AliESDtrack *track){
75d81601	318	//
	319	// Calculation of the TRD Signal via truncated mean
	320	// Method 2: Take only first 5 slices per chamber
	321	// Order them in increasing order
	322	// Cut out upper half
	323	// Now do mean with the reamining 3 slices per chamber
	324	//
	325	Double_t trdSlicesLowTime[30], trdSlicesRemaining[32];
	326	Int_t indices[30];
	327	Int_t cntLowTime=0, cntRemaining = 0;
	328	for(Int_t idet = 0; idet < 6; idet++)
	329	for(Int_t islice = 0; islice < 8; islice++){
	330	if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
	331	if(islice < 5){
e3fc062d	332	AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
75d81601	333	trdSlicesLowTime[cntLowTime++] = track->GetTRDslice(idet, islice);
75d81601	334	} else{
e3fc062d	335	AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
75d81601	336	trdSlicesRemaining[cntRemaining++] = track->GetTRDslice(idet, islice);
	337	}
	338	}
	339	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
	340	TMath::Sort(cntLowTime, trdSlicesLowTime, indices, kFALSE);
	341	// Fill the second array with the lower half of the first time bins
	342	for(Int_t ien = 0; ien < static_cast<Int_t>(static_cast<Double_t>(cntLowTime) * 0.5); ien++)
	343	trdSlicesRemaining[cntRemaining++] = trdSlicesLowTime[indices[ien]];
	344	Double_t trdSignal = TMath::Mean(cntRemaining, trdSlicesRemaining);
	345	Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
e3fc062d	346	AliDebug(3, Form("PID Meth. 2: p[%f], TRDSignal[%f]", mom, trdSignal));
75d81601	347	return trdSignal;
75d81601	348	}