Updates in order to enable the '2D' PID for the TRD developed by Daniel Lohner.
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliTRDPIDResponse.cxx
index edff1a1..fe69bca 100644 (file)
 
 #include "AliLog.h"
 
-#include "AliTRDPIDParams.h"
-#include "AliTRDPIDReference.h"
+#include "AliTRDPIDResponseObject.h"
+//#include "AliTRDPIDParams.h"
+//#include "AliTRDPIDReference.h"
 #include "AliTRDPIDResponse.h"
+#include "AliTRDTKDInterpolator.h"
 
 ClassImp(AliTRDPIDResponse)
 
 //____________________________________________________________
 AliTRDPIDResponse::AliTRDPIDResponse():
   TObject()
-  ,fkPIDReference(NULL)
-  ,fkPIDParams(NULL)
+  ,fkPIDResponseObject(NULL)
   ,fGainNormalisationFactor(1.)
   ,fPIDmethod(kLQ1D)
 {
@@ -58,8 +59,7 @@ AliTRDPIDResponse::AliTRDPIDResponse():
 //____________________________________________________________
 AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
   TObject(ref)
-  ,fkPIDReference(ref.fkPIDReference)
-  ,fkPIDParams(ref.fkPIDParams)
+  ,fkPIDResponseObject(NULL)
   ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
   ,fPIDmethod(ref.fPIDmethod)
 {
@@ -78,8 +78,7 @@ AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
   // Make copy
   TObject::operator=(ref);
   fGainNormalisationFactor = ref.fGainNormalisationFactor;
-  fkPIDReference = ref.fkPIDReference;
-  fkPIDParams = ref.fkPIDParams;
+  fkPIDResponseObject = ref.fkPIDResponseObject;
   fPIDmethod = ref.fPIDmethod;
   
   return *this;
@@ -90,11 +89,11 @@ AliTRDPIDResponse::~AliTRDPIDResponse(){
   //
   // 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
   //
@@ -110,11 +109,11 @@ 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;
 }
 
@@ -139,73 +138,110 @@ Bool_t AliTRDPIDResponse::GetResponse(Int_t n, const Double_t * const dedx, cons
   //   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){
-      AliDebug(2, Form("Null all species prob layer[%d].", il));
-      continue;
+       AliDebug(2, "Null total prob.");
+       return kFALSE;
     }
-    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 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));
-  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(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;
 }
 
@@ -214,10 +250,10 @@ 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);
 }
 
 //____________________________________________________________
@@ -228,42 +264,49 @@ Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Doub
        //
   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++) 
-      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];
+      }
+      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:
-    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){
+    //
+    // 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(!fkPIDParams->GetThresholdParameters(nTracklets, level, params)){
+  if(!fkPIDResponseObject->GetThresholdParameters(nTracklets, level, params,fPIDmethod)){
     AliError("No Params found for the given configuration");
     return kTRUE;
   }