[u/mrichter/AliRoot.git] / STEER / STEERBase / AliTRDPIDResponse.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.                  *
**************************************************************************/
//
// PID Response class for the TRD detector
// Based on 1D Likelihood approach
// Calculation of probabilities using Bayesian approach
// Attention: This method is only used to separate electrons from pions
//
// Authors:
//  Markus Fasel <M.Fasel@gsi.de>
//  Anton Andronic <A.Andronic@gsi.de>
//
#include <TAxis.h>
#include <TClass.h>
#include <TDirectory.h>
#include <TFile.h>
#include <TH1.h>
#include <TKey.h>
#include <TMath.h>
#include <TObjArray.h>
#include <TROOT.h> 
#include <TString.h>
#include <TSystem.h>

#include "AliLog.h"

#include "AliTRDPIDResponseObject.h"
#include "AliTRDPIDResponse.h"
#include "AliTRDTKDInterpolator.h"

ClassImp(AliTRDPIDResponse)

//____________________________________________________________
AliTRDPIDResponse::AliTRDPIDResponse():
  TObject()
  ,fkPIDResponseObject(NULL)
  ,fGainNormalisationFactor(1.)
{
  //
  // Default constructor
  //
}

//____________________________________________________________
AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
  TObject(ref)
  ,fkPIDResponseObject(NULL)
  ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
{
  //
  // Copy constructor
  //
}

//____________________________________________________________
AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
  //
  // Assignment operator
  //
  if(this == &ref) return *this;
  
  // Make copy
  TObject::operator=(ref);
  fGainNormalisationFactor = ref.fGainNormalisationFactor;
  fkPIDResponseObject = ref.fkPIDResponseObject;
  
  return *this;
}

//____________________________________________________________
AliTRDPIDResponse::~AliTRDPIDResponse(){
  //
  // Destructor
  //
    if(IsOwner()) delete fkPIDResponseObject;
}

//____________________________________________________________
Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
  //
  // Load References into the toolkit
  //
  AliDebug(1, "Loading reference histos from root file");
  TDirectory *owd = gDirectory;// store old working directory
  
  if(!filename)
    filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT"));
  TFile *in = TFile::Open(filename);
  if(!in){
    AliError("Ref file not available.");
    return kFALSE;
  }
  
  gROOT->cd();
  fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(in->Get("TRDPIDResponse")->Clone());
  in->Close(); delete in;
  owd->cd();
  SetBit(kIsOwner, kTRUE);
  AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDResponseObject->GetNumberOfMomentumBins()));
  return kTRUE;
}

//____________________________________________________________
Bool_t AliTRDPIDResponse::GetResponse(Int_t n, const Double_t * const dedx, const Float_t * const p, Double_t prob[AliPID::kSPECIES],ETRDPIDMethod PIDmethod,Bool_t kNorm) const
{
  //
  // Calculate TRD likelihood values for the track based on dedx and 
  // momentum values. The likelihoods are calculated by query the 
  // reference data depending on the PID method selected
  //
  // Input parameter :
  //   n - number of dedx slices/chamber
  //   dedx - array of dedx slices organized layer wise
  //   p - array of momentum measurements organized layer wise
  // 
  // Return parameters
  //   prob - probabilities allocated by TRD for particle specis
  //   kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
  // 
  // Return value
  //   true if calculation success
    //
    AliDebug(3,Form(" Response for PID method: %d",PIDmethod));


    if(!fkPIDResponseObject){
	AliWarning("Missing reference data. PID calculation not possible.");
	return kFALSE;
    }

    for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
    Double_t prLayer[AliPID::kSPECIES];
    Double_t dE[10], s(0.);
    for(Int_t il(kNlayer); il--;){
	memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
	if(!CookdEdx(n, &dedx[il*n], &dE[0],PIDmethod)) continue;

	s=0.;
        Bool_t filled=kTRUE;
	for(Int_t is(AliPID::kSPECIES); is--;){
	    if((PIDmethod==kLQ2D)&&(!(is==0||is==2)))continue;
	    if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], &dE[0],PIDmethod);
	    AliDebug(3, Form("Probability for Species %d in Layer %d: %e", is, il, prLayer[is]));
	    if(prLayer[is]<1.e-30){
		AliDebug(2, Form("Null for species %d species prob layer[%d].",is,il));
		filled=kFALSE;
		break;
	    }
	    s+=prLayer[is];
	}
	if(!filled){
	    continue;
	}
	if(s<1.e-30){
	    AliDebug(2, Form("Null all species prob layer[%d].", il));
	    continue;
	}
	for(Int_t is(AliPID::kSPECIES); is--;){
	    if(kNorm) prLayer[is] /= s;
	    prob[is] *= prLayer[is];
	}
    }
    if(!kNorm) return kTRUE;

    s=0.;
    for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
    if(s<1.e-30){
	AliDebug(2, "Null total prob.");
	return kFALSE;
    }
    for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
    return kTRUE;
}

