#include <TROOT.h>
#include <TString.h>
#include <TSystem.h>
-#include <TVectorT.h>
#include "AliLog.h"
-#include "AliTRDPIDReference.h"
+#include "AliTRDPIDResponseObject.h"
#include "AliTRDPIDResponse.h"
+#include "AliTRDTKDInterpolator.h"
ClassImp(AliTRDPIDResponse)
//____________________________________________________________
AliTRDPIDResponse::AliTRDPIDResponse():
TObject()
- ,fkPIDReference(NULL)
- ,fkPIDParams(NULL)
+ ,fkPIDResponseObject(NULL)
,fGainNormalisationFactor(1.)
- ,fPIDmethod(kLQ1D)
{
//
// Default constructor
//____________________________________________________________
AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
TObject(ref)
- ,fkPIDReference(ref.fkPIDReference)
- ,fkPIDParams(ref.fkPIDParams)
+ ,fkPIDResponseObject(NULL)
,fGainNormalisationFactor(ref.fGainNormalisationFactor)
- ,fPIDmethod(ref.fPIDmethod)
{
//
// Copy constructor
// Make copy
TObject::operator=(ref);
fGainNormalisationFactor = ref.fGainNormalisationFactor;
- fkPIDReference = ref.fkPIDReference;
- fkPIDParams = ref.fkPIDParams;
- fPIDmethod = ref.fPIDmethod;
+ fkPIDResponseObject = ref.fkPIDResponseObject;
return *this;
}
//
// Destructor
//
- if(IsOwner()) delete fkPIDReference;
+ if(IsOwner()) delete fkPIDResponseObject;
}
//____________________________________________________________
-Bool_t AliTRDPIDResponse::Load(const Char_t * filename, const Char_t *refName){
+Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
//
// Load References into the toolkit
//
}
gROOT->cd();
- fkPIDReference = dynamic_cast<const AliTRDPIDReference *>(in->Get(refName)->Clone());
+ 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", fkPIDReference->GetNumberOfMomentumBins()));
+ 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], Bool_t kNorm) const
+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
//
// Return value
// true if calculation success
- //
+ //
+ AliDebug(3,Form(" Response for PID method: %d",PIDmethod));
- if(!fkPIDReference){
- 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])) continue;
+ if(!fkPIDResponseObject){
+ AliDebug(3,"Missing reference data. PID calculation not possible.");
+ return kFALSE;
+ }
- s=0.;
- for(Int_t is(AliPID::kSPECIES); is--;){
- if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], dE[0]);
- AliDebug(3, Form("Probability for Species %d in Layer %d: %f", is, il, prLayer[is]));
- s+=prLayer[is];
+ 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, 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];
+ AliDebug(2, "Null total prob.");
+ return kFALSE;
}
- }
- 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;
+ 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) const {
+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));
- Float_t pLower, pUpper;
+
Double_t probLayer = 0.;
- TH1 *refUpper = dynamic_cast<TH1 *>(fkPIDReference->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper)),
- *refLower = dynamic_cast<TH1 *>(fkPIDReference->GetLowerReference((AliPID::EParticleType)species, plocal, pLower));
- // Do Interpolation exept for underflow and overflow
- if(refLower && refUpper){
- Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx));
- Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx));
-
- probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
- } else if(refLower){
- // underflow
- probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx));
- } else if(refUpper){
- // overflow
- probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx));
- } else {
- AliError("No references available");
+ 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;
}
- AliDebug(1, Form("Probability %f", probLayer));
return probLayer;
}
//
// Make Deep Copy of the Reference Histograms
//
- if(!fkPIDReference || IsOwner()) return;
- const AliTRDPIDReference *tmp = fkPIDReference;
- fkPIDReference = dynamic_cast<const AliTRDPIDReference *>(tmp->Clone());
- SetBit(kIsOwner, kTRUE);
+ 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) const
+Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out,ETRDPIDMethod PIDmethod) const
{
- //
- // Recalculate dE/dx
- //
- switch(fPIDmethod){
+ //
+ // Recalculate dE/dx
+ //
+ switch(PIDmethod){
case kNN: // NN
- break;
- case kLQ2D: // 2D LQ
- 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;
- break;
+ 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 kFALSE;
}
return kTRUE;
}
//____________________________________________________________
-Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level) 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(!fkPIDParams){
- AliError("No PID Param object available");
+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){
+ AliDebug(3,"No PID Param object available");
return kTRUE;
}
Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]);
- const TVectorD *params = GetParams(nTracklets, level);
- if(!params){
+ 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);
+ 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;
}
//____________________________________________________________
-const TVectorD* AliTRDPIDResponse::GetParams(Int_t ntracklets, Double_t level) const {
- //
- // returns the threshold for a given number of tracklets and a given efficiency level
- //tby definition the lower of step is given.
- //
- if(ntracklets > 6 || ntracklets <=0) return NULL;
- TObjArray * entry = dynamic_cast<TObjArray *>(fkPIDParams->At(ntracklets - 1));
- if(!entry) return NULL;
-
- TObjArray*cut = NULL;
- TVectorF *effLevel = NULL; const TVectorD *parameters = NULL;
- Float_t currentLower = 0.;
- TIter cutIter(entry);
- while((cut = dynamic_cast<TObjArray *>(cutIter()))){
- effLevel = static_cast<TVectorF *>(cut->At(0));
- if((*effLevel)[0] > currentLower && (*effLevel)[0] <= level){
- // New Lower entry found
- parameters = static_cast<const TVectorD *>(cut->At(1));
- currentLower = (*effLevel)[0];
- }
- }
- AliDebug(2, Form("Taking params for %d tracklets and %f electron Efficiency\n", ntracklets, currentLower));
+Bool_t AliTRDPIDResponse::SetPIDResponseObject(const AliTRDPIDResponseObject * obj){
- return parameters;
+ fkPIDResponseObject = obj;
+ if((AliLog::GetDebugLevel("",IsA()->GetName()))>0)fkPIDResponseObject->Print("");
+ return kTRUE;
}