Update master to aliroot
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliTRDPIDResponse.cxx
index edff1a1..cd2b72b 100644 (file)
 #include <TSystem.h>
 
 #include "AliLog.h"
+#include "AliVTrack.h"
 
-#include "AliTRDPIDParams.h"
-#include "AliTRDPIDReference.h"
+#include "AliTRDPIDResponseObject.h"
 #include "AliTRDPIDResponse.h"
+//#include "AliTRDTKDInterpolator.h"
+#include "AliTRDNDFast.h"
+#include "AliTRDdEdxParams.h"
 
 ClassImp(AliTRDPIDResponse)
 
 //____________________________________________________________
 AliTRDPIDResponse::AliTRDPIDResponse():
   TObject()
-  ,fkPIDReference(NULL)
-  ,fkPIDParams(NULL)
+  ,fkPIDResponseObject(NULL)
+  ,fkTRDdEdxParams(NULL)
   ,fGainNormalisationFactor(1.)
-  ,fPIDmethod(kLQ1D)
 {
   //
   // Default constructor
@@ -58,10 +61,9 @@ AliTRDPIDResponse::AliTRDPIDResponse():
 //____________________________________________________________
 AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
   TObject(ref)
-  ,fkPIDReference(ref.fkPIDReference)
-  ,fkPIDParams(ref.fkPIDParams)
+  ,fkPIDResponseObject(NULL)
+  ,fkTRDdEdxParams(NULL)
   ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
-  ,fPIDmethod(ref.fPIDmethod)
 {
   //
   // Copy constructor
@@ -78,10 +80,9 @@ AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
   // Make copy
   TObject::operator=(ref);
   fGainNormalisationFactor = ref.fGainNormalisationFactor;
-  fkPIDReference = ref.fkPIDReference;
-  fkPIDParams = ref.fkPIDParams;
-  fPIDmethod = ref.fPIDmethod;
-  
+  fkPIDResponseObject = ref.fkPIDResponseObject;
+  fkTRDdEdxParams = ref.fkTRDdEdxParams;
+
   return *this;
 }
 
@@ -90,11 +91,14 @@ AliTRDPIDResponse::~AliTRDPIDResponse(){
   //
   // Destructor
   //
-  if(IsOwner()) delete fkPIDReference;
+  if(IsOwner()) {
+    delete fkPIDResponseObject;
+    delete fkTRDdEdxParams;
+  }
 }
 
 //____________________________________________________________
-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
   //
@@ -110,16 +114,180 @@ Bool_t AliTRDPIDResponse::Load(const Char_t * filename, const Char_t *refName){
   }
   
   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;
 }
 
