[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 "AliAODTrack.h"
#include "AliAODMCParticle.h"
#include "AliESDtrack.h"
#include "AliESDpid.h"
#include "AliLog.h"
#include "AliMCParticle.h"
#include "AliPID.h"

#include "AliHFEcollection.h"
#include "AliHFEpidTRD.h"

ClassImp(AliHFEpidTRD)

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

//___________________________________________________________________
AliHFEpidTRD::AliHFEpidTRD() :
    AliHFEpidBase()
  , fMinP(1.)
  , fPIDMethod(kNN)
  , fContainer(0x0)
{
  //
  // default  constructor
  // 
  memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
}

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

//___________________________________________________________________
AliHFEpidTRD::AliHFEpidTRD(const AliHFEpidTRD &ref):
    AliHFEpidBase("")
  , fMinP(1.)
  , fPIDMethod(kLQ)
  , fContainer(0x0)
{
  //
  // 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;
  memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams);
  if(fContainer) target.fContainer = dynamic_cast<AliHFEcollection *>(fContainer->Clone());
  AliHFEpidBase::Copy(ref);
}

//___________________________________________________________________
AliHFEpidTRD::~AliHFEpidTRD(){
  //
  // Destructor
  //
  if(fContainer)
    delete fContainer;
}

//______________________________________________________
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){
  //
  // Does PID for TRD alone:
  // PID thresholds based on 90% Electron Efficiency level approximated by a linear 
  // step function
  //
  if(track->fAnalysisType == AliHFEpidObject::kESDanalysis){
      AliESDtrack *esdTrack = dynamic_cast<AliESDtrack *>(track->fRecTrack);
      if(!esdTrack) return 0;
      AliMCParticle *esdmc = dynamic_cast<AliMCParticle *>(track->fMCtrack);
      return MakePIDesd(esdTrack, esdmc);
   } else {
     AliAODTrack *aodTrack = dynamic_cast<AliAODTrack *>(track->fRecTrack);
     if(!aodTrack) return 0;
     AliAODMCParticle *aodmc = dynamic_cast<AliAODMCParticle *>(track->fMCtrack);
     return MakePIDaod(aodTrack, aodmc);
   }
}

//______________________________________________________
Int_t AliHFEpidTRD::MakePIDaod(AliAODTrack * /*aodTrack*/, AliAODMCParticle * /*aodmc*/){
  //
  // Does PID decision for AOD tracks as discribed above
  //
  AliError("AOD PID not yet implemented");
  return 0;
}