//____________________________________________________________
Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t *dEdx,ETRDPIDMethod PIDmethod) const {
  //
  // Get the non-normalized probability for a certain particle species as coming
  // from the reference histogram
  // Interpolation between momentum bins
  //
  AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));

  Double_t probLayer = 0.;
  Float_t pLower, pUpper;
	
  switch(PIDmethod){
  case kNN: // NN
      break;
  case kLQ2D: // 2D LQ
      {
	  if(species==0||species==2){ // references only for electrons and pions
	      Double_t error;
	      Double_t point[kNslicesLQ2D];
	      for(Int_t idim=0;idim<kNslicesLQ2D;idim++){point[idim]=dEdx[idim];}

	      AliTRDTKDInterpolator *refLower = dynamic_cast<AliTRDTKDInterpolator*>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ2D));

	      if(refLower){
		  refLower->Eval(point,probLayer,error);
	      }
	      else {
		  AliError("No references available");
	      }
	      AliDebug(2,Form("Eval 2D Q0 %f Q1 %f P %e Err %e",point[0],point[1],probLayer,error));
	  }
      }
      break;
  case kLQ1D: // 1D LQ
      {
	  TH1 *refUpper = dynamic_cast<TH1 *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ1D)),
	      *refLower = dynamic_cast<TH1 *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ1D));
	  // Do Interpolation exept for underflow and overflow
	  if(refLower && refUpper){
	      Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
	      Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));

	      probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
	  } else if(refLower){
	      // underflow
	      probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
	  } else if(refUpper){
	      // overflow
	      probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
	  } else {
	      AliError("No references available");
	  }
	  AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
      }
      break;
  default:
      break;
  }
  return probLayer;
}

//____________________________________________________________
void AliTRDPIDResponse::SetOwner(){
  //
  // Make Deep Copy of the Reference Histograms
  //
    if(!fkPIDResponseObject || IsOwner()) return;
    const AliTRDPIDResponseObject *tmp = fkPIDResponseObject;
    fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(tmp->Clone());
    SetBit(kIsOwner, kTRUE);
}

//____________________________________________________________
Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out,ETRDPIDMethod PIDmethod) const
{
    //
    // Recalculate dE/dx
    //
  switch(PIDmethod){
  case kNN: // NN 
      break;
  case kLQ2D: // 2D LQ
      out[0]=0;
      out[1]=0;
      for(Int_t islice = 0; islice < nSlice; islice++){
	  if(islice<fkPIDResponseObject->GetNSlicesQ0())out[0]+= in[islice];
	  else out[1]+= in[islice];
      }
      if(out[0] < 1e-6) return kFALSE;
      AliDebug(3,Form("CookdEdx Q0 %f Q1 %f",out[0],out[1]));
      break;
  case kLQ1D: // 1D LQ
      out[0]= 0.;
      for(Int_t islice = 0; islice < nSlice; islice++)
	  if(in[islice] > 0) out[0] += in[islice] * fGainNormalisationFactor;   // Protect against negative values for slices having no dE/dx information
      if(out[0] < 1e-6) return kFALSE;
      AliDebug(3,Form("CookdEdx dEdx %f",out[0]));
      break;

  default:
      return kFALSE;
  }
  return kTRUE;
}

//____________________________________________________________
Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level,Double_t centrality,ETRDPIDMethod PIDmethod) const {
    //
    // Check whether particle is identified as electron assuming a certain electron efficiency level
    // Only electron and pion hypothesis is taken into account
    //
    // Inputs:
    //         Number of tracklets
    //         Likelihood values
    //         Momentum
    //         Electron efficiency level
    //
    // If the function fails when the params are not accessible, the function returns true
    //
  if(!fkPIDResponseObject){
    AliError("No PID Param object available");
    return kTRUE;
  } 
  Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]);
  Double_t params[4];
  if(!fkPIDResponseObject->GetThresholdParameters(nTracklets, level, params,centrality,PIDmethod)){
    AliError("No Params found for the given configuration");
    return kTRUE;
  }
  Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
  if(probEle > TMath::Max(TMath::Min(threshold, 0.99), 0.2)) return kTRUE; // truncate the threshold upperwards to 0.999 and lowerwards to 0.2 and exclude unphysical values
  return kFALSE;
}

