[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 "AliAODPid.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)

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

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

//___________________________________________________________________
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);
  
  Bool_t defaultParameters = UseDefaultParameters();
  target.SetUseDefaultParameters(defaultParameters);
  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(UseDefaultParameters()){
    if(fPIDMethod == kLQ)
      InitParameters1DLQ();
    else
      InitParameters();
  }
  return kTRUE;
}

//______________________________________________________
Int_t AliHFEpidTRD::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const {
  //
  // 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) 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 AliHFEpidTRD::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 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::RenormalizeElPi(const Double_t * const likein, Double_t * const 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 AliHFEpidTRD::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 AliHFEpidTRD::GetElectronLikelihood(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const {
  //
  // 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);
    if(esdtrack) esdtrack->GetTRDpid(pidProbs);
  } else {
    const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
    if(aodtrack)fAODpid->MakeTRDPID(const_cast<AliAODTrack *>(aodtrack), pidProbs);
  }
  if(!IsRenormalizeElPi()) return pidProbs[AliPID::kElectron];
  Double_t probsNew[AliPID::kSPECIES];
  RenormalizeElPi(pidProbs, probsNew);
  return probsNew[AliPID::kElectron];
}

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

//___________________________________________________________________
Double_t AliHFEpidTRD::GetChargeLayer(const AliVParticle *track, UInt_t layer, AliHFEpidObject::AnalysisType_t anaType) const {
  //
  // Get the Charge in a single TRD layer
  //
  if(layer >= 6) return 0.;
  Double_t charge = 0.;
  if(anaType == AliHFEpidObject::kESDanalysis){
    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);
  } else {
    const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
    AliAODPid *aoddetpid = aodtrack ? aodtrack->GetDetPid() : NULL;
    if(aoddetpid)
      for(Int_t islice = 0; islice < aoddetpid->GetTRDnSlices(); islice++) charge += aoddetpid->GetTRDsignal()[layer * aoddetpid->GetTRDnSlices() + islice];
  }
  return charge;
}

//___________________________________________________________________
void AliHFEpidTRD::GetTRDmomenta(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType, Double_t *mom) const {
  //
  // Fill Array with momentum information at the TRD tracklet
  //
  if(anaType == AliHFEpidObject::kESDanalysis){
    const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
    if(esdtrack)
      for(Int_t itl = 0; itl < 6; itl++) 
        mom[itl] = esdtrack->GetTRDmomentum(itl);
  } else {
    const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
    AliAODPid *aoddetpid = aodtrack ? aodtrack->GetDetPid() : NULL;
    if(aoddetpid){
      Float_t *trdmom = aoddetpid->GetTRDmomentum();
      for(Int_t itl = 0; itl < 6; itl++){
        mom[itl] = trdmom[itl]; 
      }
    }
  }
}

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