//______________________________________________________
Int_t AliHFEpidTRD::MakePIDesd(AliESDtrack *esdTrack, AliMCParticle * /*mcTrack*/){
  //
  // Does PID decision for ESD tracks as discribed above
  //
  Double_t p = esdTrack->GetOuterParam() ? esdTrack->GetOuterParam()->P() : esdTrack->P();
  if(IsQAon())
    FillStandardQA(0, esdTrack);
  if(p < fMinP) return 0;

  Double_t pidProbs[AliPID::kSPECIES];
  esdTrack->GetTRDpid(pidProbs);
  Double_t threshold = GetTRDthresholds(0.71, p);
  AliDebug(1, Form("Threshold: %f\n", threshold));
  if(IsQAon()) fContainer->Fill("fTRDthresholds", p, threshold);
  if(pidProbs[AliPID::kElectron] > threshold){
    if(IsQAon()) 
      FillStandardQA(1, esdTrack);
    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::GetTRDSignalV1(AliESDtrack *track, Int_t mcPID){
  //
  // 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));
  if(IsQAon() && mom > 0.) FillHistogramsTRDSignal(trdSignal, mom, mcPID, 1);
  return trdSignal;
}

//___________________________________________________________________
Double_t AliHFEpidTRD::GetTRDSignalV2(AliESDtrack *track, Int_t mcPID){
  //
  // 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));
  if(IsQAon() && mom > 0.) FillHistogramsTRDSignal(trdSignal, mom, mcPID, 2);
  return trdSignal;
}

//___________________________________________________________________
void AliHFEpidTRD::FillHistogramsTRDSignal(Double_t signal, Double_t p, Int_t species, UInt_t version){
  //
  // Fill histograms for TRD Signal from Method 2 vs. p for different particle species
  // Non-standard QA content
  //
  if(version == 0 || version > 2) return;
  if(species >= AliPID::kSPECIES || species < -1) species = -1;
  Double_t content[3] = {species, p, signal};
  Char_t histname[256]; sprintf(histname, "fTRDsignalV%d", version);
  THnSparseF *hsig = dynamic_cast<THnSparseF *>(fContainer->Get(histname));
  if(hsig) hsig->Fill(content);
}

//___________________________________________________________________
void AliHFEpidTRD::AddQAhistograms(TList *l){
  //
  // Adding QA histograms for the TRD PID
  // QA histograms are:
  // + TRD Signal from Meth. 1 vs p for all species
  // + TRD Signal from Meth. 2 vs p for all species
  // + For each species
  //    - TRD Signal from Meth. 1 vs p
  //    - TRD Signal from Meth. 2 vs p
  //
  const Int_t kMomentumBins = 100;
  const Double_t kPtMin = 0.1;
  const Double_t kPtMax = 20.;
  
  if(!fContainer) fContainer = new AliHFEcollection("fQAhistosTRD", "TRD QA histos");
  fContainer->CreateTH2F("fTRDlikeBefore", "TRD Electron Likelihood before cut; p [GeV/c]; TRD Electron Likelihood", kMomentumBins, kPtMin, kPtMax, 1000, 0, 1);
  fContainer->CreateTH2F("fTRDlikeAfter", "TRD Electron Likelihood after cut; p [GeV/c]; TRD Electron Likelihood", kMomentumBins, kPtMin, kPtMax, 1000, 0, 1);
  fContainer->CreateTH2F("fTRDthresholds", "TRD Electron Thresholds; p [GeV/c]; Thresholds", kMomentumBins, kPtMin, kPtMax, 1000, 0., 1.);
  fContainer->CreateTH2F("fTPCsignalBefore", "TPC dEdx Spectrum before TRD PID; p [GeV/c]; TPC Signal [a.u.]", kMomentumBins, kPtMin, kPtMax, 100, 0., 100.);
  fContainer->CreateTH2F("fTPCsignalAfter", "TPC dEdx Spectrum before TRD PID; p [GeV/c]; TPC Signal [a.u.]", kMomentumBins, kPtMin, kPtMax, 100, 0., 100.);
  fContainer->CreateTH2F("fTPCsigmaBefore", "TPC dEdx before TRD PID; p [GeV/c]; Normalized TPC distance to the electron line [n#sigma]", kMomentumBins, kPtMin, kPtMax, 100, -10., 10.);
  fContainer->CreateTH2F("fTPCsigmaAfter", "TPC dEdx Spectrum before TRD PID; p [GeV/c]; Normalized TPC distance to the electron line [n#sigma]", kMomentumBins, kPtMin, kPtMax, 100, -10., 10.);

  // Monitor THnSparse for TRD Signal
  const Int_t kBinsTRDsignal = 3; 
  Int_t nBins[kBinsTRDsignal] = {AliPID::kSPECIES +1, kMomentumBins, 100};
  Double_t binMin[kBinsTRDsignal] = {-1, kPtMin, 0};
  Double_t binMax[kBinsTRDsignal] = {AliPID::kSPECIES, kPtMax, 1000};
  fContainer->CreateTHnSparse("fTRDsignalV1", "TRD Signal V1", kBinsTRDsignal, nBins, binMin, binMax);
  fContainer->CreateTHnSparse("fTRDsignalV2", "TRD Signal V2", kBinsTRDsignal, nBins, binMin, binMax);

  // Make logatithmic binning
  TString hnames[7] = {"fTRDlikeBefore", "fTRDlikeAfter", "fTRDthresholds", "fTRDsignalBefore", "fTRDsignalAfter", "fTPCsigmaBefore", "fTPCsigmaAfter"};
  for(Int_t ihist = 0; ihist < 7; ihist++)
    fContainer->BinLogAxis(hnames[ihist].Data(), 0);
  fContainer->BinLogAxis("fTRDsignalV1", 1);
  fContainer->BinLogAxis("fTRDsignalV2", 1);

  l->Add(fContainer->GetList());
}

//___________________________________________________________________
void AliHFEpidTRD::FillStandardQA(Int_t whenFilled, AliESDtrack *esdTrack){
  //
  // Fill Likelihood Histogram before respectively after decision
  //
  Double_t p =  esdTrack->P();
  Double_t like[AliPID::kSPECIES];
  esdTrack->GetTRDpid(like);
  TString step;
  if(whenFilled)
    step = "After";
  else
    step = "Before";
  TString histos[3] = {"fTRDlike", "fTPCsignal", "fTPCsigma"};
  for(Int_t ihist = 0; ihist < 3; ihist++) histos[ihist] += step;
  fContainer->Fill(histos[0].Data(), esdTrack->P(), like[AliPID::kElectron]);
  fContainer->Fill(histos[1].Data(), p, esdTrack->GetTPCsignal());
  if(fESDpid)
    fContainer->Fill(histos[2].Data(), p, fESDpid->NumberOfSigmasTPC(esdTrack, AliPID::kElectron));
}
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
722347d8	30	#include "AliAODTrack.h"
722347d8	31	#include "AliAODMCParticle.h"
809a4336	32	#include "AliESDtrack.h"
e3fc062d	33	#include "AliESDpid.h"
75d81601	34	#include "AliLog.h"
722347d8	35	#include "AliMCParticle.h"
809a4336	36	#include "AliPID.h"
809a4336	37
e3fc062d	38	#include "AliHFEcollection.h"
809a4336	39	#include "AliHFEpidTRD.h"
	40
	41	ClassImp(AliHFEpidTRD)
	42
75d81601	43	const Double_t AliHFEpidTRD::fgkVerySmall = 1e-12;
75d81601	44
faee3b18	45	//___________________________________________________________________
	46	AliHFEpidTRD::AliHFEpidTRD() :
	47	AliHFEpidBase()
e3fc062d	48	, fMinP(1.)
faee3b18	49	, fPIDMethod(kNN)
	50	, fContainer(0x0)
	51	{
	52	//
	53	// default constructor
	54	//
	55	memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
	56	}
	57
809a4336	58	//___________________________________________________________________
	59	AliHFEpidTRD::AliHFEpidTRD(const char* name) :
	60	AliHFEpidBase(name)
e3fc062d	61	, fMinP(1.)
809a4336	62	, fPIDMethod(kNN)
75d81601	63	, fContainer(0x0)
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("")
e3fc062d	74	, fMinP(1.)
809a4336	75	, fPIDMethod(kLQ)
75d81601	76	, fContainer(0x0)
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;
dbe3abbe	105	memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams);
e3fc062d	106	if(fContainer) target.fContainer = dynamic_cast<AliHFEcollection *>(fContainer->Clone());
809a4336	107	AliHFEpidBase::Copy(ref);
	108	}
	109
	110	//___________________________________________________________________
	111	AliHFEpidTRD::~AliHFEpidTRD(){
	112	//
	113	// Destructor
	114	//
e3fc062d	115	if(fContainer)
75d81601	116	delete fContainer;
809a4336	117	}
	118
	119	//______________________________________________________
	120	Bool_t AliHFEpidTRD::InitializePID(){
	121	//
	122	// InitializePID: Load TRD thresholds and create the electron efficiency axis
	123	// to navigate
	124	//
e3fc062d	125	if(fPIDMethod == kLQ)
	126	InitParameters1DLQ();
	127	else
	128	InitParameters();
809a4336	129	return kTRUE;
	130	}
	131
	132	//______________________________________________________
