- fast TRD PID implementation with muon constraints (Daniel Lohner)
authorpachmay <Yvonne.Chiara.Pachmayer@cern.ch>
Thu, 13 Mar 2014 13:20:11 +0000 (14:20 +0100)
committershahoian <ruben.shahoyan@cern.ch>
Thu, 13 Mar 2014 15:54:48 +0000 (16:54 +0100)
STEER/CMakelibSTEERBase.pkg
STEER/STEERBase/AliPIDResponse.cxx
STEER/STEERBase/AliPIDResponse.h
STEER/STEERBase/AliTRDNDFast.cxx [new file with mode: 0644]
STEER/STEERBase/AliTRDNDFast.h [new file with mode: 0644]
STEER/STEERBase/AliTRDPIDParams.cxx
STEER/STEERBase/AliTRDPIDResponse.cxx
STEER/STEERBase/AliTRDPIDResponse.h
STEER/STEERBase/AliTRDPIDResponseObject.cxx
STEER/STEERBaseLinkDef.h

index 8bcf4b5..f1980f8 100644 (file)
@@ -63,6 +63,7 @@ set ( SRCS
     STEERBase/AliTRDPIDReference.cxx 
     STEERBase/AliTRDPIDParams.cxx 
     STEERBase/AliTRDPIDResponseObject.cxx
+    STEERBase/AliTRDNDFast.cxx
     STEERBase/AliTRDTKDInterpolator.cxx
     STEERBase/AliPIDResponse.cxx 
     STEERBase/AliITSPIDResponse.cxx 
index 38b3e17..9a44642 100644 (file)
@@ -1643,15 +1643,23 @@ void AliPIDResponse::InitializeHMPIDResponse(){
 }
 
 //______________________________________________________________________________
-Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel,Double_t centrality,AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const {
+Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack,Double_t efficiencyLevel,Double_t centrality,AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const {
+    // old function for compatibility
+    Int_t ntracklets=0;
+    return IdentifiedAsElectronTRD(vtrack,ntracklets,efficiencyLevel,centrality,PIDmethod);
+}
+
+//______________________________________________________________________________
+Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Int_t &ntracklets,Double_t efficiencyLevel,Double_t centrality,AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const {
   //
   // Check whether track is identified as electron under a given electron efficiency hypothesis
     //
+    // ntracklets is the number of tracklets that has been used to calculate the PID signal
 
   Double_t probs[AliPID::kSPECIES];
-  ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs,PIDmethod);
 
-  Int_t ntracklets = vtrack->GetTRDntrackletsPID();
+  ntracklets =CalculateTRDResponse(vtrack,probs,PIDmethod);
+
   // Take mean of the TRD momenta in the given tracklets
   Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
   Int_t nmomenta = 0;
@@ -2339,8 +2347,31 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTOFProbability  (const A
   
   return kDetPidOk;
 }
+
+Int_t AliPIDResponse::CalculateTRDResponse(const AliVTrack *track,Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const
+{
+    // new function for backward compatibility
+    // returns number of tracklets PID
+
+    UInt_t TRDslicesForPID[2];
+    SetTRDSlices(TRDslicesForPID,PIDmethod);
+
+    Float_t mom[6]={0.};
+    Double_t dedx[48]={0.};  // Allocate space for the maximum number of TRD slices
+    Int_t nslices = TRDslicesForPID[1] - TRDslicesForPID[0] + 1;
+    AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d",  TRDslicesForPID[0], TRDslicesForPID[1], nslices));
+    for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
+       mom[ilayer] = track->GetTRDmomentum(ilayer);
+       for(UInt_t islice = TRDslicesForPID[0]; islice <= TRDslicesForPID[1]; islice++){
+           dedx[ilayer*nslices+islice-TRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
+       }
+    }
+
+    return fTRDResponse.GetResponse(nslices, dedx, mom, p,PIDmethod);
+
+}
 //______________________________________________________________________________
-AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod/*=AliTRDPIDResponse::kLQ1D*/) const
+AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const
 {
   //
   // Compute PID probabilities for the TRD
@@ -2352,21 +2383,8 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTRDProbability  (const A
   const EDetPidStatus pidStatus=GetTRDPIDStatus(track);
   if (pidStatus!=kDetPidOk) return pidStatus;
 
-  UInt_t TRDslicesForPID[2];
-  SetTRDSlices(TRDslicesForPID,PIDmethod);
-  
-  Float_t mom[6]={0.};
-  Double_t dedx[48]={0.};  // Allocate space for the maximum number of TRD slices
-  Int_t nslices = TRDslicesForPID[1] - TRDslicesForPID[0] + 1;
-  AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d",  TRDslicesForPID[0], TRDslicesForPID[1], nslices));
-  for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
-    mom[ilayer] = track->GetTRDmomentum(ilayer);
-    for(UInt_t islice = TRDslicesForPID[0]; islice <= TRDslicesForPID[1]; islice++){
-      dedx[ilayer*nslices+islice-TRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
-    }
-  }
-  
-  fTRDResponse.GetResponse(nslices, dedx, mom, p,PIDmethod);
+  CalculateTRDResponse(track,p,PIDmethod);
+
   return kDetPidOk;
 }
 
index ff24386..ac7ae84 100644 (file)
@@ -102,6 +102,8 @@ public:
   virtual Float_t NumberOfSigmasEMCAL(const AliVParticle *track, AliPID::EParticleType type) const;
 
   Bool_t IdentifiedAsElectronTRD(const AliVTrack *track, Double_t efficiencyLevel,Double_t centrality=-1,AliTRDPIDResponse::ETRDPIDMethod PIDmethod=AliTRDPIDResponse::kLQ1D) const;
+  Bool_t IdentifiedAsElectronTRD(const AliVTrack *track, Int_t &ntracklets, Double_t efficiencyLevel,Double_t centrality=-1,AliTRDPIDResponse::ETRDPIDMethod PIDmethod=AliTRDPIDResponse::kLQ1D) const;
+  
 
   // Signal delta
   EDetPidStatus GetSignalDelta(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type, Double_t &val, Bool_t ratio=kFALSE) const;
@@ -195,7 +197,8 @@ protected:
   //unbuffered PID calculation
   virtual Float_t GetNumberOfSigmasTOFold  (const AliVParticle */*track*/, AliPID::EParticleType /*type*/) const {return 0;}
   virtual Float_t GetSignalDeltaTOFold(const AliVParticle */*track*/, AliPID::EParticleType /*type*/, Bool_t /*ratio*/=kFALSE) const {return -9999.;}
-  
+
+  Int_t CalculateTRDResponse(const AliVTrack *track, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const;
   EDetPidStatus GetComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod=AliTRDPIDResponse::kLQ1D) const;
   EDetPidStatus GetTOFPIDStatus(const AliVTrack *track) const;
 
diff --git a/STEER/STEERBase/AliTRDNDFast.cxx b/STEER/STEERBase/AliTRDNDFast.cxx
new file mode 100644 (file)
index 0000000..99ee898
--- /dev/null
@@ -0,0 +1,341 @@
+// Author: Daniel.Lohner@cern.ch
+
+#include "AliTRDNDFast.h"
+#include "TCanvas.h"
+
+extern Double_t langaufunc(Double_t *x, Double_t *par) {
+
+    // needs protection, parameter [3]>0!!!!!
+    // needs protection, parameter [4]>0!!!!!
+
+    //Fit parameters:
+    //par[0]=Width (scale) parameter of Landau density
+    //par[1]=Most Probable (MP, location) parameter of Landau density
+    //par[2]=Total area (integral -inf to inf, normalization constant)
+    //par[3]=Width (sigma) of convoluted Gaussian function
+    //par[4]=Exponential Slope Parameter
+    //
+    //In the Landau distribution (represented by the CERNLIB approximation),
+    //the maximum is located at x=-0.22278298 with the location parameter=0.
+    //This shift is corrected within this function, so that the actual
+    //maximum is identical to the MP parameter.
+
+    // Numeric constants
+    Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
+    Double_t mpshift  = -0.22278298;       // Landau maximum location
+
+    // Control constants
+    Double_t np = 100.0;      // number of convolution stpdf
+    Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
+
+    // Variables
+    Double_t xx;
+    Double_t mpc;
+    Double_t fland;
+    Double_t sum = 0.0;
+    Double_t xlow,xupp;
+    Double_t step;
+    Double_t i;
+
+    // MP shift correction
+    mpc = par[1] - mpshift * par[0];
+
+    // Range of convolution integral
+    xlow = x[0] - sc * par[3];
+    xupp = x[0] + sc * par[3];
+
+    if(xlow<0)xlow=0;
+    if(xupp<xlow)cout<<"ERROR: Convolution around Negative MPV"<<endl;
+
+    step = (xupp-xlow) / np;
+
+    // Convolution integral of Landau and Gaussian by sum
+    for(i=1.0; i<=np/2; i++) {
+       xx = xlow + (i-.5) * step;
+       fland = TMath::Landau(xx,mpc,par[0])*TMath::Exp(-par[4]*xx) / par[0];   // mpc taken out
+       sum += fland * TMath::Gaus(x[0],xx,par[3]);
+
+       xx = xupp - (i-.5) * step;
+       fland = TMath::Landau(xx,mpc,par[0])*TMath::Exp(-par[4]*xx) / par[0];   // mpc taken out
+       sum += fland * TMath::Gaus(x[0],xx,par[3]);
+    }
+
+    return (par[2] * step * sum * invsq2pi / par[3]);
+}
+
+
+
+ClassImp(AliTRDNDFast);
+
+AliTRDNDFast::AliTRDNDFast(): TObject(),
+    fNDim(0),
+    fTitle(""),
+    fFunc(NULL),
+    fHistos(NULL),
+    fPars()
+{
+    // default constructor
+}
+
+AliTRDNDFast::AliTRDNDFast(const char *name,Int_t ndim,Int_t nbins,Double_t xlow,Double_t xup): TObject(),
+fNDim(ndim),
+fTitle(name),
+fFunc(NULL),
+fHistos(NULL),
+fPars()
+{
+    Init();
+    for(Int_t idim=0;idim<fNDim;idim++){
+       fHistos[idim]=new TH1F(Form("%s_axis_%d",fTitle.Data(),idim),Form("%s_axis_%d",fTitle.Data(),idim),nbins,xlow,xup);
+       fHistos[idim]->SetDirectory(0);
+    }
+}
+
+AliTRDNDFast::AliTRDNDFast(const AliTRDNDFast &ref) : TObject(ref),
+fNDim(ref.fNDim),
+fTitle(ref.fTitle),
+fFunc(NULL),
+fHistos(NULL),
+fPars()
+{
+    Cleanup();
+    Init();
+    for(Int_t idim=0;idim<fNDim;idim++){
+       fHistos[idim]=(TH1F*)ref.fHistos[idim]->Clone(Form("%s_axis_%d",GetName(),idim));
+       fHistos[idim]->SetDirectory(0);
+       for(Int_t ipar=0;ipar<kNpar;ipar++)fPars[idim][ipar]=ref.fPars[idim][ipar];
+    }
+}
+
+AliTRDNDFast &AliTRDNDFast::operator=(const AliTRDNDFast &ref){
+    //
+    // Assignment operator
+    //
+    if(this != &ref){
+       TObject::operator=(ref);
+       fNDim=ref.fNDim;
+       fTitle=ref.fTitle;
+       fFunc = ref.fFunc;
+       fHistos=ref.fHistos;
+       for(Int_t idim=0;idim<fNDim;idim++){
+           for(Int_t ipar=0;ipar<kNpar;ipar++)fPars[idim][ipar]=ref.fPars[idim][ipar];
+       }
+    }
+    return *this;
+}
+
+AliTRDNDFast::~AliTRDNDFast(){
+    Cleanup();
+
+}
+
+void AliTRDNDFast::Init(){
+
+    for(Int_t ipar=0;ipar<kNpar;ipar++)fPars[ipar].Set(fNDim);
+    fFunc=new TF1*[fNDim];
+    fHistos=new TH1F*[fNDim];
+    for(Int_t idim=0;idim<fNDim;idim++){
+        fHistos[idim]=NULL;
+       for(Int_t ipar=0;ipar<kNpar;ipar++)fPars[ipar].SetAt(0,idim);
+       fFunc[idim]=NULL;
+    }
+}
+
+void AliTRDNDFast::Cleanup(){
+    if(fHistos){
+       for(Int_t idim=0;idim<fNDim;idim++){
+           if(fHistos[idim]){
+               delete fHistos[idim];
+               fHistos[idim]=NULL;
+           }
+       }
+       delete[] fHistos;
+       fHistos=NULL;
+    }
+    if(fPars){
+       for(Int_t ipar=0;ipar<kNpar;ipar++){
+           fPars[ipar].Reset();
+       }
+    }
+    if(fFunc){
+       for(Int_t idim=0;idim<fNDim;idim++){
+           if(fFunc[idim]){
+               delete fFunc[idim];
+               fFunc[idim]=NULL;
+           }
+       }
+       delete[] fFunc;
+       fFunc=NULL;
+    }
+}
+
+void AliTRDNDFast::PrintPars(){
+    for(Int_t idim=0;idim<fNDim;idim++){
+       for(Int_t ipar=0;ipar<kNpar;ipar++)cout<<idim<<" "<<ipar<<" "<<GetParam(idim,ipar)<<endl;
+    }
+}
+
+void AliTRDNDFast::ScaleLangauFun(TF1 *func,Double_t mpv){
+
+    Double_t start[kNpar],low[kNpar],high[kNpar];
+    for(Int_t ii=0;ii<kNpar;ii++){
+       start[ii]=func->GetParameter(ii);
+       func->GetParLimits(ii,low[ii],high[ii]);
+    }
+
+    Double_t scalefactor=mpv/start[1];
+
+    for(Int_t ii=0;ii<kNpar;ii++){
+       Double_t scale=1;
+       if(ii==0||ii==1||ii==3)scale=scalefactor;
+       if(ii==4)scale=1./scalefactor;
+       start[ii]*=scale;
+       low[ii]*=scale;
+       high[ii]*=scale;
+       func->SetParLimits(ii, low[ii], high[ii]);
+    }
+    func->SetParameters(start);
+}
+
+//---------------------------------------------------------------
+//---------------------------------------------------------------
+TF1 *AliTRDNDFast::GetLangauFun(TString funcname,Double_t range[2],Double_t scalefactor){
+
+    Double_t start[kNpar] = {120, 1000, 1, 100, 1e-5};
+    Double_t low[kNpar] = {10, 300, 0.01, 1, 0.};
+    Double_t high[kNpar] = {1000, 3000, 100, 1000, 1.};
+
+    TF1 *hlandauE=new TF1(funcname.Data(),langaufunc,0,range[1],kNpar);
+    hlandauE->SetParameters(start);
+    hlandauE->SetParNames("Width","MP","Area","GSigma","delta");
+    for (int i=0; i<kNpar; i++) {
+       hlandauE->SetParLimits(i, low[i], high[i]);
+    }
+
+    if(scalefactor!=1){ScaleLangauFun(hlandauE,scalefactor*start[1]);}
+
+    return hlandauE;
+}
+
+TF1 *AliTRDNDFast::FitLandau(TString name,TH1F *htemp,Double_t range[2],TString option){
+
+    TF1 *fitlandau1D=GetLangauFun(name,range);
+    TF1 fit("land","landau");
+    Double_t max=htemp->GetXaxis()->GetBinCenter(htemp->GetMaximumBin());
+    fit.SetParameter(1,max);
+    fit.SetParLimits(1,0,htemp->GetXaxis()->GetXmax());
+    fit.SetParameter(2,0.3*max); // MPV may never become negative!!!!!
+    htemp->Fit("land",option.Data(),"",range[0],range[1]);
+    ScaleLangauFun(fitlandau1D,fit.GetParameter(1));
+    htemp->Fit(fitlandau1D,option.Data(),"",range[0],range[1]); // range for references
+    return fitlandau1D;
+}
+
+void AliTRDNDFast::BuildHistos(){
+
+    for(Int_t idim=0;idim<fNDim;idim++){
+        fHistos[idim]->Reset();
+       // Fill Histo
+       Double_t pars[kNpar];
+       for(Int_t ipar=0;ipar<kNpar;ipar++)pars[ipar]=GetParam(idim,ipar);
+        // Also Fill overflow bins!!!
+       for(Int_t ii=0;ii<=fHistos[idim]->GetNbinsX()+1;ii++){
+           Double_t xval=fHistos[idim]->GetXaxis()->GetBinCenter(ii);
+           Double_t val=langaufunc(&xval,&pars[0]);
+           //Double_t val=fFunc[idim]->Eval(xval);
+           fHistos[idim]->SetBinContent(ii,val);
+       }
+       // Normalization
+       fHistos[idim]->Scale(1./fHistos[idim]->Integral(1,fHistos[idim]->GetNbinsX(),"width"));
+    }
+}
+
+void AliTRDNDFast::Build(Double_t **pars){
+    // pars[ndim][npar]
+    for(Int_t idim=0;idim<fNDim;idim++){
+       for(Int_t ipar=0;ipar<kNpar;ipar++){
+           fPars[ipar].SetAt(pars[idim][ipar],idim);
+       }
+    }
+    BuildHistos();
+}
+
+
+void AliTRDNDFast::Build(TH1F **hdEdx,TString path){
+
+    Double_t range[2];
+
+    TCanvas *CANV=new TCanvas("a","a");
+    if(fNDim==2)CANV->Divide(2,1);
+    if(fNDim==3)CANV->Divide(2,2);
+    if(fNDim>6)CANV->Divide(3,3);
+    // Fit and Extract Parameters
+    for(Int_t idim=0;idim<fNDim;idim++){
+       CANV->cd(idim+1);
+        gPad->SetLogy();
+       range[0]=hdEdx[idim]->GetXaxis()->GetXmin();
+       range[1]=hdEdx[idim]->GetXaxis()->GetXmax();
+       // Norm Histogram
+       hdEdx[idim]->Scale(1./hdEdx[idim]->Integral(1,hdEdx[idim]->GetNbinsX(),"width"));
+        // Fit Histogram
+       fFunc[idim]=FitLandau(Form("fit%d",idim),hdEdx[idim],range,"RMB0");
+       // Norm Landau
+       fFunc[idim]->SetParameter(2,fFunc[idim]->GetParameter(2)/fFunc[idim]->Integral(range[0],range[1]));
+       hdEdx[idim]->DrawCopy();
+        fFunc[idim]->DrawCopy("same");
+       // Set Pars
+       for(Int_t ipar=0;ipar<kNpar;ipar++)fPars[ipar].SetAt(fFunc[idim]->GetParameter(ipar),idim);
+    }
+    if(path!="")CANV->Print(Form("%s/%s_Build.pdf",path.Data(),fTitle.Data()));
+    delete CANV;
+
+    BuildHistos();
+}
+
+Double_t AliTRDNDFast::Eval(Double_t *point) const{
+    Double_t val=1;
+    for(Int_t idim=0;idim<fNDim;idim++){
+       Int_t bin=fHistos[idim]->GetXaxis()->FindBin(point[idim]);
+       val*=fHistos[idim]->GetBinContent(bin);
+    }
+    return val;
+}
+
+void AliTRDNDFast::Random(Double_t *point) const{
+    for(Int_t idim=0;idim<fNDim;idim++){
+       point[idim]=fHistos[idim]->GetRandom();
+    }
+}
+
+void AliTRDNDFast::Random(Double_t *point,AliTRDNDFast *nd0,AliTRDNDFast *nd1,Double_t w0,Double_t w1){
+    for(Int_t idim=0;idim<nd0->fNDim;idim++){
+       point[idim]=GetRandomInterpolation(nd0->fHistos[idim],nd1->fHistos[idim],w0,w1);
+    }
+}
+
+Int_t AliTRDNDFast::BinarySearchInterpolation(Int_t start,Int_t end,Double_t *a0,Double_t *a1,Double_t w0,Double_t w1,Double_t val){
+
+    if((end-start)==1)return start;
+    Int_t mid=Int_t((end+start)/2);
+    Double_t valmid=(w0*a0[mid]+w1*a1[mid])/(w0+w1);
+    if(val>=valmid)return BinarySearchInterpolation(mid,end,a0,a1,w0,w1,val);
+    return BinarySearchInterpolation(start,mid,a0,a1,w0,w1,val);
+}
+
+Double_t AliTRDNDFast::GetRandomInterpolation(TH1F *hist0,TH1F *hist1,Double_t w0,Double_t w1){
+
+    // Draw Random Variable
+    Double_t rand=gRandom->Rndm();
+
+    // Get Integral Arrays
+    Double_t *integral0=hist0->GetIntegral();
+    Double_t *integral1=hist1->GetIntegral();
+
+    Int_t nbinsX=hist0->GetNbinsX();
+
+    Int_t ibin=BinarySearchInterpolation(1,nbinsX+1,integral0,integral1,w0,w1,rand);
+    return hist0->GetXaxis()->GetBinCenter(ibin);
+}
+
+
+
diff --git a/STEER/STEERBase/AliTRDNDFast.h b/STEER/STEERBase/AliTRDNDFast.h
new file mode 100644 (file)
index 0000000..16d85ac
--- /dev/null
@@ -0,0 +1,68 @@
+// Author: Daniel.Lohner@cern.ch
+
+#ifndef ALIROOT_AliTRDNDFast
+#define ALIROOT_AliTRDNDFast
+
+#ifndef ROOT_TH1
+#include "TH1F.h"
+#endif
+#ifndef ROOT_TArrayF
+#include "TArrayF.h"
+#endif
+#ifndef ROOT_TF2
+#include "TF1.h"
+#endif
+#ifndef ROOT_TMath
+#include "TMath.h"
+#endif
+#ifndef ROOT_TRandom
+#include "TRandom.h"
+#endif
+
+using namespace std;
+
+extern Double_t langaufun(Double_t *x,Double_t *par);
+
+class AliTRDNDFast : public TObject {
+
+public:
+    static const Int_t kNpar = 5;
+
+    AliTRDNDFast();
+    AliTRDNDFast(const char *name,Int_t ndim,Int_t nbins,Double_t xlow,Double_t xup);
+    AliTRDNDFast(const AliTRDNDFast&);
+    AliTRDNDFast &operator=(const AliTRDNDFast &ref);
+    virtual ~AliTRDNDFast();
+    
+    TF1 *FitLandau(TString name,TH1F *htemp,Double_t range[2],TString option);
+
+    void Build(TH1F **hdEdx,TString path="");
+    void Build(Double_t **pars);
+    Double_t Eval(Double_t *point) const;
+    void Random(Double_t *point) const;
+    Int_t GetNDim(){return fNDim;};
+    Double_t GetParam(Int_t dim,Int_t par){if((dim>=0)&&(dim<fNDim)&&(par>=0)&&(par<kNpar)){return fPars[par].GetAt(dim);}else{return 0;}};
+    void PrintPars();
+    static void Random(Double_t *point,AliTRDNDFast *nd0,AliTRDNDFast *nd1,Double_t w0,Double_t w1);
+
+private:
+
+    void ScaleLangauFun(TF1 *func,Double_t mpv);
+    TF1 *GetLangauFun(TString funcname,Double_t range[2],Double_t scalefactor=1);
+    void BuildHistos();
+    void Init();
+    void Cleanup();
+    
+    static Int_t BinarySearchInterpolation(Int_t start,Int_t end,Double_t *a0,Double_t *a1,Double_t w0,Double_t w1,Double_t val);
+    static Double_t GetRandomInterpolation(TH1F *hist0,TH1F *hist1,Double_t w0,Double_t w1);
+
+    Int_t fNDim; // Dimensions
+    TString fTitle; //title
+    TF1 **fFunc; //! functions, do not store
+    TH1F **fHistos; //[fNDim] Histograms
+    TArrayF fPars[kNpar]; // parameters
+
+    ClassDef(AliTRDNDFast,1)  //Fast TRD ND class
+};
+
+#endif
index a421a48..7377bbd 100644 (file)
@@ -99,6 +99,10 @@ AliTRDPIDParams::AliTRDPIDCentrality *AliTRDPIDParams::FindCentrality(Double_t v
       break;
     }
   }
+  if(!tmp){
+      AliDebug(10,"Using default centrality class");
+      return FindCentrality(-1);
+  }
   return tmp;
 }
 
index 1b2e64a..0db79d9 100644 (file)
@@ -38,7 +38,8 @@
 
 #include "AliTRDPIDResponseObject.h"
 #include "AliTRDPIDResponse.h"
-#include "AliTRDTKDInterpolator.h"
+//#include "AliTRDTKDInterpolator.h"
+#include "AliTRDNDFast.h"
 
 ClassImp(AliTRDPIDResponse)
 
@@ -113,7 +114,7 @@ Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
 }
 
 //____________________________________________________________