//____________________________________________________________
Bool_t AliTRDPIDResponse::SetPIDResponseObject(const AliTRDPIDResponseObject * obj){

    fkPIDResponseObject = obj;
    if((AliLog::GetDebugLevel("",IsA()->GetName()))>0)fkPIDResponseObject->Print("");
    return kTRUE;
}
Commit	Line	Data
ffb1ee30	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15	//
	16	// PID Response class for the TRD detector
	17	// Based on 1D Likelihood approach
	18	// Calculation of probabilities using Bayesian approach
	19	// Attention: This method is only used to separate electrons from pions
	20	//
	21	// Authors:
	22	// Markus Fasel <M.Fasel@gsi.de>
	23	// Anton Andronic <A.Andronic@gsi.de>
	24	//
ffb1ee30	25	#include <TAxis.h>
ea235c90	26	#include <TClass.h>
ea235c90	27	#include <TDirectory.h>
ffb1ee30	28	#include <TFile.h>
ffb1ee30	29	#include <TH1.h>
ce5d6908	30	#include <TKey.h>
ea235c90	31	#include <TMath.h>
ffb1ee30	32	#include <TObjArray.h>
ce5d6908	33	#include <TROOT.h>
ea235c90	34	#include <TString.h>
ffb1ee30	35	#include <TSystem.h>
	36
	37	#include "AliLog.h"
	38
db0e2c5f	39	#include "AliTRDPIDResponseObject.h"
ffb1ee30	40	#include "AliTRDPIDResponse.h"
db0e2c5f	41	#include "AliTRDTKDInterpolator.h"
ffb1ee30	42
	43	ClassImp(AliTRDPIDResponse)
	44
ffb1ee30	45	//____________________________________________________________
e0de37e9	46	AliTRDPIDResponse::AliTRDPIDResponse():
e0de37e9	47	TObject()
db0e2c5f	48	,fkPIDResponseObject(NULL)
e0de37e9	49	,fGainNormalisationFactor(1.)
ffb1ee30	50	{
e0de37e9	51	//
	52	// Default constructor
	53	//
ffb1ee30	54	}
	55
	56	//____________________________________________________________
	57	AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
e0de37e9	58	TObject(ref)
db0e2c5f	59	,fkPIDResponseObject(NULL)
e0de37e9	60	,fGainNormalisationFactor(ref.fGainNormalisationFactor)
ffb1ee30	61	{
e0de37e9	62	//
e0de37e9	63	// Copy constructor
e0de37e9	64	//
ffb1ee30	65	}
	66
	67	//____________________________________________________________
	68	AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
e0de37e9	69	//
e0de37e9	70	// Assignment operator
e0de37e9	71	//
	72	if(this == &ref) return *this;
	73
	74	// Make copy
	75	TObject::operator=(ref);
51a0ce25	76	fGainNormalisationFactor = ref.fGainNormalisationFactor;
db0e2c5f	77	fkPIDResponseObject = ref.fkPIDResponseObject;
e0de37e9	78
e0de37e9	79	return *this;
ffb1ee30	80	}
	81
	82	//____________________________________________________________
	83	AliTRDPIDResponse::~AliTRDPIDResponse(){
e0de37e9	84	//
	85	// Destructor
	86	//
db0e2c5f	87	if(IsOwner()) delete fkPIDResponseObject;
ffb1ee30	88	}
	89
	90	//____________________________________________________________
db0e2c5f	91	Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
e0de37e9	92	//
	93	// Load References into the toolkit
	94	//
	95	AliDebug(1, "Loading reference histos from root file");
2ad33025	96	TDirectory *owd = gDirectory;// store old working directory
e0de37e9	97
	98	if(!filename)
	99	filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT"));
	100	TFile *in = TFile::Open(filename);
	101	if(!in){
	102	AliError("Ref file not available.");
	103	return kFALSE;
	104	}
	105
	106	gROOT->cd();
db0e2c5f	107	fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(in->Get("TRDPIDResponse")->Clone());
51a0ce25	108	in->Close(); delete in;
e0de37e9	109	owd->cd();
e0de37e9	110	SetBit(kIsOwner, kTRUE);
db0e2c5f	111	AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDResponseObject->GetNumberOfMomentumBins()));
e0de37e9	112	return kTRUE;
ffb1ee30	113	}
	114
	115	//____________________________________________________________