722347d8	133	Int_t AliHFEpidTRD::IsSelected(AliHFEpidObject *track){
809a4336	134	//
	135	// Does PID for TRD alone:
	136	// PID thresholds based on 90% Electron Efficiency level approximated by a linear
	137	// step function
	138	//
722347d8	139	if(track->fAnalysisType == AliHFEpidObject::kESDanalysis){
	140	AliESDtrack esdTrack = dynamic_cast<AliESDtrack >(track->fRecTrack);
	141	if(!esdTrack) return 0;
	142	AliMCParticle esdmc = dynamic_cast<AliMCParticle >(track->fMCtrack);
	143	return MakePIDesd(esdTrack, esdmc);
	144	} else {
	145	AliAODTrack aodTrack = dynamic_cast<AliAODTrack >(track->fRecTrack);
	146	if(!aodTrack) return 0;
	147	AliAODMCParticle aodmc = dynamic_cast<AliAODMCParticle >(track->fMCtrack);
	148	return MakePIDaod(aodTrack, aodmc);
	149	}
	150	}
	151
	152	//______________________________________________________
	153	Int_t AliHFEpidTRD::MakePIDaod(AliAODTrack * /aodTrack/, AliAODMCParticle * /aodmc/){
50685501	154	//
	155	// Does PID decision for AOD tracks as discribed above
	156	//
722347d8	157	AliError("AOD PID not yet implemented");
	158	return 0;
	159	}
	160
	161	//______________________________________________________
	162	Int_t AliHFEpidTRD::MakePIDesd(AliESDtrack esdTrack, AliMCParticle /mcTrack/){
50685501	163	//
	164	// Does PID decision for ESD tracks as discribed above
	165	//
809a4336	166	Double_t p = esdTrack->GetOuterParam() ? esdTrack->GetOuterParam()->P() : esdTrack->P();
e3fc062d	167	if(IsQAon())
	168	FillStandardQA(0, esdTrack);
	169	if(p < fMinP) return 0;
809a4336	170
809a4336	171	Double_t pidProbs[AliPID::kSPECIES];
809a4336	172	esdTrack->GetTRDpid(pidProbs);
e3fc062d	173	Double_t threshold = GetTRDthresholds(0.71, p);
0792aa82	174	AliDebug(1, Form("Threshold: %f\n", threshold));
e3fc062d	175	if(IsQAon()) fContainer->Fill("fTRDthresholds", p, threshold);
0792aa82	176	if(pidProbs[AliPID::kElectron] > threshold){
e3fc062d	177	if(IsQAon())
e3fc062d	178	FillStandardQA(1, esdTrack);
0792aa82	179	return 11;
	180	}
	181	return 211;
809a4336	182	}
	183
	184	//___________________________________________________________________