-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
+Int_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 
@@ -130,27 +131,27 @@ Bool_t AliTRDPIDResponse::GetResponse(Int_t n, const Double_t * const dedx, cons
   //   kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
   // 
   // Return value
-  //   true if calculation success
+  //   number of tracklets used for PID, 0 if no PID
     //
     AliDebug(3,Form(" Response for PID method: %d",PIDmethod));
 
 
     if(!fkPIDResponseObject){
        AliDebug(3,"Missing reference data. PID calculation not possible.");
-       return kFALSE;
+       return 0;
     }
 
     for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
     Double_t prLayer[AliPID::kSPECIES];
     Double_t dE[10], s(0.);
+    Int_t ntrackletsPID=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.;
+        s=0.;
         Bool_t filled=kTRUE;
        for(Int_t is(AliPID::kSPECIES); is--;){
-           if((PIDmethod==kLQ2D)&&(!(is==0||is==2)))continue;
+           //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){
@@ -171,17 +172,18 @@ Bool_t AliTRDPIDResponse::GetResponse(Int_t n, const Double_t * const dedx, cons
            if(kNorm) prLayer[is] /= s;
            prob[is] *= prLayer[is];
        }
+       ntrackletsPID++;
     }
-    if(!kNorm) return kTRUE;
+    if(!kNorm) return ntrackletsPID;
 
     s=0.;
     for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
     if(s<1.e-30){
        AliDebug(2, "Null total prob.");
-       return kFALSE;
+       return 0;
     }
     for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
-    return kTRUE;
+    return ntrackletsPID;
 }
 
 //____________________________________________________________