bd58d4b9	116	Bool_t AliTRDPIDResponse::GetResponse(Int_t n, const Double_t * const dedx, const Float_t * const p, Double_t prob[AliPID::kSPECIES],ETRDPIDMethod PIDmethod,Bool_t kNorm) const
ffb1ee30	117	{
e0de37e9	118	//
	119	// Calculate TRD likelihood values for the track based on dedx and
	120	// momentum values. The likelihoods are calculated by query the
	121	// reference data depending on the PID method selected
	122	//
	123	// Input parameter :
	124	// n - number of dedx slices/chamber
	125	// dedx - array of dedx slices organized layer wise
	126	// p - array of momentum measurements organized layer wise
	127	//
	128	// Return parameters
	129	// prob - probabilities allocated by TRD for particle specis
	130	// kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
	131	//
	132	// Return value
	133	// true if calculation success
bd58d4b9	134	//
	135	AliDebug(3,Form(" Response for PID method: %d",PIDmethod));
	136
ffb1ee30	137
db0e2c5f	138	if(!fkPIDResponseObject){
	139	AliWarning("Missing reference data. PID calculation not possible.");
	140	return kFALSE;
	141	}
ffb1ee30	142
db0e2c5f	143	for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
	144	Double_t prLayer[AliPID::kSPECIES];
	145	Double_t dE[10], s(0.);
	146	for(Int_t il(kNlayer); il--;){
	147	memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
bd58d4b9	148	if(!CookdEdx(n, &dedx[il*n], &dE[0],PIDmethod)) continue;
ffb1ee30	149
db0e2c5f	150	s=0.;
	151	Bool_t filled=kTRUE;
	152	for(Int_t is(AliPID::kSPECIES); is--;){
bd58d4b9	153	if((PIDmethod==kLQ2D)&&(!(is==0\|\|is==2)))continue;
bd58d4b9	154	if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], &dE[0],PIDmethod);
f2762b1c	155	AliDebug(3, Form("Probability for Species %d in Layer %d: %e", is, il, prLayer[is]));
db0e2c5f	156	if(prLayer[is]<1.e-30){
	157	AliDebug(2, Form("Null for species %d species prob layer[%d].",is,il));
	158	filled=kFALSE;
	159	break;
	160	}
	161	s+=prLayer[is];
	162	}
	163	if(!filled){
	164	continue;
	165	}
	166	if(s<1.e-30){
	167	AliDebug(2, Form("Null all species prob layer[%d].", il));
	168	continue;
	169	}
	170	for(Int_t is(AliPID::kSPECIES); is--;){
	171	if(kNorm) prLayer[is] /= s;
	172	prob[is] *= prLayer[is];
	173	}
e0de37e9	174	}
db0e2c5f	175	if(!kNorm) return kTRUE;
	176
	177	s=0.;
	178	for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
e0de37e9	179	if(s<1.e-30){
db0e2c5f	180	AliDebug(2, "Null total prob.");
db0e2c5f	181	return kFALSE;
e0de37e9	182	}
db0e2c5f	183	for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
db0e2c5f	184	return kTRUE;
ffb1ee30	185	}
ffb1ee30	186
ffb1ee30	187	//____________________________________________________________
bd58d4b9	188	Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t *dEdx,ETRDPIDMethod PIDmethod) const {
e0de37e9	189	//
	190	// Get the non-normalized probability for a certain particle species as coming
	191	// from the reference histogram
	192	// Interpolation between momentum bins
	193	//
	194	AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
db0e2c5f	195
51a0ce25	196	Double_t probLayer = 0.;
db0e2c5f	197	Float_t pLower, pUpper;
db0e2c5f	198
bd58d4b9	199	switch(PIDmethod){
db0e2c5f	200	case kNN: // NN
	201	break;
	202	case kLQ2D: // 2D LQ
	203	{
f2762b1c	204	if(species==0\|\|species==2){ // references only for electrons and pions
	205	Double_t error;
	206	Double_t point[kNslicesLQ2D];
	207	for(Int_t idim=0;idim<kNslicesLQ2D;idim++){point[idim]=dEdx[idim];}
db0e2c5f	208
f2762b1c	209	AliTRDTKDInterpolator refLower = dynamic_cast<AliTRDTKDInterpolator>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ2D));
db0e2c5f	210
f2762b1c	211	if(refLower){
	212	refLower->Eval(point,probLayer,error);
	213	}
	214	else {
	215	AliError("No references available");
	216	}
	217	AliDebug(2,Form("Eval 2D Q0 %f Q1 %f P %e Err %e",point[0],point[1],probLayer,error));
db0e2c5f	218	}
	219	}
	220	break;
	221	case kLQ1D: // 1D LQ
	222	{
	223	TH1 refUpper = dynamic_cast<TH1 >(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ1D)),
	224	refLower = dynamic_cast<TH1 >(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ1D));
	225	// Do Interpolation exept for underflow and overflow
	226	if(refLower && refUpper){
	227	Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
	228	Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
	229
	230	probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
	231	} else if(refLower){
	232	// underflow
	233	probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
	234	} else if(refUpper){
	235	// overflow
	236	probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
	237	} else {
	238	AliError("No references available");
	239	}