dbe3abbe	185	Double_t AliHFEpidTRD::GetTRDthresholds(Double_t electronEff, Double_t p){
	186	//
	187	// Return momentum dependent and electron efficiency dependent TRD thresholds
	188	//
	189	Double_t params[4];
	190	GetParameters(electronEff, params);
0792aa82	191	Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
0792aa82	192	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	193	}
	194
	195	//___________________________________________________________________
dbe3abbe	196	void AliHFEpidTRD::InitParameters(){
	197	//
	198	// Fill the Parameters into an array
	199	//
	200
e3fc062d	201	AliDebug(2, "Loading threshold Parameter");
dbe3abbe	202	// Parameters for 6 Layers
	203	fThreshParams[0] = -0.001839; // 0.7 electron eff
	204	fThreshParams[1] = 0.000276;
	205	fThreshParams[2] = 0.044902;
	206	fThreshParams[3] = 1.726751;
	207	fThreshParams[4] = -0.002405; // 0.75 electron eff
	208	fThreshParams[5] = 0.000372;
	209	fThreshParams[6] = 0.061775;
	210	fThreshParams[7] = 1.739371;
	211	fThreshParams[8] = -0.003178; // 0.8 electron eff
	212	fThreshParams[9] = 0.000521;
	213	fThreshParams[10] = 0.087585;
	214	fThreshParams[11] = 1.749154;
	215	fThreshParams[12] = -0.004058; // 0.85 electron eff
	216	fThreshParams[13] = 0.000748;
	217	fThreshParams[14] = 0.129583;
	218	fThreshParams[15] = 1.782323;
	219	fThreshParams[16] = -0.004967; // 0.9 electron eff
	220	fThreshParams[17] = 0.001216;
	221	fThreshParams[18] = 0.210128;
	222	fThreshParams[19] = 1.807665;
	223	fThreshParams[20] = -0.000996; // 0.95 electron eff
	224	fThreshParams[21] = 0.002627;
	225	fThreshParams[22] = 0.409099;
	226	fThreshParams[23] = 1.787076;
809a4336	227	}
809a4336	228
e3fc062d	229	//___________________________________________________________________
	230	void AliHFEpidTRD::InitParameters1DLQ(){
	231	//
	232	// Init Parameters for 1DLQ PID (M. Fasel, Sept. 6th, 2010)
	233	//
	234
	235	// Parameters for 6 Layers
	236	AliDebug(2, Form("Loading threshold parameter for Method 1DLQ"));
	237	fThreshParams[0] = -0.02241; // 0.7 electron eff
	238	fThreshParams[1] = 0.05043;
	239	fThreshParams[2] = 0.7925;
	240	fThreshParams[3] = 2.625;
	241	fThreshParams[4] = 0.07438; // 0.75 electron eff
	242	fThreshParams[5] = 0.05158;
	243	fThreshParams[6] = 2.864;
	244	fThreshParams[7] = 4.356;
	245	fThreshParams[8] = 0.1977; // 0.8 electron eff
	246	fThreshParams[9] = 0.05956;
	247	fThreshParams[10] = 2.853;
	248	fThreshParams[11] = 3.713;
	249	fThreshParams[12] = 0.5206; // 0.85 electron eff
	250	fThreshParams[13] = 0.03077;
	251	fThreshParams[14] = 2.966;
	252	fThreshParams[15] = 4.07;
	253	fThreshParams[16] = 0.8808; // 0.9 electron eff
	254	fThreshParams[17] = 0.002092;
	255	fThreshParams[18] = 1.17;
	256	fThreshParams[19] = 4.506;
	257	fThreshParams[20] = 1.; // 0.95 electron eff
	258	fThreshParams[21] = 0.;
	259	fThreshParams[22] = 0.;
	260	fThreshParams[23] = 0.;
	261
	262	}
	263