@@ -194,10 +196,35 @@ Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t pl
   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
+  AliTRDNDFast *refUpper = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,PIDmethod)),
+      *refLower = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,PIDmethod));
+
+  // Do Interpolation exept for underflow and overflow
+  if(refLower && refUpper){
+      Double_t probLower = refLower->Eval(dEdx);
+      Double_t probUpper = refUpper->Eval(dEdx);
+
+      probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
+  } else if(refLower){
+     // underflow
+      probLayer = refLower->Eval(dEdx);
+  } else if(refUpper){
+      // overflow
+      probLayer = refUpper->Eval(dEdx);
+  } else {
+      AliError("No references available");
+  }
+  AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
+
+  return probLayer;
+
+/* old implementation
+
+switch(PIDmethod){
+case kNN: // NN
       break;
   case kLQ2D: // 2D LQ
       {
@@ -206,12 +233,21 @@ Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t pl
              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){
+             AliTRDTKDInterpolator *refUpper = dynamic_cast<AliTRDTKDInterpolator *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ2D)),
+             *refLower = dynamic_cast<AliTRDTKDInterpolator *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ2D));
+             // Do Interpolation exept for underflow and overflow
+             if(refLower && refUpper){
+                  Double_t probLower=0,probUpper=0;
+                 refLower->Eval(point,probLower,error);
+                  refUpper->Eval(point,probUpper,error);
+                 probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
+             } else if(refLower){
+                 // underflow
                  refLower->Eval(point,probLayer,error);
-             }
-             else {
+                 } else if(refUpper){
+                 // overflow
+                 refUpper->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));
@@ -242,8 +278,10 @@ Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t pl
       break;
   default:
       break;
-  }
-  return probLayer;
+      }
+      return probLayer;
+      */
+
 }
 
 //____________________________________________________________
