[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 <TList.h>
#include <TString.h>

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

#include "AliHFEpidTRD.h"

ClassImp(AliHFEpidTRD)

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

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

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

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

//______________________________________________________
Bool_t AliHFEpidTRD::InitializePID(){
  //
  // InitializePID: Load TRD thresholds and create the electron efficiency axis
  // to navigate 
  //
  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(p < 2.0) return 0;

  Double_t pidProbs[AliPID::kSPECIES];
  esdTrack->GetTRDpid(pidProbs);
  if(IsQAon())FillHistogramsLikelihood(0, p, pidProbs[AliPID::kElectron]);
  Double_t threshold = GetTRDthresholds(0.91, p);
  AliDebug(1, Form("Threshold: %f\n", threshold));
  if(IsQAon()) (dynamic_cast<TH2F *>(fContainer->At(kHistTRDthresholds)))->Fill(p, threshold);
  if(pidProbs[AliPID::kElectron] > threshold){
    if(IsQAon()) FillHistogramsLikelihood(1, p, pidProbs[AliPID::kElectron]);
    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
  //

  // 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::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(2, 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(1, Form("PID Meth. 1: p[%f], TRDSignal[%f]", mom, trdSignal));
  if(IsQAon() && mom > 0.) FillHistogramsTRDSignalV1(trdSignal, mom, mcPID);
  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(2, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
        trdSlicesLowTime[cntLowTime++] = track->GetTRDslice(idet, islice);
      } else{
        AliDebug(2, 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(1, Form("PID Meth. 2: p[%f], TRDSignal[%f]", mom, trdSignal));
  if(IsQAon() && mom > 0.) FillHistogramsTRDSignalV2(trdSignal, mom, mcPID);
  return trdSignal;
}

//___________________________________________________________________
void AliHFEpidTRD::FillHistogramsTRDSignalV1(Double_t signal, Double_t mom, Int_t species){
  //
  // Fill histograms for TRD Signal from Method 1 vs. p for different particle species
  //
  (dynamic_cast<TH2F *>(fContainer->At(kHistTRDSigV1)))->Fill(mom, signal);
  if(species < 0 || species >= AliPID::kSPECIES) return;
  (dynamic_cast<TH2F *>(fContainer->At(kHistOverallSpecies + species)))->Fill(mom, signal);
}

//___________________________________________________________________
void AliHFEpidTRD::FillHistogramsTRDSignalV2(Double_t signal, Double_t mom, Int_t species){
  //
  // Fill histograms for TRD Signal from Method 2 vs. p for different particle species
  //
  (dynamic_cast<TH2F *>(fContainer->At(kHistTRDSigV2)))->Fill(mom, signal);  
  if(species < 0 || species >= AliPID::kSPECIES) return;
  (dynamic_cast<TH2F *>(fContainer->At(kHistOverallSpecies + AliPID::kSPECIES + species)))->Fill(mom, signal);
}

//___________________________________________________________________
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 = 41;
  const Double_t kPtMin = 0.1;
  const Double_t kPtMax = 10.;
  const Int_t kSigBinsMeth1 = 100; 
  const Int_t kSigBinsMeth2 = 100;
  const Double_t kMinSig = 0.;
  const Double_t kMaxSigMeth1 = 10000.;
  const Double_t kMaxSigMeth2 = 10000.;
  
  if(!fContainer) fContainer = new TList;
  fContainer->SetName("fTRDqaHistograms");

  Double_t momentumBins[kMomentumBins];
  for(Int_t ibin = 0; ibin < kMomentumBins; ibin++)
    momentumBins[ibin] = static_cast<Double_t>(TMath::Power(10,TMath::Log10(kPtMin) + (TMath::Log10(kPtMax)-TMath::Log10(kPtMin))/(kMomentumBins-1)*static_cast<Double_t>(ibin)));
  // Likelihood Histograms
  fContainer->AddAt(new TH2F("fTRDlikeBefore", "TRD Electron Likelihood before cut", kMomentumBins - 1, momentumBins, 1000, 0, 1), kHistTRDlikeBefore);
  fContainer->AddAt(new TH2F("fTRDlikeAfter", "TRD Electron Likelihood after cut", kMomentumBins - 1, momentumBins, 1000, 0, 1), kHistTRDlikeAfter);
  fContainer->AddAt(new TH2F("fTRDthesholds", "TRD Electron thresholds", kMomentumBins - 1, momentumBins, 1000, 0, 1), kHistTRDthresholds);
  // Signal Histograms
  fContainer->AddAt(new TH2F("fTRDSigV1all", "TRD Signal (all particles, Method 1)", kMomentumBins - 1, momentumBins, kSigBinsMeth1, kMinSig, kMaxSigMeth1), kHistTRDSigV1);
  fContainer->AddAt(new TH2F("fTRDSigV2all", "TRD Signal (all particles, Method 2)", kMomentumBins - 1, momentumBins, kSigBinsMeth2, kMinSig, kMaxSigMeth2), kHistTRDSigV2);
  for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
    fContainer->AddAt(new TH2F(Form("fTRDSigV1%s", AliPID::ParticleName(ispec)), Form("TRD Signal (%s, Method 1)", AliPID::ParticleName(ispec)), kMomentumBins - 1, momentumBins, kSigBinsMeth1, kMinSig, kMaxSigMeth1), kHistOverallSpecies + ispec);
    fContainer->AddAt(new TH2F(Form("fTRDSigV2%s", AliPID::ParticleName(ispec)), Form("TRD Signal (%s, Method 2)", AliPID::ParticleName(ispec)), kMomentumBins - 1, momentumBins, kSigBinsMeth1, kMinSig, kMaxSigMeth2), kHistOverallSpecies + AliPID::kSPECIES + ispec);
  }
  l->AddLast(fContainer);
}

//___________________________________________________________________
void AliHFEpidTRD::FillHistogramsLikelihood(Int_t whenFilled, Float_t p, Float_t elProb){
  //
  // Fill Likelihood Histogram before respectively after decision
  //
  TH2F *histo = NULL;
  if(whenFilled)
    histo = dynamic_cast<TH2F *>(fContainer->At(kHistTRDlikeAfter));
  else
    histo = dynamic_cast<TH2F *>(fContainer->At(kHistTRDlikeBefore));
  if(!histo){
    AliError("QA histograms not found");
    return;
  }
  histo->Fill(p, elProb);
}
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>
50685501	26	#include <TList.h>
809a4336	27	#include <TString.h>
809a4336	28
722347d8	29	#include "AliAODTrack.h"
722347d8	30	#include "AliAODMCParticle.h"
809a4336	31	#include "AliESDtrack.h"
75d81601	32	#include "AliLog.h"
722347d8	33	#include "AliMCParticle.h"
809a4336	34	#include "AliPID.h"
	35
	36	#include "AliHFEpidTRD.h"
	37
	38	ClassImp(AliHFEpidTRD)
	39
75d81601	40	const Double_t AliHFEpidTRD::fgkVerySmall = 1e-12;
75d81601	41
809a4336	42	//___________________________________________________________________
	43	AliHFEpidTRD::AliHFEpidTRD(const char* name) :
	44	AliHFEpidBase(name)
809a4336	45	, fPIDMethod(kNN)
75d81601	46	, fContainer(0x0)
809a4336	47	{
	48	//
	49	// default constructor
	50	//
dbe3abbe	51	memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
809a4336	52	}
	53
	54	//___________________________________________________________________
	55	AliHFEpidTRD::AliHFEpidTRD(const AliHFEpidTRD &ref):
	56	AliHFEpidBase("")
809a4336	57	, fPIDMethod(kLQ)
75d81601	58	, fContainer(0x0)
809a4336	59	{
	60	//
	61	// Copy constructor
	62	//
dbe3abbe	63	memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
809a4336	64	ref.Copy(*this);
	65	}
	66
	67	//___________________________________________________________________
	68	AliHFEpidTRD &AliHFEpidTRD::operator=(const AliHFEpidTRD &ref){
	69	//
	70	// Assignment operator
	71	//
	72	if(this != &ref){
	73	ref.Copy(*this);
	74	}
	75	return *this;
	76	}
	77
	78	//___________________________________________________________________
	79	void AliHFEpidTRD::Copy(TObject &ref) const {
	80	//
	81	// Performs the copying of the object
	82	//
	83	AliHFEpidTRD &target = dynamic_cast<AliHFEpidTRD &>(ref);
	84
809a4336	85	target.fPIDMethod = fPIDMethod;
dbe3abbe	86	memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams);
75d81601	87	if(fContainer) target.fContainer = dynamic_cast<TList *>(fContainer->Clone());
809a4336	88	AliHFEpidBase::Copy(ref);
	89	}
	90
	91	//___________________________________________________________________
	92	AliHFEpidTRD::~AliHFEpidTRD(){
	93	//
	94	// Destructor
	95	//
75d81601	96	if(fContainer){
	97	fContainer->Clear();
	98	delete fContainer;
	99	}
809a4336	100	}
	101
	102	//______________________________________________________
	103	Bool_t AliHFEpidTRD::InitializePID(){
	104	//
	105	// InitializePID: Load TRD thresholds and create the electron efficiency axis
	106	// to navigate
	107	//
dbe3abbe	108	InitParameters();
809a4336	109	return kTRUE;
	110	}
	111
	112	//______________________________________________________
722347d8	113	Int_t AliHFEpidTRD::IsSelected(AliHFEpidObject *track){
809a4336	114	//
	115	// Does PID for TRD alone:
	116	// PID thresholds based on 90% Electron Efficiency level approximated by a linear
	117	// step function
	118	//
722347d8	119	if(track->fAnalysisType == AliHFEpidObject::kESDanalysis){
	120	AliESDtrack esdTrack = dynamic_cast<AliESDtrack >(track->fRecTrack);
	121	if(!esdTrack) return 0;
	122	AliMCParticle esdmc = dynamic_cast<AliMCParticle >(track->fMCtrack);
	123	return MakePIDesd(esdTrack, esdmc);
	124	} else {
	125	AliAODTrack aodTrack = dynamic_cast<AliAODTrack >(track->fRecTrack);
	126	if(!aodTrack) return 0;
	127	AliAODMCParticle aodmc = dynamic_cast<AliAODMCParticle >(track->fMCtrack);
	128	return MakePIDaod(aodTrack, aodmc);
	129	}
	130	}
	131
	132	//______________________________________________________
	133	Int_t AliHFEpidTRD::MakePIDaod(AliAODTrack * /aodTrack/, AliAODMCParticle * /aodmc/){
50685501	134	//
	135	// Does PID decision for AOD tracks as discribed above
	136	//
722347d8	137	AliError("AOD PID not yet implemented");
	138	return 0;
	139	}
	140
	141	//______________________________________________________
	142	Int_t AliHFEpidTRD::MakePIDesd(AliESDtrack esdTrack, AliMCParticle /mcTrack/){
50685501	143	//
	144	// Does PID decision for ESD tracks as discribed above
	145	//
809a4336	146	Double_t p = esdTrack->GetOuterParam() ? esdTrack->GetOuterParam()->P() : esdTrack->P();
722347d8	147	if(p < 2.0) return 0;
809a4336	148
809a4336	149	Double_t pidProbs[AliPID::kSPECIES];
809a4336	150	esdTrack->GetTRDpid(pidProbs);
0792aa82	151	if(IsQAon())FillHistogramsLikelihood(0, p, pidProbs[AliPID::kElectron]);
	152	Double_t threshold = GetTRDthresholds(0.91, p);
	153	AliDebug(1, Form("Threshold: %f\n", threshold));
	154	if(IsQAon()) (dynamic_cast<TH2F *>(fContainer->At(kHistTRDthresholds)))->Fill(p, threshold);
	155	if(pidProbs[AliPID::kElectron] > threshold){
	156	if(IsQAon()) FillHistogramsLikelihood(1, p, pidProbs[AliPID::kElectron]);
	157	return 11;
	158	}
	159	return 211;
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
	179	// Parameters for 6 Layers
	180	fThreshParams[0] = -0.001839; // 0.7 electron eff
	181	fThreshParams[1] = 0.000276;
	182	fThreshParams[2] = 0.044902;
	183	fThreshParams[3] = 1.726751;
	184	fThreshParams[4] = -0.002405; // 0.75 electron eff
	185	fThreshParams[5] = 0.000372;
	186	fThreshParams[6] = 0.061775;
	187	fThreshParams[7] = 1.739371;
	188	fThreshParams[8] = -0.003178; // 0.8 electron eff
	189	fThreshParams[9] = 0.000521;
	190	fThreshParams[10] = 0.087585;
	191	fThreshParams[11] = 1.749154;
	192	fThreshParams[12] = -0.004058; // 0.85 electron eff
	193	fThreshParams[13] = 0.000748;
	194	fThreshParams[14] = 0.129583;
	195	fThreshParams[15] = 1.782323;
	196	fThreshParams[16] = -0.004967; // 0.9 electron eff
	197	fThreshParams[17] = 0.001216;
	198	fThreshParams[18] = 0.210128;
	199	fThreshParams[19] = 1.807665;
	200	fThreshParams[20] = -0.000996; // 0.95 electron eff
	201	fThreshParams[21] = 0.002627;
	202	fThreshParams[22] = 0.409099;
	203	fThreshParams[23] = 1.787076;
809a4336	204	}
	205
	206	//___________________________________________________________________
dbe3abbe	207	void AliHFEpidTRD::GetParameters(Double_t electronEff, Double_t *parameters){
	208	//
	209	// return parameter set for the given efficiency bin
	210	//
0792aa82	211	Int_t effbin = static_cast<Int_t>((electronEff - 0.7)/0.05);
dbe3abbe	212	memcpy(parameters, fThreshParams + effbin * 4, sizeof(Double_t) * 4);
809a4336	213	}
75d81601	214
75d81601	215	//___________________________________________________________________
722347d8	216	Double_t AliHFEpidTRD::GetTRDSignalV1(AliESDtrack *track, Int_t mcPID){
75d81601	217	//
	218	// Calculation of the TRD Signal via truncated mean
	219	// Method 1: Take all Slices available
	220	// cut out 0s
	221	// Order them in increasing order
	222	// Cut out the upper third
	223	// Calculate mean over the last 2/3 slices
	224	//
	225	const Int_t kNSlices = 48;
	226	AliDebug(2, Form("Number of Tracklets: %d\n", track->GetTRDpidQuality()));
	227	Double_t trdSlices[kNSlices], tmp[kNSlices];
	228	Int_t indices[48];
	229	Int_t icnt = 0;
	230	for(Int_t idet = 0; idet < 6; idet++)
	231	for(Int_t islice = 0; islice < 8; islice++){
	232	AliDebug(2, Form("Chamber[%d], Slice[%d]: TRDSlice = %f", idet, islice, track->GetTRDslice(idet, islice)));
	233	if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
	234	trdSlices[icnt++] = track->GetTRDslice(idet, islice);
	235	}
	236	AliDebug(1, Form("Number of Slices: %d\n", icnt));
	237	if(icnt < 6) return 0.; // We need at least 6 Slices for the truncated mean
	238	TMath::Sort(icnt, trdSlices, indices, kFALSE);
	239	memcpy(tmp, trdSlices, sizeof(Double_t) * icnt);
	240	for(Int_t ien = 0; ien < icnt; ien++)
	241	trdSlices[ien] = tmp[indices[ien]];
	242	Double_t trdSignal = TMath::Mean(static_cast<Int_t>(static_cast<Double_t>(icnt) * 2./3.), trdSlices);
	243	Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
	244	AliDebug(1, Form("PID Meth. 1: p[%f], TRDSignal[%f]", mom, trdSignal));
722347d8	245	if(IsQAon() && mom > 0.) FillHistogramsTRDSignalV1(trdSignal, mom, mcPID);
75d81601	246	return trdSignal;
	247	}
	248
	249	//___________________________________________________________________
722347d8	250	Double_t AliHFEpidTRD::GetTRDSignalV2(AliESDtrack *track, Int_t mcPID){
75d81601	251	//
	252	// Calculation of the TRD Signal via truncated mean
	253	// Method 2: Take only first 5 slices per chamber
	254	// Order them in increasing order
	255	// Cut out upper half
	256	// Now do mean with the reamining 3 slices per chamber
	257	//
	258	Double_t trdSlicesLowTime[30], trdSlicesRemaining[32];
	259	Int_t indices[30];
	260	Int_t cntLowTime=0, cntRemaining = 0;
	261	for(Int_t idet = 0; idet < 6; idet++)
	262	for(Int_t islice = 0; islice < 8; islice++){
	263	if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
	264	if(islice < 5){
	265	AliDebug(2, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
	266	trdSlicesLowTime[cntLowTime++] = track->GetTRDslice(idet, islice);
	267	} else{
	268	AliDebug(2, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
	269	trdSlicesRemaining[cntRemaining++] = track->GetTRDslice(idet, islice);
	270	}
	271	}
	272	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
	273	TMath::Sort(cntLowTime, trdSlicesLowTime, indices, kFALSE);
	274	// Fill the second array with the lower half of the first time bins
	275	for(Int_t ien = 0; ien < static_cast<Int_t>(static_cast<Double_t>(cntLowTime) * 0.5); ien++)
	276	trdSlicesRemaining[cntRemaining++] = trdSlicesLowTime[indices[ien]];
	277	Double_t trdSignal = TMath::Mean(cntRemaining, trdSlicesRemaining);
	278	Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
	279	AliDebug(1, Form("PID Meth. 2: p[%f], TRDSignal[%f]", mom, trdSignal));
722347d8	280	if(IsQAon() && mom > 0.) FillHistogramsTRDSignalV2(trdSignal, mom, mcPID);
75d81601	281	return trdSignal;
	282	}
	283
	284	//___________________________________________________________________
	285	void AliHFEpidTRD::FillHistogramsTRDSignalV1(Double_t signal, Double_t mom, Int_t species){
	286	//
	287	// Fill histograms for TRD Signal from Method 1 vs. p for different particle species
	288	//
	289	(dynamic_cast<TH2F *>(fContainer->At(kHistTRDSigV1)))->Fill(mom, signal);
	290	if(species < 0 \|\| species >= AliPID::kSPECIES) return;
	291	(dynamic_cast<TH2F *>(fContainer->At(kHistOverallSpecies + species)))->Fill(mom, signal);
	292	}
	293
	294	//___________________________________________________________________
	295	void AliHFEpidTRD::FillHistogramsTRDSignalV2(Double_t signal, Double_t mom, Int_t species){
	296	//
	297	// Fill histograms for TRD Signal from Method 2 vs. p for different particle species
	298	//
	299	(dynamic_cast<TH2F *>(fContainer->At(kHistTRDSigV2)))->Fill(mom, signal);
	300	if(species < 0 \|\| species >= AliPID::kSPECIES) return;
	301	(dynamic_cast<TH2F *>(fContainer->At(kHistOverallSpecies + AliPID::kSPECIES + species)))->Fill(mom, signal);
	302	}
	303
75d81601	304	//___________________________________________________________________
	305	void AliHFEpidTRD::AddQAhistograms(TList *l){
	306	//
	307	// Adding QA histograms for the TRD PID
	308	// QA histograms are:
	309	// + TRD Signal from Meth. 1 vs p for all species
	310	// + TRD Signal from Meth. 2 vs p for all species
	311	// + For each species
	312	// - TRD Signal from Meth. 1 vs p
	313	// - TRD Signal from Meth. 2 vs p
	314	//
	315	const Int_t kMomentumBins = 41;
	316	const Double_t kPtMin = 0.1;
	317	const Double_t kPtMax = 10.;
	318	const Int_t kSigBinsMeth1 = 100;
	319	const Int_t kSigBinsMeth2 = 100;
	320	const Double_t kMinSig = 0.;
	321	const Double_t kMaxSigMeth1 = 10000.;
	322	const Double_t kMaxSigMeth2 = 10000.;
	323
	324	if(!fContainer) fContainer = new TList;
	325	fContainer->SetName("fTRDqaHistograms");
	326
	327	Double_t momentumBins[kMomentumBins];
	328	for(Int_t ibin = 0; ibin < kMomentumBins; ibin++)
	329	momentumBins[ibin] = static_cast<Double_t>(TMath::Power(10,TMath::Log10(kPtMin) + (TMath::Log10(kPtMax)-TMath::Log10(kPtMin))/(kMomentumBins-1)*static_cast<Double_t>(ibin)));
0792aa82	330	// Likelihood Histograms
	331	fContainer->AddAt(new TH2F("fTRDlikeBefore", "TRD Electron Likelihood before cut", kMomentumBins - 1, momentumBins, 1000, 0, 1), kHistTRDlikeBefore);
	332	fContainer->AddAt(new TH2F("fTRDlikeAfter", "TRD Electron Likelihood after cut", kMomentumBins - 1, momentumBins, 1000, 0, 1), kHistTRDlikeAfter);
	333	fContainer->AddAt(new TH2F("fTRDthesholds", "TRD Electron thresholds", kMomentumBins - 1, momentumBins, 1000, 0, 1), kHistTRDthresholds);
	334	// Signal Histograms
75d81601	335	fContainer->AddAt(new TH2F("fTRDSigV1all", "TRD Signal (all particles, Method 1)", kMomentumBins - 1, momentumBins, kSigBinsMeth1, kMinSig, kMaxSigMeth1), kHistTRDSigV1);
	336	fContainer->AddAt(new TH2F("fTRDSigV2all", "TRD Signal (all particles, Method 2)", kMomentumBins - 1, momentumBins, kSigBinsMeth2, kMinSig, kMaxSigMeth2), kHistTRDSigV2);
	337	for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
	338	fContainer->AddAt(new TH2F(Form("fTRDSigV1%s", AliPID::ParticleName(ispec)), Form("TRD Signal (%s, Method 1)", AliPID::ParticleName(ispec)), kMomentumBins - 1, momentumBins, kSigBinsMeth1, kMinSig, kMaxSigMeth1), kHistOverallSpecies + ispec);
	339	fContainer->AddAt(new TH2F(Form("fTRDSigV2%s", AliPID::ParticleName(ispec)), Form("TRD Signal (%s, Method 2)", AliPID::ParticleName(ispec)), kMomentumBins - 1, momentumBins, kSigBinsMeth1, kMinSig, kMaxSigMeth2), kHistOverallSpecies + AliPID::kSPECIES + ispec);
	340	}
	341	l->AddLast(fContainer);
	342	}
	343
0792aa82	344	//___________________________________________________________________
	345	void AliHFEpidTRD::FillHistogramsLikelihood(Int_t whenFilled, Float_t p, Float_t elProb){
	346	//
	347	// Fill Likelihood Histogram before respectively after decision
	348	//
	349	TH2F *histo = NULL;
	350	if(whenFilled)
	351	histo = dynamic_cast<TH2F *>(fContainer->At(kHistTRDlikeAfter));
	352	else
	353	histo = dynamic_cast<TH2F *>(fContainer->At(kHistTRDlikeBefore));
	354	if(!histo){
	355	AliError("QA histograms not found");
	356	return;
	357	}
	358	histo->Fill(p, elProb);
	359	}