+
+
+//____________________________________________________________
+Double_t AliTRDPIDResponse::GetNumberOfSigmas(const AliVTrack *track, AliPID::EParticleType type) const
+{
+  //
+  //calculate the TRD nSigma
+  //
+
+  const Double_t badval = -9999;
+  Double_t info[5]; for(int i=0; i<5; i++){info[i]=badval;}
+
+  const Double_t delta = GetSignalDelta(track, type, kFALSE, info);
+
+  const Double_t mean = info[0];
+  const Double_t res = info[1];
+  if(res<0){
+    return badval;
+  }
+
+  const Double_t sigma = mean*res;
+  const Double_t eps = 1e-12;
+  return delta/(sigma + eps);
+}
+
+//____________________________________________________________
+Double_t AliTRDPIDResponse::GetSignalDelta( const AliVTrack* track, AliPID::EParticleType type, Bool_t ratio/*=kFALSE*/, Double_t *info/*=0x0*/) const
+{
+  //
+  //calculate the TRD signal difference w.r.t. the expected
+  //output other information in info[]
+  //
+
+  const Double_t badval = -9999;
+
+  if(!track){
+    return badval;
+  }
+
+  Double_t pTRD = -999;
+  for(Int_t ich=0; ich<6; ich++){
+    pTRD = track->GetTRDmomentum(ich);
+    if(pTRD>0)
+      break;
+  }
+  if(pTRD<0){
+    return badval;
+  }
+
+  if(!fkTRDdEdxParams){
+    AliError("fkTRDdEdxParams null");
+    return -99999;
+  }
+
+  const Double_t nch = track->GetTRDNchamberdEdx();
+  const Double_t ncls = track->GetTRDNclusterdEdx();
+
+  const TVectorF meanvec = fkTRDdEdxParams->GetMeanParameter(type, nch, ncls);
+  if(meanvec.GetNrows()==0){
+    return badval;
+  }
+
+  const TVectorF resvec  = fkTRDdEdxParams->GetSigmaParameter(type, nch, ncls);
+  if(resvec.GetNrows()==0){
+    return badval;
+  }
+
+  const Float_t *meanpar = meanvec.GetMatrixArray();
+  const Float_t *respar  = resvec.GetMatrixArray();
+
+  //============================================================================================<<<<<<<<<<<<<
+
+  const Double_t bg = pTRD/AliPID::ParticleMass(type);
+  const Double_t expsig = MeandEdxTR(&bg, meanpar);
+
+  if(info){
+    info[0]= expsig;
+    info[1]= ResolutiondEdxTR(&ncls, respar);
+  }
+
+  const Double_t eps = 1e-10;
+
+  if(ratio){
+    return track->GetTRDsignal()/(expsig +  eps);
+  }
+  else{
+    return track->GetTRDsignal() - expsig;
+  }
+}
+
+
+Double_t AliTRDPIDResponse::ResolutiondEdxTR(const Double_t * xx,  const Float_t * par)
+{
+  //
+  //resolution 
+  //npar=3
+  //
+
+  const Double_t ncls = xx[0];
+
+  return par[0]+par[1]*TMath::Power(ncls, par[2]);
+}
+
+Double_t AliTRDPIDResponse::MeandEdxTR(const Double_t * xx,  const Float_t * pin)
+{
+  //
+  //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
+  //npar = 8 = 3+5
+  //
+
+  Float_t ptr[4]={0};
+  for(int ii=0; ii<3; ii++){
+    ptr[ii+1]=pin[ii];
+  }
+  return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
+}
+
+Double_t AliTRDPIDResponse::MeanTR(const Double_t * xx,  const Float_t * par)
+{
+  //
+  //LOGISTIC parametrization for TR, in unit of MIP
+  //npar = 4
+  //
+                                                                                                                                                                                                                                                        
+  const Double_t bg = xx[0];
+  const Double_t gamma = sqrt(1+bg*bg);
+
+  const Double_t p0 = TMath::Abs(par[1]);
+  const Double_t p1 = TMath::Abs(par[2]);
+  const Double_t p2 = TMath::Abs(par[3]);
+
+  const Double_t zz = TMath::Log(gamma);
+  const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
+
+  return par[0]+tryield;
+}
+
+Double_t AliTRDPIDResponse::MeandEdx(const Double_t * xx,  const Float_t * par)
+{
+  //
+  //ALEPH parametrization for dEdx
+  //npar = 5
+  //
+
+  const Double_t bg = xx[0];
+  const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
+
+  const Double_t p0 = TMath::Abs(par[0]);
+  const Double_t p1 = TMath::Abs(par[1]);
+
+  const Double_t p2 = TMath::Abs(par[2]);
+
+  const Double_t p3 = TMath::Abs(par[3]);
+  const Double_t p4 = TMath::Abs(par[4]);
+
+  const Double_t aa = TMath::Power(beta, p3);
+
+  const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
+
+  return (p1-aa-bb)*p0/aa;
+
+}
+
+
 //____________________________________________________________
-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
+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 
@@ -136,77 +304,157 @@ 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(!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 0;
+    }
 
-    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.);
+    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.;
+        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];
+       }
+       ntrackletsPID++;
     }
+    if(!kNorm) return ntrackletsPID;
+
+    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;
+       AliDebug(2, "Null total prob.");
+       return 0;
     }
-    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;
+    for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
+    return ntrackletsPID;
 }
 
 //____________________________________________________________
-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));
+
+  Float_t pLower, pUpper;
+       
+  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->GetBinContent(refLower->GetXaxis()->FindBin(dEdx));
-    Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx));
-  
-    probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
+      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->GetBinContent(refLower->GetXaxis()->FindBin(dEdx));
+     // underflow
+      probLayer = refLower->Eval(dEdx);
   } else if(refUpper){
-    // overflow
-    probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx));
+      // overflow
+      probLayer = refUpper->Eval(dEdx);
   } else {
-    AliError("No references available");
+      AliError("No references available");
   }
-  AliDebug(1, Form("Probability %f", probLayer));
+  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
+      {
+         if(species==0||species==2){ // references only for electrons and pions
+             Double_t error = 0.;
+             Double_t point[kNslicesLQ2D];
+             for(Int_t idim=0;idim<kNslicesLQ2D;idim++){point[idim]=dEdx[idim];}
+
+             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 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));
+         }
+      }
+      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;
+      */
+
 }
 
 //____________________________________________________________
@@ -214,56 +462,71 @@ void AliTRDPIDResponse::SetOwner(){
   //
   // 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(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;
+  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
+      out[0]*=1./Double_t(nSlice);
+      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]);
   Double_t params[4];
-  if(!fkPIDParams->GetThresholdParameters(nTracklets, level, params)){
+  if(!fkPIDResponseObject->GetThresholdParameters(nTracklets, level, params,centrality,PIDmethod)){
     AliError("No Params found for the given configuration");
     return kTRUE;
   }
@@ -272,3 +535,18 @@ Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t
   return kFALSE;
 }
 
+//____________________________________________________________
+Bool_t AliTRDPIDResponse::SetPIDResponseObject(const AliTRDPIDResponseObject * obj){
+
+    fkPIDResponseObject = obj;
+    if((AliLog::GetDebugLevel("",IsA()->GetName()))>0)fkPIDResponseObject->Print("");
+    return kTRUE;
+}
+
+
+//____________________________________________________________    
+Bool_t AliTRDPIDResponse::SetdEdxParams(const AliTRDdEdxParams * par)
+{
+  fkTRDdEdxParams = par;
+  return kTRUE;
+}