809a4336	264	//___________________________________________________________________
dbe3abbe	265	void AliHFEpidTRD::GetParameters(Double_t electronEff, Double_t *parameters){
	266	//
	267	// return parameter set for the given efficiency bin
	268	//
0792aa82	269	Int_t effbin = static_cast<Int_t>((electronEff - 0.7)/0.05);
dbe3abbe	270	memcpy(parameters, fThreshParams + effbin * 4, sizeof(Double_t) * 4);
809a4336	271	}
75d81601	272
75d81601	273	//___________________________________________________________________
722347d8	274	Double_t AliHFEpidTRD::GetTRDSignalV1(AliESDtrack *track, Int_t mcPID){
75d81601	275	//
	276	// Calculation of the TRD Signal via truncated mean
	277	// Method 1: Take all Slices available
	278	// cut out 0s
	279	// Order them in increasing order
	280	// Cut out the upper third
	281	// Calculate mean over the last 2/3 slices
	282	//
	283	const Int_t kNSlices = 48;
e3fc062d	284	AliDebug(3, Form("Number of Tracklets: %d\n", track->GetTRDpidQuality()));
75d81601	285	Double_t trdSlices[kNSlices], tmp[kNSlices];
	286	Int_t indices[48];
	287	Int_t icnt = 0;
	288	for(Int_t idet = 0; idet < 6; idet++)
	289	for(Int_t islice = 0; islice < 8; islice++){
	290	AliDebug(2, Form("Chamber[%d], Slice[%d]: TRDSlice = %f", idet, islice, track->GetTRDslice(idet, islice)));
	291	if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
	292	trdSlices[icnt++] = track->GetTRDslice(idet, islice);
	293	}
	294	AliDebug(1, Form("Number of Slices: %d\n", icnt));
	295	if(icnt < 6) return 0.; // We need at least 6 Slices for the truncated mean
	296	TMath::Sort(icnt, trdSlices, indices, kFALSE);
	297	memcpy(tmp, trdSlices, sizeof(Double_t) * icnt);
	298	for(Int_t ien = 0; ien < icnt; ien++)
	299	trdSlices[ien] = tmp[indices[ien]];
	300	Double_t trdSignal = TMath::Mean(static_cast<Int_t>(static_cast<Double_t>(icnt) * 2./3.), trdSlices);
	301	Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
e3fc062d	302	AliDebug(3, Form("PID Meth. 1: p[%f], TRDSignal[%f]", mom, trdSignal));
e3fc062d	303	if(IsQAon() && mom > 0.) FillHistogramsTRDSignal(trdSignal, mom, mcPID, 1);
75d81601	304	return trdSignal;
	305	}
	306
	307	//___________________________________________________________________