@@ -270,9 +308,13 @@ Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Doub
       out[0]=0;
       out[1]=0;
       for(Int_t islice = 0; islice < nSlice; islice++){
+         if(in[islice]<=0){out[0]=0;out[1]=0;return kFALSE;}  // Require that all slices are filled
          if(islice<fkPIDResponseObject->GetNSlicesQ0())out[0]+= in[islice];
          else out[1]+= in[islice];
       }
+      // normalize signal to number of slices
+      out[0]*=1./Double_t(fkPIDResponseObject->GetNSlicesQ0());
+      out[1]*=1./Double_t(nSlice-fkPIDResponseObject->GetNSlicesQ0());
       if(out[0] < 1e-6) return kFALSE;
       AliDebug(3,Form("CookdEdx Q0 %f Q1 %f",out[0],out[1]));
       break;
@@ -280,6 +322,7 @@ Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Doub
       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
+      out[0]*=1./Double_t(nSlice);
       if(out[0] < 1e-6) return kFALSE;
       AliDebug(3,Form("CookdEdx dEdx %f",out[0]));
       break;
index cb85fa9..b9c1db6 100644 (file)
@@ -39,7 +39,7 @@ class AliTRDPIDResponse : public TObject {
     };
     enum ETRDPIDResponseDef {
        kNlayer = 6
-       ,kNPBins = 6
+       ,kNPBins = 6
     };
     enum ETRDPIDMethod {
        kNN   = 0,
@@ -51,15 +51,15 @@ class AliTRDPIDResponse : public TObject {
     };
     enum ETRDNslices {
        kNslicesLQ1D = 1,
-      kNslicesLQ2D = 2,
-      kNslicesNN = 7
+       kNslicesLQ2D = 2,
+       kNslicesNN = 7
     };
     AliTRDPIDResponse();
     AliTRDPIDResponse(const AliTRDPIDResponse &ref);
     AliTRDPIDResponse& operator=(const AliTRDPIDResponse &ref);
     ~AliTRDPIDResponse();
     
-    Bool_t    GetResponse(Int_t n, const Double_t * const dedx, const Float_t * const p, Double_t prob[AliPID::kSPECIES],ETRDPIDMethod PIDmethod=kLQ1D, Bool_t kNorm=kTRUE) const;
+    Int_t    GetResponse(Int_t n, const Double_t * const dedx, const Float_t * const p, Double_t prob[AliPID::kSPECIES],ETRDPIDMethod PIDmethod=kLQ1D, Bool_t kNorm=kTRUE) const;
     inline ETRDNslices  GetNumberOfSlices(ETRDPIDMethod PIDmethod=kLQ1D) const;
     
     Bool_t    IsOwner() const {return TestBit(kIsOwner);}
index a491650..b7f0445 100644 (file)
@@ -200,6 +200,7 @@ Bool_t AliTRDPIDResponseObject::GetThresholdParameters(Int_t ntracklets, Double_
     if(fPIDParams[method]){
        return fPIDParams[method]->GetThresholdParameters(ntracklets,efficiency,params,centrality);
     }
+    AliError("TRD Threshold Container does not exist");
     return kFALSE;
 }
 
index f519afd..ce1cb5d 100644 (file)
@@ -93,6 +93,7 @@
 #pragma link C++ class AliTRDPIDParams::AliTRDPIDCentrality+;
 /* #endif */
 #pragma link C++ class AliTRDPIDResponseObject+;
+#pragma link C++ class AliTRDNDFast+;
 #pragma link C++ class AliTRDTKDInterpolator+;
 #pragma link C++ class AliTRDTKDInterpolator::AliTRDTKDNodeInfo+;
 #pragma link C++ class AliITSPidParams+;