Updates in order to enable the '2D' PID for the TRD developed by Daniel Lohner.

[u/mrichter/AliRoot.git] / STEER / STEERBase / AliTRDPIDResponse.cxx

diff --git a/STEER/STEERBase/AliTRDPIDResponse.cxx b/STEER/STEERBase/AliTRDPIDResponse.cxx
index 2f403369eb707acb19eb963ac8cac682e7671fb1..fe69bca7258e9e73930db1ae6f8749c8e8d0aae0 100644 (file)
--- a/STEER/STEERBase/AliTRDPIDResponse.cxx
+++ b/STEER/STEERBase/AliTRDPIDResponse.cxx
@@ -33,20 +33,21 @@
 #include <TROOT.h> 
 #include <TString.h>
 #include <TSystem.h>
 #include <TROOT.h> 
 #include <TString.h>
 #include <TSystem.h>

-#include <TVectorT.h>

 
 #include "AliLog.h"
 
 
 #include "AliLog.h"
 

-#include "AliTRDPIDReference.h"
+#include "AliTRDPIDResponseObject.h"
+//#include "AliTRDPIDParams.h"
+//#include "AliTRDPIDReference.h"

 #include "AliTRDPIDResponse.h"
 #include "AliTRDPIDResponse.h"

+#include "AliTRDTKDInterpolator.h"

 
 ClassImp(AliTRDPIDResponse)
 
 //____________________________________________________________
 AliTRDPIDResponse::AliTRDPIDResponse():
   TObject()
 
 ClassImp(AliTRDPIDResponse)
 
 //____________________________________________________________
 AliTRDPIDResponse::AliTRDPIDResponse():
   TObject()

-  ,fkPIDReference(NULL)
-  ,fkPIDParams(NULL)
+  ,fkPIDResponseObject(NULL)

   ,fGainNormalisationFactor(1.)
   ,fPIDmethod(kLQ1D)
 {
   ,fGainNormalisationFactor(1.)
   ,fPIDmethod(kLQ1D)
 {


@@ -58,8 +59,7 @@ AliTRDPIDResponse::AliTRDPIDResponse():
 //____________________________________________________________
 AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
   TObject(ref)
 //____________________________________________________________
 AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
   TObject(ref)

-  ,fkPIDReference(ref.fkPIDReference)
-  ,fkPIDParams(ref.fkPIDParams)
+  ,fkPIDResponseObject(NULL)

   ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
   ,fPIDmethod(ref.fPIDmethod)
 {
   ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
   ,fPIDmethod(ref.fPIDmethod)
 {


@@ -78,8 +78,7 @@ AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
   // Make copy
   TObject::operator=(ref);
   fGainNormalisationFactor = ref.fGainNormalisationFactor;
   // Make copy
   TObject::operator=(ref);
   fGainNormalisationFactor = ref.fGainNormalisationFactor;

-  fkPIDReference = ref.fkPIDReference;
-  fkPIDParams = ref.fkPIDParams;
+  fkPIDResponseObject = ref.fkPIDResponseObject;

   fPIDmethod = ref.fPIDmethod;
   
   return *this;
   fPIDmethod = ref.fPIDmethod;
   
   return *this;


@@ -90,11 +89,11 @@ AliTRDPIDResponse::~AliTRDPIDResponse(){
   //
   // Destructor
   //
   //
   // 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
   //
   //
   // Load References into the toolkit
   //


@@ -110,11 +109,11 @@ Bool_t AliTRDPIDResponse::Load(const Char_t * filename, const Char_t *refName){
   }
   
   gROOT->cd();
   }
   
   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);
   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;
 }
 
   return kTRUE;
 }
 


@@ -139,73 +138,110 @@ Bool_t AliTRDPIDResponse::GetResponse(Int_t n, const Double_t * const dedx, cons
   //   true if calculation success
   // 
 
   //   true if calculation success
   // 
 

-  if(!fkPIDReference){
-    AliWarning("Missing reference data. PID calculation not possible.");
-    return kFALSE;
-  }
+    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])) continue;
+    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;

 
 

-    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];
+       s=0.;
+        Bool_t filled=kTRUE;
+       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]));
+           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){
     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) 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));
   //
   // 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.;
   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(fPIDmethod){
+  case kNN: // NN
+      break;
+  case kLQ2D: // 2D LQ
+      {
+         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");
+         }
+      }
+      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("Probability %f", probLayer));
+      }
+      break;
+  default:
+      break;

   }
   }

-  AliDebug(1, Form("Probability %f", probLayer));

   return probLayer;
 }
 
   return probLayer;
 }
 


@@ -214,10 +250,10 @@ void AliTRDPIDResponse::SetOwner(){
   //
   // Make Deep Copy of the Reference Histograms
   //
   //
   // 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);

 }
 
 //____________________________________________________________
 }
 
 //____________________________________________________________


@@ -228,71 +264,54 @@ Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Doub
        //
   switch(fPIDmethod){
   case kNN: // NN 
        //
   switch(fPIDmethod){
   case kNN: // NN 

-    break;
-  case kLQ2D: // 2D LQ 
-    break;
-  case kLQ1D: // 1D LQ 
-    out[0]= 0.;
-    for(Int_t islice = 0; islice < nSlice; islice++) out[0] += in[islice] * fGainNormalisationFactor;
-    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];
+      }
+      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;
+

   default:
   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 {
   }
   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){
+    //
+    // 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]);
     AliError("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,fPIDmethod)){

     AliError("No Params found for the given configuration");
     return kTRUE;
   }
     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;
 }
 
   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];
-    }
-  }  
-
-  return parameters;
-}