722347d8	308	Double_t AliHFEpidTRD::GetTRDSignalV2(AliESDtrack *track, Int_t mcPID){
75d81601	309	//
	310	// Calculation of the TRD Signal via truncated mean
	311	// Method 2: Take only first 5 slices per chamber
	312	// Order them in increasing order
	313	// Cut out upper half
	314	// Now do mean with the reamining 3 slices per chamber
	315	//
	316	Double_t trdSlicesLowTime[30], trdSlicesRemaining[32];
	317	Int_t indices[30];
	318	Int_t cntLowTime=0, cntRemaining = 0;
	319	for(Int_t idet = 0; idet < 6; idet++)
	320	for(Int_t islice = 0; islice < 8; islice++){
	321	if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
	322	if(islice < 5){
e3fc062d	323	AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
75d81601	324	trdSlicesLowTime[cntLowTime++] = track->GetTRDslice(idet, islice);
75d81601	325	} else{
e3fc062d	326	AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
75d81601	327	trdSlicesRemaining[cntRemaining++] = track->GetTRDslice(idet, islice);
	328	}
	329	}
	330	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
	331	TMath::Sort(cntLowTime, trdSlicesLowTime, indices, kFALSE);
	332	// Fill the second array with the lower half of the first time bins
	333	for(Int_t ien = 0; ien < static_cast<Int_t>(static_cast<Double_t>(cntLowTime) * 0.5); ien++)
	334	trdSlicesRemaining[cntRemaining++] = trdSlicesLowTime[indices[ien]];
	335	Double_t trdSignal = TMath::Mean(cntRemaining, trdSlicesRemaining);
	336	Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
e3fc062d	337	AliDebug(3, Form("PID Meth. 2: p[%f], TRDSignal[%f]", mom, trdSignal));
e3fc062d	338	if(IsQAon() && mom > 0.) FillHistogramsTRDSignal(trdSignal, mom, mcPID, 2);
75d81601	339	return trdSignal;
	340	}
	341
	342	//___________________________________________________________________