f2762b1c	240	AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
db0e2c5f	241	}
	242	break;
	243	default:
	244	break;
e0de37e9	245	}
51a0ce25	246	return probLayer;
ffb1ee30	247	}
	248
	249	//____________________________________________________________
	250	void AliTRDPIDResponse::SetOwner(){
e0de37e9	251	//
	252	// Make Deep Copy of the Reference Histograms
	253	//
db0e2c5f	254	if(!fkPIDResponseObject \|\| IsOwner()) return;
	255	const AliTRDPIDResponseObject *tmp = fkPIDResponseObject;
	256	fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(tmp->Clone());
	257	SetBit(kIsOwner, kTRUE);
ffb1ee30	258	}
	259
	260	//____________________________________________________________
bd58d4b9	261	Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out,ETRDPIDMethod PIDmethod) const
ffb1ee30	262	{
f2762b1c	263	//
	264	// Recalculate dE/dx
	265	//
bd58d4b9	266	switch(PIDmethod){
e0de37e9	267	case kNN: // NN
db0e2c5f	268	break;
	269	case kLQ2D: // 2D LQ
	270	out[0]=0;
	271	out[1]=0;
	272	for(Int_t islice = 0; islice < nSlice; islice++){
	273	if(islice<fkPIDResponseObject->GetNSlicesQ0())out[0]+= in[islice];
	274	else out[1]+= in[islice];
	275	}
f2762b1c	276	if(out[0] < 1e-6) return kFALSE;
f2762b1c	277	AliDebug(3,Form("CookdEdx Q0 %f Q1 %f",out[0],out[1]));
db0e2c5f	278	break;
	279	case kLQ1D: // 1D LQ
	280	out[0]= 0.;
	281	for(Int_t islice = 0; islice < nSlice; islice++)
	282	if(in[islice] > 0) out[0] += in[islice] * fGainNormalisationFactor; // Protect against negative values for slices having no dE/dx information
	283	if(out[0] < 1e-6) return kFALSE;
f2762b1c	284	AliDebug(3,Form("CookdEdx dEdx %f",out[0]));
db0e2c5f	285	break;
db0e2c5f	286
e0de37e9	287	default:
db0e2c5f	288	return kFALSE;
e0de37e9	289	}
e0de37e9	290	return kTRUE;
ffb1ee30	291	}
ffb1ee30	292
ea235c90	293	//____________________________________________________________
bd58d4b9	294	Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level,Double_t centrality,ETRDPIDMethod PIDmethod) const {
db0e2c5f	295	//
	296	// Check whether particle is identified as electron assuming a certain electron efficiency level
	297	// Only electron and pion hypothesis is taken into account
	298	//
	299	// Inputs:
	300	// Number of tracklets
	301	// Likelihood values
	302	// Momentum
	303	// Electron efficiency level
	304	//
	305	// If the function fails when the params are not accessible, the function returns true
	306	//
	307	if(!fkPIDResponseObject){
ea235c90	308	AliError("No PID Param object available");
	309	return kTRUE;
	310	}
	311	Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]);
ce487a7f	312	Double_t params[4];
bd58d4b9	313	if(!fkPIDResponseObject->GetThresholdParameters(nTracklets, level, params,centrality,PIDmethod)){
ea235c90	314	AliError("No Params found for the given configuration");
	315	return kTRUE;
	316	}
ce487a7f	317	Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
ea235c90	318	if(probEle > TMath::Max(TMath::Min(threshold, 0.99), 0.2)) return kTRUE; // truncate the threshold upperwards to 0.999 and lowerwards to 0.2 and exclude unphysical values
	319	return kFALSE;
	320	}
	321
f2762b1c	322	//____________________________________________________________
	323	Bool_t AliTRDPIDResponse::SetPIDResponseObject(const AliTRDPIDResponseObject * obj){
	324
	325	fkPIDResponseObject = obj;
	326	if((AliLog::GetDebugLevel("",IsA()->GetName()))>0)fkPIDResponseObject->Print("");
	327	return kTRUE;
	328	}