e3fc062d	343	void AliHFEpidTRD::FillHistogramsTRDSignal(Double_t signal, Double_t p, Int_t species, UInt_t version){
75d81601	344	//
75d81601	345	// Fill histograms for TRD Signal from Method 2 vs. p for different particle species
e3fc062d	346	// Non-standard QA content
	347	//
	348	if(version == 0 \|\| version > 2) return;
	349	if(species >= AliPID::kSPECIES \|\| species < -1) species = -1;
	350	Double_t content[3] = {species, p, signal};
	351	Char_t histname[256]; sprintf(histname, "fTRDsignalV%d", version);
	352	THnSparseF hsig = dynamic_cast<THnSparseF >(fContainer->Get(histname));
	353	if(hsig) hsig->Fill(content);
75d81601	354	}
75d81601	355
75d81601	356	//___________________________________________________________________
	357	void AliHFEpidTRD::AddQAhistograms(TList *l){
	358	//
	359	// Adding QA histograms for the TRD PID
	360	// QA histograms are:
	361	// + TRD Signal from Meth. 1 vs p for all species
	362	// + TRD Signal from Meth. 2 vs p for all species
	363	// + For each species
	364	// - TRD Signal from Meth. 1 vs p
	365	// - TRD Signal from Meth. 2 vs p
	366	//
e3fc062d	367	const Int_t kMomentumBins = 100;
75d81601	368	const Double_t kPtMin = 0.1;
e3fc062d	369	const Double_t kPtMax = 20.;
75d81601	370
e3fc062d	371	if(!fContainer) fContainer = new AliHFEcollection("fQAhistosTRD", "TRD QA histos");
	372	fContainer->CreateTH2F("fTRDlikeBefore", "TRD Electron Likelihood before cut; p [GeV/c]; TRD Electron Likelihood", kMomentumBins, kPtMin, kPtMax, 1000, 0, 1);
	373	fContainer->CreateTH2F("fTRDlikeAfter", "TRD Electron Likelihood after cut; p [GeV/c]; TRD Electron Likelihood", kMomentumBins, kPtMin, kPtMax, 1000, 0, 1);
	374	fContainer->CreateTH2F("fTRDthresholds", "TRD Electron Thresholds; p [GeV/c]; Thresholds", kMomentumBins, kPtMin, kPtMax, 1000, 0., 1.);
	375	fContainer->CreateTH2F("fTPCsignalBefore", "TPC dEdx Spectrum before TRD PID; p [GeV/c]; TPC Signal [a.u.]", kMomentumBins, kPtMin, kPtMax, 100, 0., 100.);
	376	fContainer->CreateTH2F("fTPCsignalAfter", "TPC dEdx Spectrum before TRD PID; p [GeV/c]; TPC Signal [a.u.]", kMomentumBins, kPtMin, kPtMax, 100, 0., 100.);
	377	fContainer->CreateTH2F("fTPCsigmaBefore", "TPC dEdx before TRD PID; p [GeV/c]; Normalized TPC distance to the electron line [n#sigma]", kMomentumBins, kPtMin, kPtMax, 100, -10., 10.);
	378	fContainer->CreateTH2F("fTPCsigmaAfter", "TPC dEdx Spectrum before TRD PID; p [GeV/c]; Normalized TPC distance to the electron line [n#sigma]", kMomentumBins, kPtMin, kPtMax, 100, -10., 10.);
75d81601	379
e3fc062d	380	// Monitor THnSparse for TRD Signal
	381	const Int_t kBinsTRDsignal = 3;
	382	Int_t nBins[kBinsTRDsignal] = {AliPID::kSPECIES +1, kMomentumBins, 100};
	383	Double_t binMin[kBinsTRDsignal] = {-1, kPtMin, 0};
	384	Double_t binMax[kBinsTRDsignal] = {AliPID::kSPECIES, kPtMax, 1000};
	385	fContainer->CreateTHnSparse("fTRDsignalV1", "TRD Signal V1", kBinsTRDsignal, nBins, binMin, binMax);
	386	fContainer->CreateTHnSparse("fTRDsignalV2", "TRD Signal V2", kBinsTRDsignal, nBins, binMin, binMax);
	387
	388	// Make logatithmic binning
	389	TString hnames[7] = {"fTRDlikeBefore", "fTRDlikeAfter", "fTRDthresholds", "fTRDsignalBefore", "fTRDsignalAfter", "fTPCsigmaBefore", "fTPCsigmaAfter"};
	390	for(Int_t ihist = 0; ihist < 7; ihist++)
	391	fContainer->BinLogAxis(hnames[ihist].Data(), 0);
	392	fContainer->BinLogAxis("fTRDsignalV1", 1);
	393	fContainer->BinLogAxis("fTRDsignalV2", 1);
	394
	395	l->Add(fContainer->GetList());
75d81601	396	}
75d81601	397
0792aa82	398	//___________________________________________________________________
e3fc062d	399	void AliHFEpidTRD::FillStandardQA(Int_t whenFilled, AliESDtrack *esdTrack){
0792aa82	400	//
	401	// Fill Likelihood Histogram before respectively after decision
	402	//
e3fc062d	403	Double_t p = esdTrack->P();
	404	Double_t like[AliPID::kSPECIES];
	405	esdTrack->GetTRDpid(like);
	406	TString step;
0792aa82	407	if(whenFilled)
e3fc062d	408	step = "After";
0792aa82	409	else
e3fc062d	410	step = "Before";
	411	TString histos[3] = {"fTRDlike", "fTPCsignal", "fTPCsigma"};
	412	for(Int_t ihist = 0; ihist < 3; ihist++) histos[ihist] += step;
	413	fContainer->Fill(histos[0].Data(), esdTrack->P(), like[AliPID::kElectron]);
	414	fContainer->Fill(histos[1].Data(), p, esdTrack->GetTPCsignal());
	415	if(fESDpid)
	416	fContainer->Fill(histos[2].Data(), p, fESDpid->NumberOfSigmasTPC(esdTrack, AliPID::kElectron));
0792aa82	417	}