Updates in order to enable the '2D' PID for the TRD developed by Daniel Lohner.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 27 Jun 2012 15:21:03 +0000 (15:21 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 27 Jun 2012 15:21:03 +0000 (15:21 +0000)
15 files changed:
ANALYSIS/TenderSupplies/AliTRDTenderSupply.cxx
STEER/CMakelibSTEERBase.pkg
STEER/STEERBase/AliPIDResponse.cxx
STEER/STEERBase/AliPIDResponse.h
STEER/STEERBase/AliTRDPIDParams.cxx
STEER/STEERBase/AliTRDPIDParams.h
STEER/STEERBase/AliTRDPIDReference.cxx
STEER/STEERBase/AliTRDPIDResponse.cxx
STEER/STEERBase/AliTRDPIDResponse.h
STEER/STEERBase/AliTRDPIDResponseObject.cxx [new file with mode: 0644]
STEER/STEERBase/AliTRDPIDResponseObject.h [new file with mode: 0644]
STEER/STEERBase/AliTRDTKDInterpolator.cxx [new file with mode: 0644]
STEER/STEERBase/AliTRDTKDInterpolator.h [new file with mode: 0644]
STEER/STEERBaseLinkDef.h
TRD/AliTRDcalibDB.cxx

index d4b3763..0f798f7 100644 (file)
@@ -47,7 +47,7 @@
 #include <AliESDInputHandler.h>
 #include <AliAnalysisManager.h>
 #include <AliTrackerBase.h>
-#include <AliTRDPIDReference.h>
+#include <AliTRDPIDResponseObject.h>
 #include <AliTRDPIDResponse.h>
 #include "AliTRDCalChamberStatus.h"
 #include <AliTender.h>
@@ -276,15 +276,15 @@ void AliTRDTenderSupply::LoadReferences(){
     // Get new references
     TIter refs(arr);
     TObject *o = NULL;
-    AliTRDPIDReference *ref = NULL;
+    AliTRDPIDResponseObject *ref = NULL;
     while((o = refs())){
-      if(!TString(o->IsA()->GetName()).CompareTo("AliTRDPIDReference")){
-        ref = dynamic_cast<AliTRDPIDReference *>(o);
+      if(!TString(o->IsA()->GetName()).CompareTo("AliTRDPIDResponseObject")){
+        ref = dynamic_cast<AliTRDPIDResponseObject *>(o);
         break;
       }
     }
     if(ref){
-      fESDpid->GetTRDResponse().Load(ref);
+      fESDpid->GetTRDResponse().SetPIDResponseObject(ref);
       fHasReferences = kTRUE;
       AliDebug(1, "Reference distributions loaded into the PID Response");
     } else {
index d7916e0..ef4d53d 100644 (file)
@@ -62,11 +62,13 @@ set ( SRCS
     STEERBase/AliITSPidParams.cxx 
     STEERBase/AliTRDPIDReference.cxx 
     STEERBase/AliTRDPIDParams.cxx 
-    STEERBase/AliTOFPIDParams.cxx 
+    STEERBase/AliTRDPIDResponseObject.cxx
+    STEERBase/AliTRDTKDInterpolator.cxx
     STEERBase/AliPIDResponse.cxx 
     STEERBase/AliITSPIDResponse.cxx 
     STEERBase/AliTPCPIDResponse.cxx 
     STEERBase/AliTOFPIDResponse.cxx 
+    STEERBase/AliTOFPIDParams.cxx 
     STEERBase/AliTRDPIDResponse.cxx 
     STEERBase/AliEMCALPIDResponse.cxx 
     STEERBase/AliPIDCombined.cxx
index 511475d..57ed3ff 100644 (file)
 #include <TF1.h>
 #include <TSpline.h>
 #include <TFile.h>
+#include <TArrayF.h>
 
 #include <AliVEvent.h>
 #include <AliVTrack.h>
 #include <AliLog.h>
 #include <AliPID.h>
 #include <AliOADBContainer.h>
-#include <AliTRDPIDParams.h>
-#include <AliTRDPIDReference.h>
+#include <AliTRDPIDResponseObject.h>
 #include <AliTOFPIDParams.h>
 
 #include "AliPIDResponse.h"
@@ -67,8 +67,7 @@ fRun(0),
 fOldRun(0),
 fArrPidResponseMaster(0x0),
 fResolutionCorrection(0x0),
-fTRDPIDParams(0x0),
-fTRDPIDReference(0x0),
+fTRDPIDResponseObject(0x0),
 fTOFtail(1.1),
 fTOFPIDParams(0x0),
 fEMCALPIDParams(0x0),
@@ -91,9 +90,8 @@ AliPIDResponse::~AliPIDResponse()
   //
   // dtor
   //
-  delete fArrPidResponseMaster;
-  delete fTRDPIDParams;
-  delete fTRDPIDReference;
+    delete fArrPidResponseMaster;
+    delete fTRDPIDResponseObject;
   if (fTOFPIDParams) delete fTOFPIDParams;
 }
 
@@ -120,8 +118,7 @@ fRun(0),
 fOldRun(0),
 fArrPidResponseMaster(0x0),
 fResolutionCorrection(0x0),
-fTRDPIDParams(0x0),
-fTRDPIDReference(0x0),
+fTRDPIDResponseObject(0x0),
 fTOFtail(1.1),
 fTOFPIDParams(0x0),
 fEMCALPIDParams(0x0),
@@ -163,8 +160,7 @@ AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
     fOldRun=0;
     fArrPidResponseMaster=0x0;
     fResolutionCorrection=0x0;
-    fTRDPIDParams=0x0;
-    fTRDPIDReference=0x0;
+    fTRDPIDResponseObject=0x0;
     fEMCALPIDParams=0x0;
     memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
     fTOFtail=1.1;
@@ -846,20 +842,20 @@ void AliPIDResponse::SetTRDPidResponseMaster()
   //
   // Load the TRD pid params and references from the OADB
   //
-  if(fTRDPIDParams) return;
+  if(fTRDPIDResponseObject) return;
   AliOADBContainer contParams("contParams"); 
 
-  Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()), "AliTRDPIDParams");
-  if(statusPars){
-    AliError("Failed initializing PID Params from OADB");
+  Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject");
+  if(statusResponse){
+    AliError("Failed initializing PID Response Object from OADB");
   } else {
-    AliInfo(Form("Loading TRD Params from %s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()));
-    fTRDPIDParams = dynamic_cast<AliTRDPIDParams *>(contParams.GetObject(fRun));
-    if(!fTRDPIDParams){
-      AliError(Form("TRD Params not found in run %d", fRun));
+    AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()));
+    fTRDPIDResponseObject = dynamic_cast<AliTRDPIDResponseObject *>(contParams.GetObject(fRun));
+    if(!fTRDPIDResponseObject){
+      AliError(Form("TRD Response not found in run %d", fRun));
     }
   }
-
+  /*
   AliOADBContainer contRefs("contRefs");
   Int_t statusRefs = contRefs.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()), "AliTRDPIDReference");
   if(statusRefs){
@@ -870,7 +866,8 @@ void AliPIDResponse::SetTRDPidResponseMaster()
     if(!fTRDPIDReference){
       AliError(Form("TRD References not found in OADB Container for run %d", fRun));
     }
-  }
+    }
+    */
 }
 
 //______________________________________________________________________________
@@ -878,12 +875,11 @@ void AliPIDResponse::InitializeTRDResponse(){
   //
   // Set PID Params and references to the TRD PID response
   // 
-  fTRDResponse.SetPIDParams(fTRDPIDParams);
-  fTRDResponse.Load(fTRDPIDReference);
-  if(fLHCperiod == "LHC10b" || fLHCperiod == "LHC10c" || fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
-    fTRDslicesForPID[0] = 0;
-    fTRDslicesForPID[1] = 7;
-  }
+    fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject);
+    if(fLHCperiod == "LHC10b" || fLHCperiod == "LHC10c" || fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
+       fTRDslicesForPID[0] = 0;
+       fTRDslicesForPID[1] = 7;
+    }
 }
 
 //______________________________________________________________________________
index 6f11e9f..d715730 100644 (file)
@@ -25,6 +25,7 @@
 
 class AliVEvent;
 class TF1;
+class AliTRDPIDResponseObject; 
 
 class AliPIDResponse : public TNamed {
 public:
@@ -130,8 +131,7 @@ private:
   TObjArray *fArrPidResponseMaster;    //!  TPC pid splines
   TF1       *fResolutionCorrection;    //! TPC resolution correction
 
-  AliTRDPIDParams *fTRDPIDParams;       //! TRD PID Params
-  AliTRDPIDReference *fTRDPIDReference; //! TRD PID References
+  AliTRDPIDResponseObject *fTRDPIDResponseObject; //! TRD PID Response Object
   UInt_t fTRDslicesForPID[2];           //! TRD PID slices
 
   Float_t fTOFtail;                    //! TOF tail effect used in TOF probability
index 5f8e60c..950ac42 100644 (file)
@@ -49,6 +49,30 @@ AliTRDPIDParams::AliTRDPIDParams(const char *name) :
   fEntries = new TSortedList;
 }
 
+//____________________________________________________________
+AliTRDPIDParams::AliTRDPIDParams(const AliTRDPIDParams &ref):
+TNamed(ref),
+fEntries(NULL)
+{
+    //
+    // Copy constructor
+    //
+
+    fEntries=(TSortedList*)ref.fEntries->Clone();
+}
+
+//____________________________________________________________
+AliTRDPIDParams::AliTRDPIDParams(const AliTRDPIDParams &ref):
+TNamed(ref),
+fEntries(NULL)
+{
+    //
+    // Copy constructor
+    //
+
+    fEntries=(TSortedList*)ref.fEntries->Clone();
+}
+
 //____________________________________________________________
 AliTRDPIDParams::~AliTRDPIDParams(){
   //
index bbadcae..0d46a9f 100644 (file)
@@ -28,6 +28,7 @@ class AliTRDPIDParams : public TNamed{
   public:
     AliTRDPIDParams();
     AliTRDPIDParams(const char *name);
+    AliTRDPIDParams(const AliTRDPIDParams &);
     virtual ~AliTRDPIDParams();
     virtual void Print(Option_t *) const;
 
@@ -36,35 +37,34 @@ class AliTRDPIDParams : public TNamed{
 
   private:
     class AliTRDPIDThresholds : public TObject{
-      public:
-        AliTRDPIDThresholds();
-        AliTRDPIDThresholds(Int_t nTracklets, Double_t effMin, Double_t effMax, Double_t *params = NULL);
-        AliTRDPIDThresholds(Int_t nTracklets, Double_t eff, Double_t *params = NULL);
-        AliTRDPIDThresholds(const AliTRDPIDThresholds &);
-        AliTRDPIDThresholds &operator=(const AliTRDPIDThresholds &);
-        virtual ~AliTRDPIDThresholds() {}
-        
+    public:
+       AliTRDPIDThresholds();
+       AliTRDPIDThresholds(Int_t nTracklets, Double_t effMin, Double_t effMax, Double_t *params = NULL);
+       AliTRDPIDThresholds(Int_t nTracklets, Double_t eff, Double_t *params = NULL);
+       AliTRDPIDThresholds(const AliTRDPIDThresholds &);
+       AliTRDPIDThresholds &operator=(const AliTRDPIDThresholds &);
+       virtual ~AliTRDPIDThresholds() {}
+
         Int_t GetNTracklets() const { return fNTracklets; }
-        Double_t GetElectronEfficiency(Int_t step = 0) const { if(step == 0) return fEfficiency[0]; else return fEfficiency[1]; }
-        const Double_t *GetThresholdParams() const { return fParams; }
+       Double_t GetElectronEfficiency(Int_t step = 0) const { if(step == 0) return fEfficiency[0]; else return fEfficiency[1]; }
+       const Double_t *GetThresholdParams() const { return fParams; }
 
-        virtual Bool_t IsSortable() const { return kTRUE; }
-        virtual Int_t Compare(const TObject *ref) const;
+       virtual Bool_t IsSortable() const { return kTRUE; }
+       virtual Int_t Compare(const TObject *ref) const;
         virtual Bool_t IsEqual(const TObject *ref) const;
 
-      private:
-        Int_t fNTracklets;          //
-        Double_t fEfficiency[2];    //
-        Double_t fParams[4];        //
+    private:
+       Int_t fNTracklets;          //
+       Double_t fEfficiency[2];    //
+       Double_t fParams[4];        //
 
-        ClassDef(AliTRDPIDThresholds, 1);
+       ClassDef(AliTRDPIDThresholds, 1);
     };
 
+    AliTRDPIDParams &operator=(const AliTRDPIDParams &);
+  
     static const Double_t kVerySmall;
 
-    AliTRDPIDParams(const AliTRDPIDParams &);
-    AliTRDPIDParams &operator=(const AliTRDPIDParams &);
-    
     TSortedList *fEntries; //
 
     ClassDef(AliTRDPIDParams, 1);
index e3e6379..c50707d 100644 (file)
@@ -58,14 +58,26 @@ AliTRDPIDReference::AliTRDPIDReference(const Char_t *name):
 //____________________________________________________________
 AliTRDPIDReference::AliTRDPIDReference(const AliTRDPIDReference &ref):
                  TNamed(ref),
-                 fRefContainer(ref.fRefContainer),
+                 fRefContainer(NULL),
                  fMomentumBins(ref.fMomentumBins)
 {
-       //
-       // Copy constructor
-       // Only copies poiters, object is not the owner of the references
-       //
-       SetBit(kIsOwner, kFALSE);
+    //
+    // Copy constructor
+    //
+    fRefContainer = new TObjArray(fMomentumBins.GetSize() * AliPID::kSPECIES);
+    fRefContainer->SetOwner();
+
+    for(Int_t ip = 0; ip < GetNumberOfMomentumBins(); ip++){
+       for(Int_t is = 0; is < 5; is++){
+           Int_t ent=is * fMomentumBins.GetSize() + ip;
+           TObject *obj=ref.fRefContainer->At(ent);
+           if(obj){
+               fRefContainer->AddAt(obj->Clone(),ent);
+           }
+       }
+    }
+
+    SetBit(kIsOwner, kTRUE);
 }
 
 //____________________________________________________________
@@ -124,7 +136,7 @@ void AliTRDPIDReference::AddReference(TObject *ref, AliPID::EParticleType spec,
                return;
        }
        AliDebug(1, Form("Adding object with address %p to position %d", ref, spec * fMomentumBins.GetSize() + pbin));
-       fRefContainer->AddAt(ref, spec * fMomentumBins.GetSize() + pbin);
+       fRefContainer->AddAt(ref->Clone(), spec * fMomentumBins.GetSize() + pbin);
 }
 
 //____________________________________________________________
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;
   }
index 53d90de..f1c079b 100644 (file)
 
 class TObjArray;
 class AliVTrack;
-class AliTRDPIDParams;
-class AliTRDPIDReference;
+class AliTRDPIDResponseObject;
+
 class AliTRDPIDResponse : public TObject {
   public:
     enum ETRDPIDResponseStatus {
       kIsOwner = BIT(14)
     };
     enum ETRDPIDResponseDef {
-      kNlayer = 6
-      ,kNPBins = 6
+       kNlayer = 6
+       ,kNPBins = 6
     };
     enum ETRDPIDMethod {
       kNN   = 0,
       kLQ2D = 1,
       kLQ1D = 2
     };
+    enum ETRDPIDNMethod {
+       kNMethod=3
+    };
     enum ETRDNslices {
       kNslicesLQ1D = 1,
       kNslicesLQ2D = 2,
@@ -65,19 +68,18 @@ class AliTRDPIDResponse : public TObject {
     void      SetOwner();
     void      SetPIDmethod(ETRDPIDMethod m) {fPIDmethod=m;}
     void      SetGainNormalisationFactor(Double_t gainFactor) { fGainNormalisationFactor = gainFactor; }
-    void      SetPIDParams(const AliTRDPIDParams * params) { fkPIDParams = params; }
 
-    Bool_t    Load(const Char_t *filename = NULL, const Char_t *refName = "RefTRDLQ1D");
-    Bool_t    Load(const AliTRDPIDReference *ref) { fkPIDReference = ref; return kTRUE; }
+    Bool_t    SetPIDResponseObject(const AliTRDPIDResponseObject * obj) { fkPIDResponseObject = obj; return kTRUE;}
+
+    Bool_t    Load(const Char_t *filename = NULL);
   
     Bool_t    IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level) const;
   
   private:
     Bool_t    CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out) const;
-    Double_t  GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t dEdx) const;
+    Double_t  GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t *dEdx) const;
     
-    const AliTRDPIDReference *fkPIDReference;   // PID References
-    const AliTRDPIDParams *fkPIDParams;         // PID Params
+    const AliTRDPIDResponseObject *fkPIDResponseObject;   // PID References and Params
     Double_t  fGainNormalisationFactor;         // Gain normalisation factor
     ETRDPIDMethod   fPIDmethod;                 // PID method selector
       
@@ -96,3 +98,4 @@ AliTRDPIDResponse::ETRDNslices AliTRDPIDResponse::GetNumberOfSlices() const {
   return slices;
 }
 #endif
+
diff --git a/STEER/STEERBase/AliTRDPIDResponseObject.cxx b/STEER/STEERBase/AliTRDPIDResponseObject.cxx
new file mode 100644 (file)
index 0000000..566fa90
--- /dev/null
@@ -0,0 +1,219 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+//
+// Container class for the reference distributions for TRD PID
+// The class contains the reference distributions and the momentum steps
+// the references are taken at. Mapping is done inside. To derive references,
+// the functions GetUpperReference and GetLowerReference return the next
+// reference distribution object and the momentum step above respectively below
+// the tracklet momentum.
+//
+// Authors:
+//    Markus Fasel <M.Fasel@gsi.de>
+//    Daniel Lohner <Daniel.Lohner@cern.ch>
+
+#include "AliLog.h"
+
+#include "AliTRDPIDResponseObject.h"
+
+#ifndef AliTRDPIDREFERENCE_H
+#include "AliTRDPIDReference.h"
+#endif
+
+#ifndef AliTRDPIDPARAMS_H
+#include "AliTRDPIDParams.h"
+#endif
+
+
+ClassImp(AliTRDPIDResponseObject)
+
+//____________________________________________________________
+AliTRDPIDResponseObject::AliTRDPIDResponseObject():
+    TNamed(),
+    fNSlicesQ0(4)
+{
+    //
+    // Dummy constructor
+    //
+    SetBit(kIsOwner, kTRUE);
+
+    for(Int_t method=0;method<AliTRDPIDResponse::kNMethod;method++){
+       fkPIDParams[method]=NULL;
+       fkPIDReference[method]=NULL;
+    }
+}
+
+//____________________________________________________________
+AliTRDPIDResponseObject::AliTRDPIDResponseObject(const Char_t *name):
+TNamed(name, "TRD PID Response Object"),
+fNSlicesQ0(4)
+{
+       //
+       // Default constructor
+       //
+       SetBit(kIsOwner, kTRUE);
+
+       for(Int_t method=0;method<AliTRDPIDResponse::kNMethod;method++){
+           fkPIDParams[method]=NULL;
+           fkPIDReference[method]=NULL;
+       }
+}
+
+//____________________________________________________________
+AliTRDPIDResponseObject::AliTRDPIDResponseObject(const AliTRDPIDResponseObject &ref):
+TNamed(ref),
+fNSlicesQ0(ref.fNSlicesQ0)
+{
+    //
+    // Copy constructor
+    // Only copies pointers, object is not the owner of the references
+    //
+    SetBit(kIsOwner, kFALSE);
+
+    for(Int_t method=0;method<AliTRDPIDResponse::kNMethod;method++){
+       fkPIDParams[method]=ref.fkPIDParams[method];       // new Object is not owner, copy only pointer
+       fkPIDReference[method]=ref.fkPIDReference[method];    // new Object is not owner, copy only pointer
+    }
+}
+//____________________________________________________________
+AliTRDPIDResponseObject &AliTRDPIDResponseObject::operator=(const AliTRDPIDResponseObject &ref){
+       //
+       // Assginment operator
+       // Only copies poiters, object is not the owner of the references
+       //
+       if(this != &ref){
+           TNamed::operator=(ref);
+           fNSlicesQ0=ref.fNSlicesQ0;
+           for(Int_t method=0;method<AliTRDPIDResponse::kNMethod;method++){
+               if(TestBit(kIsOwner) && fkPIDParams[method])delete fkPIDParams[method];
+                if(TestBit(kIsOwner) && fkPIDReference[method])delete fkPIDReference[method];
+
+               fkPIDParams[method]=ref.fkPIDParams[method];       // new Object is not owner, copy only pointer
+               fkPIDReference[method]=ref.fkPIDReference[method];    // new Object is not owner, copy only pointer
+           }
+           SetBit(kIsOwner, kFALSE);
+       }
+       return *this;
+}
+
+//____________________________________________________________
+AliTRDPIDResponseObject::~AliTRDPIDResponseObject(){
+       //
+       // Destructor
+       // references are deleted if the object is the owner
+       //
+    if(fkPIDParams && TestBit(kIsOwner)) delete fkPIDParams;
+    if(fkPIDReference && TestBit(kIsOwner)) delete fkPIDReference;
+
+}
+
+//____________________________________________________________
+void AliTRDPIDResponseObject::SetPIDParams(AliTRDPIDParams *params,AliTRDPIDResponse::ETRDPIDMethod method){
+
+    if(Int_t(method)>=Int_t(AliTRDPIDResponse::kNMethod)||Int_t(method)<0){
+       AliError("Method does not exist");
+       return;
+    }
+    if(fkPIDParams[method]){
+       delete fkPIDParams[method];
+        fkPIDParams[method]=NULL;
+    }
+
+    fkPIDParams[method]=new AliTRDPIDParams(*params);
+}
+
+//____________________________________________________________
+void AliTRDPIDResponseObject::SetPIDReference(AliTRDPIDReference *reference,AliTRDPIDResponse::ETRDPIDMethod method){
+
+    if(Int_t(method)>=Int_t(AliTRDPIDResponse::kNMethod)||Int_t(method)<0){
+        AliError("Method does not exist");
+       return;
+    }
+    if(fkPIDReference[method]){
+       delete fkPIDReference[method];
+       fkPIDReference[method]=NULL;
+    }
+    fkPIDReference[method]=new AliTRDPIDReference(*reference);
+}
+
+//____________________________________________________________
+TObject *AliTRDPIDResponseObject::GetUpperReference(AliPID::EParticleType spec, Float_t p, Float_t &pUpper,AliTRDPIDResponse::ETRDPIDMethod method) const{
+
+    if(Int_t(method)>=Int_t(AliTRDPIDResponse::kNMethod)||Int_t(method)<0){
+       AliError("Method does not exist");
+       return NULL;
+    }
+   
+    if(fkPIDReference[method]){
+       return fkPIDReference[method]->GetUpperReference(spec,p,pUpper);
+    }
+    return NULL;
+}
+
+
+//____________________________________________________________
+TObject *AliTRDPIDResponseObject::GetLowerReference(AliPID::EParticleType spec, Float_t p, Float_t &pLower,AliTRDPIDResponse::ETRDPIDMethod method) const{
+
+    if(Int_t(method)>=Int_t(AliTRDPIDResponse::kNMethod)||Int_t(method)<0){
+       AliError("Method does not exist");
+       return NULL;
+    }
+
+     if(fkPIDReference[method]){
+        return fkPIDReference[method]->GetLowerReference(spec,p,pLower);
+     }
+    return NULL;
+}
+
+//____________________________________________________________
+Bool_t AliTRDPIDResponseObject::GetThresholdParameters(Int_t ntracklets, Double_t efficiency, Double_t *params,AliTRDPIDResponse::ETRDPIDMethod method) const{
+
+    if(Int_t(method)>=Int_t(AliTRDPIDResponse::kNMethod)||Int_t(method)<0){
+       AliError("Method does not exist");
+       return kFALSE;
+    }
+
+    if(fkPIDParams[method]){
+       return fkPIDParams[method]->GetThresholdParameters(ntracklets,efficiency,params);
+    }
+    return kFALSE;
+}
+
+//____________________________________________________________
+Int_t AliTRDPIDResponseObject::GetNumberOfMomentumBins(AliTRDPIDResponse::ETRDPIDMethod method) const{
+
+    if(Int_t(method)>=Int_t(AliTRDPIDResponse::kNMethod)||Int_t(method)<0){
+       AliError("Method does not exist");
+       return 0;
+    }
+
+    if(fkPIDReference[method]){
+       return fkPIDReference[method]->GetNumberOfMomentumBins();
+    }
+    return 0;
+}
+
+//____________________________________________________________
+void AliTRDPIDResponseObject::Print(const Option_t* opt) const{
+       //
+       // Print content of the PID object
+       //
+    printf("Content of AliTRDPIDResponseObject \n\n");
+   
+    for(Int_t method=0;method<AliTRDPIDResponse::kNMethod;method++){
+       if(fkPIDReference[method])fkPIDReference[method]->Print(opt);
+       if(fkPIDParams[method])printf("+ Threshold Parameters \n");
+    }
+}
diff --git a/STEER/STEERBase/AliTRDPIDResponseObject.h b/STEER/STEERBase/AliTRDPIDResponseObject.h
new file mode 100644 (file)
index 0000000..701dc5a
--- /dev/null
@@ -0,0 +1,72 @@
+/**************************************************************************
+* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+*                                                                        *
+* Author: The ALICE Off-line Project.                                    *
+* Contributors are mentioned in the code where appropriate.              *
+*                                                                        *
+* Permission to use, copy, modify and distribute this software and its   *
+* documentation strictly for non-commercial purposes is hereby granted   *
+* without fee, provided that the above copyright notice appears in all   *
+* copies and that both the copyright notice and this permission notice   *
+* appear in the supporting documentation. The authors make no claims     *
+* about the suitability of this software for any purpose. It is          *
+* provided "as is" without express or implied warranty.                  *
+**************************************************************************/
+//
+// Container for TRD PID Response Objects stored in the OADB
+//
+#ifndef ALITRDPIDRESPONSEOBJECT_H
+#define ALITRDPIDRESPONSEOBJECT_H
+
+#ifndef ROOT_TNamed
+#include <TNamed.h>
+#endif
+
+#ifndef AliTRDPIDRESPONSE_H
+#include "AliTRDPIDResponse.h"
+#endif
+
+
+class AliTRDPIDParams;
+class AliTRDPIDReference;
+class AliTRDPIDResponse;
+
+class AliTRDPIDResponseObject : public TNamed{
+public:
+    enum ETRDPIDResponseObjectStatus {
+       kIsOwner = BIT(14)
+    };
+
+    AliTRDPIDResponseObject();
+    AliTRDPIDResponseObject(const char *name);
+    AliTRDPIDResponseObject(const AliTRDPIDResponseObject &ref);
+    AliTRDPIDResponseObject &operator=(const AliTRDPIDResponseObject &ref);
+    virtual ~AliTRDPIDResponseObject();
+
+    virtual void Print(Option_t *opt) const;
+
+    void SetPIDParams(AliTRDPIDParams *params,AliTRDPIDResponse::ETRDPIDMethod method=AliTRDPIDResponse::kLQ1D);
+    void SetPIDReference(AliTRDPIDReference *params,AliTRDPIDResponse::ETRDPIDMethod method=AliTRDPIDResponse::kLQ1D);
+
+    // Derive reference
+    TObject *GetLowerReference(AliPID::EParticleType spec, Float_t p, Float_t &pLower,AliTRDPIDResponse::ETRDPIDMethod method=AliTRDPIDResponse::kLQ1D) const;
+    TObject *GetUpperReference(AliPID::EParticleType spec, Float_t p, Float_t &pUpper,AliTRDPIDResponse::ETRDPIDMethod method=AliTRDPIDResponse::kLQ1D) const;
+
+    Int_t GetNumberOfMomentumBins(AliTRDPIDResponse::ETRDPIDMethod method=AliTRDPIDResponse::kLQ1D) const;
+
+    // Derive threshold params
+    Bool_t GetThresholdParameters(Int_t ntracklets, Double_t efficiency, Double_t *params,AliTRDPIDResponse::ETRDPIDMethod method=AliTRDPIDResponse::kLQ1D) const;
+
+    // Number of SlicesQ0
+    Int_t GetNSlicesQ0() const{return fNSlicesQ0;}
+    void SetNSlicesQ0(Int_t nsl){fNSlicesQ0=nsl;}
+
+private:
+
+      const AliTRDPIDParams *fkPIDParams[AliTRDPIDResponse::kNMethod]; // Contains Thresholds
+      const AliTRDPIDReference *fkPIDReference[AliTRDPIDResponse::kNMethod]; // Contains References
+      Int_t fNSlicesQ0; // Number of Slices for Q0
+
+    ClassDef(AliTRDPIDResponseObject, 1);
+};
+#endif
diff --git a/STEER/STEERBase/AliTRDTKDInterpolator.cxx b/STEER/STEERBase/AliTRDTKDInterpolator.cxx
new file mode 100644 (file)
index 0000000..d2842ed
--- /dev/null
@@ -0,0 +1,588 @@
+#include "AliTRDTKDInterpolator.h"
+
+#include "TClonesArray.h"
+#include "TLinearFitter.h"
+#include "TMath.h"
+#include "TRandom.h"
+
+#include "iostream"
+using namespace std;
+
+ClassImp(AliTRDTKDInterpolator)
+
+//_________________________________________________________________
+AliTRDTKDInterpolator::AliTRDTKDInterpolator() :
+TKDTreeIF(),
+fNDataNodes(0),
+fNodes(NULL),
+fLambda(0),
+fNPointsI(0),
+fUseHelperNodes(kFALSE),
+fUseWeights(kFALSE),
+fPDFMode(kInterpolation)
+{
+  // default constructor
+}
+
+//_________________________________________________________________
+AliTRDTKDInterpolator::AliTRDTKDInterpolator(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data) :
+TKDTreeIF(npoints, ndim, bsize, data),
+fNDataNodes(0),
+fNodes(NULL),
+fLambda(1 + ndim + (ndim*(ndim+1)>>1)),
+fNPointsI(100),
+fUseHelperNodes(kFALSE),
+fUseWeights(kFALSE),
+fPDFMode(kInterpolation)
+{
+}
+
+//_________________________________________________________________
+AliTRDTKDInterpolator::~AliTRDTKDInterpolator()
+{
+    if(fNodes){
+       fNodes->Delete();
+       delete fNodes;
+       fNodes=NULL;
+    }
+}
+//_________________________________________________________________
+AliTRDTKDInterpolator::AliTRDTKDInterpolator(const AliTRDTKDInterpolator &ref):
+TKDTreeIF(),
+fNDataNodes(ref.fNDataNodes),
+fNodes(ref.fNodes),
+fLambda(ref.fLambda),
+fNPointsI(ref.fNPointsI),
+fUseHelperNodes(ref.fUseHelperNodes),
+fUseWeights(ref.fUseWeights),
+fPDFMode(ref.fPDFMode)
+{
+    // Copy constructor
+}
+
+//____________________________________________________________
+
+AliTRDTKDInterpolator &AliTRDTKDInterpolator::operator=(const AliTRDTKDInterpolator &ref){
+    //
+    // Assignment operator
+    //
+    if(this == &ref) return *this;
+
+    // Make copy
+    TObject::operator=(ref);
+
+    return *this;
+}
+
+//_________________________________________________________________
+Bool_t AliTRDTKDInterpolator::Build()
+{
+    TKDTreeIF::Build();
+    if(!fBoundaries) MakeBoundaries();
+
+    // allocate interpolation nodes
+    fNDataNodes = fNPoints/fBucketSize + ((fNPoints%fBucketSize)?1:0);/*TKDTreeIF::GetNTNodes();*/
+
+    if(fNodes){
+       Warning("AliTRDTKDInterpolator::Build()", "Data already allocated.");
+       fNodes->Delete();
+    } else {
+       fNodes = new TClonesArray("AliTRDTKDInterpolator::AliTRDTKDNodeInfo", fNDataNodes);
+       fNodes->SetOwner();
+    }
+    for(int in=0; in<fNDataNodes; in++) new ((*fNodes)[in]) AliTRDTKDNodeInfo(fNDim);
+
+    // Set Interpolator nodes
+
+    for(int inode=0, tnode = fNNodes; inode<fNDataNodes-1; inode++, tnode++){
+       AliTRDTKDNodeInfo *node =GetNodeInfo(inode);
+       memcpy(node->fBounds,GetBoundary(tnode),2*fNDim*sizeof(Float_t));
+       node->fVal[0] =  Float_t(fBucketSize)/fNPoints;
+       for(int idim=0; idim<fNDim; idim++) node->fVal[0] /= (node->fBounds[2*idim+1] - node->fBounds[2*idim]);
+       node->fVal[1] =  node->fVal[0]/TMath::Sqrt(float(fBucketSize));
+
+       Int_t *indexPoints = GetPointsIndexes(tnode);
+       // loop points in this terminal node
+       for(int idim=0; idim<fNDim; idim++){
+           node->fData[idim] = 0.;
+           for(int ip = 0; ip<fBucketSize; ip++) node->fData[idim] += fData[idim][indexPoints[ip]];
+           node->fData[idim] /= fBucketSize;
+       }
+    }
+
+    // Analyze last (incomplete) terminal node
+
+    Int_t counts = fNPoints%fBucketSize;
+    counts = counts ? counts : fBucketSize;
+    Int_t inode = fNDataNodes - 1, tnode = inode + fNNodes;
+    AliTRDTKDNodeInfo *ftnode = GetNodeInfo(inode);
+    ftnode->fVal[0] =  Float_t(counts)/fNPoints;
+    memcpy(ftnode->fBounds,GetBoundary(tnode),2*fNDim*sizeof(Float_t));
+    for(int idim=0; idim<fNDim; idim++){
+       Float_t dx = ftnode->fBounds[2*idim+1]-ftnode->fBounds[2*idim];
+       if(dx < 1.e-30){
+           Warning("AliTRDTKDInterpolator::Build()", "Terminal bucket index[%d] too narrow on the %d dimension.", inode, idim);
+           continue;
+       }
+       ftnode->fVal[0] /= (ftnode->fBounds[2*idim+1] - ftnode->fBounds[2*idim]);
+    }
+    ftnode->fVal[1] =  ftnode->fVal[0]/TMath::Sqrt(float(counts));
+
+    // loop points in this terminal node
+    Int_t *indexPoints = GetPointsIndexes(tnode);
+    for(int idim=0; idim<fNDim; idim++){
+       ftnode->fData[idim] = 0.;
+       for(int ip = 0; ip<counts; ip++) ftnode->fData[idim] += fData[idim][indexPoints[ip]];
+       ftnode->fData[idim] /= counts;
+    }
+
+    delete [] fBoundaries;
+    fBoundaries = NULL;
+    // Add Helper Nodes
+    if(fUseHelperNodes){BuildBoundaryNodes();}
+
+    if(fNPointsI>GetNTNodes()){fNPointsI=GetNTNodes();}
+
+    BuildInterpolation();
+
+    return kTRUE;
+}
+
+
+//_________________________________________________________________
+Bool_t AliTRDTKDInterpolator::Eval(const Double_t *point, Double_t &result, Double_t &error)
+{
+    Float_t pointF[fNDim]; // local Float_t conversion for "point"
+    for(int idim=0; idim<fNDim; idim++) pointF[idim] = (Float_t)point[idim];
+    Int_t nodeIndex = GetNodeIndex(pointF);
+    if(nodeIndex<0){
+       Error("AliTRDTKDInterpolator::Eval()", "Can not retrive node for data point.");
+       result = 0.;
+       error = 1.E10;
+       return kFALSE;
+    }
+    AliTRDTKDNodeInfo *node =GetNodeInfo(nodeIndex);
+
+    switch(fPDFMode){
+
+    case kInterpolation:
+       return node->CookPDF(point, result, error);
+    case kMinError:
+       node->CookPDF(point, result, error);
+       if(error<node->fVal[1]){
+           return kTRUE;
+       }
+       error=node->fVal[1];
+       result=node->fVal[0];
+       return kTRUE;
+    case kNodeVal:
+       error=node->fVal[1];
+       result=node->fVal[0];
+       return kTRUE;
+    default:
+       return kFALSE;
+    }
+}
+
+//__________________________________________________________________
+void AliTRDTKDInterpolator::Print(const Option_t */*opt*/) const
+{
+    for(Int_t i=GetNTNodes(); i--;){
+       printf("Node %d of %d: \n",i,GetNTNodes());
+       GetNodeInfo(i)->Print();
+    }
+
+}
+
+//__________________________________________________________________
+Int_t AliTRDTKDInterpolator::GetNodeIndex(const Float_t *p)
+{
+    Int_t inode=FindNode(p)-fNDataNodes+1;
+    if(GetNodeInfo(inode)->Has(p))return inode;
+
+    // Search extra nodes
+
+    for(inode=fNDataNodes;inode<GetNTNodes();inode++){
+       if(GetNodeInfo(inode)->Has(p)){return inode;}
+    }
+
+    // Search for nearest neighbor
+    Float_t dist;
+    Float_t closestdist=10000;
+    inode=-1;
+    for(Int_t ii=0;ii<GetNTNodes();ii++){
+       AliTRDTKDNodeInfo *node=GetNodeInfo(ii);
+       dist=0;
+       for(Int_t idim=0;idim<fNDim;idim++){
+           dist+=TMath::Power((node->fData[idim]-p[idim]),2);
+       }
+       dist=TMath::Sqrt(dist);
+       if(dist<closestdist){closestdist=dist;inode=ii;}
+    }
+    return inode;
+}
+
+//_________________________________________________________________
+AliTRDTKDInterpolator::AliTRDTKDNodeInfo* AliTRDTKDInterpolator::GetNodeInfo(Int_t inode) const
+{
+  if(!fNodes || inode >= GetNTNodes()) return NULL;
+  return (AliTRDTKDNodeInfo*)(*fNodes)[inode];
+}
+
+//_________________________________________________________________
+Int_t AliTRDTKDInterpolator::GetNTNodes() const 
+{
+  return fNodes?fNodes->GetEntriesFast():0;
+}
+
+//_________________________________________________________________
+Bool_t AliTRDTKDInterpolator::GetRange(Int_t idim,Float_t range[2]) const
+{
+    if(!fNodes) return kFALSE;
+    if(idim<0 || idim>=fNDim){
+    range[0]=0.; range[1]=0.;
+    return kFALSE;
+    }
+    range[0]=1.e10; range[1]=-1.e10;
+    for(Int_t in=GetNTNodes(); in--; ){
+       AliTRDTKDNodeInfo *node = GetNodeInfo(in);
+
+       if(node->fBounds[2*idim]<range[0]) range[0] = node->fBounds[2*idim];
+       if(node->fBounds[2*idim+1]>range[1]) range[1] = node->fBounds[2*idim+1];
+    }
+
+    return kTRUE;
+}
+
+//_________________________________________________________________
+TH2Poly *AliTRDTKDInterpolator::Projection(Int_t xdim,Int_t ydim)
+{
+    Float_t rangex[2],rangey[2];
+    GetRange(xdim,rangex);
+    GetRange(ydim,rangey);
+
+    TH2Poly* h2 = new TH2Poly("hTKDnodes","hTKDnodes", rangex[0],rangex[1],rangey[0],rangey[1]);
+    h2->GetXaxis()->SetTitle(Form("Q_{%d}", xdim));
+    h2->GetYaxis()->SetTitle(Form("Q_{%d}", ydim));
+
+    for(Int_t inode=0;inode<GetNTNodes();inode++){
+
+       AliTRDTKDNodeInfo* node=GetNodeInfo(inode);
+       h2->AddBin(node->fBounds[2*xdim],node->fBounds[2*ydim],node->fBounds[2*xdim+1],node->fBounds[2*ydim+1]);
+       h2->SetBinContent(inode+1, node->fVal[0]);
+    }
+    return h2;
+}
+
+//_________________________________________________________________
+void AliTRDTKDInterpolator::BuildInterpolation()
+{
+
+    // Calculate Interpolation
+
+    Double_t buffer[fLambda];
+
+    Float_t **refPoints = new Float_t*[fNDim];
+    for(int id=0; id<fNDim; id++){
+       refPoints[id] = new Float_t[GetNTNodes()];
+       for(int in=0; in<GetNTNodes(); in++) refPoints[id][in] = GetNodeInfo(in)->fData[id];
+    }
+    TKDTreeIF *KDhelper = new TKDTreeIF(GetNTNodes(), fNDim, 30, refPoints);
+    KDhelper->Build();
+    KDhelper->MakeBoundariesExact();
+
+    Float_t dist[fNPointsI];
+    Int_t ind[fNPointsI];
+
+    TLinearFitter fitter(fLambda, Form("hyp%d", fLambda-1));
+
+    Int_t nodeIndex(0); Float_t param[6], *pp(NULL);
+    nodeIndex=GetNTNodes(); pp=&param[0];
+    while(nodeIndex--){
+
+       fitter.ClearPoints();
+
+       AliTRDTKDNodeInfo *node = GetNodeInfo(nodeIndex);
+       // find nearest neighbors
+       KDhelper->FindNearestNeighbors(node->fData,fNPointsI, &ind[0], &dist[0]);
+
+       for(int in=0;in<fNPointsI;in++){
+           AliTRDTKDNodeInfo *nnode = GetNodeInfo(ind[in]);
+
+            Float_t w=1; //weight
+           // calculate tri-cubic weighting function
+           if(fUseWeights){
+               Float_t d = dist[in]/dist[fNPointsI-1];
+               Float_t w0 = (1. - d*d*d);
+               w = w0*w0*w0;
+               if(w<1.e-30) continue;
+           }
+           Int_t ipar=0;
+           for(int idim=0; idim<fNDim; idim++){
+               buffer[ipar++] = nnode->fData[idim];
+               for(int jdim=idim; jdim<fNDim; jdim++) buffer[ipar++] = nnode->fData[idim]*nnode->fData[jdim];
+           }
+           fitter.AddPoint(buffer,nnode->fVal[0], nnode->fVal[1]/w);
+
+           // Ensure Boundary Condition
+           for(Int_t kdim=0;kdim<fNDim;kdim++){
+               if(node->fBounds[2*kdim]==0){
+                   Float_t zdata[fNDim];
+                    memcpy(&zdata[0],node->fData,fNDim*sizeof(Float_t));
+                   zdata[kdim]=0;
+                   ipar=0;
+                   for(int idim=0; idim<fNDim; idim++){
+                       buffer[ipar++] = zdata[idim];
+                       for(int jdim=idim; jdim<fNDim; jdim++) buffer[ipar++] = zdata[idim]*zdata[jdim];
+                   }
+                   fitter.AddPoint(buffer,0,1);
+               }
+           }
+       }
+
+       fitter.Eval();
+
+       // retrive fitter results
+       TMatrixD cov(fLambda, fLambda);
+       TVectorD par(fLambda);
+       fitter.GetCovarianceMatrix(cov);
+       fitter.GetParameters(par);
+
+       // store results
+       node->Store(&par,&cov);
+    }
+
+    delete KDhelper;
+    for(int id=0; id<fNDim; id++){
+       delete refPoints[id];
+    }
+    delete[] refPoints;
+}
+
+//_________________________________________________________________
+void AliTRDTKDInterpolator::BuildBoundaryNodes(){
+
+    Int_t nnew=0;
+
+    Float_t treebounds[2*fNDim];
+    for(Int_t idim=0;idim<fNDim;idim++){
+       GetRange(idim,&treebounds[2*idim]);
+    }
+
+    for(int inode=0; inode<GetNTNodes(); inode++){
+
+       AliTRDTKDNodeInfo *node=GetNodeInfo(inode);
+
+       for(Int_t vdim=0;vdim<fNDim;vdim++){
+
+           // Try expansion to lower and higher values
+           for(Int_t iter=0;iter<2;iter++){
+               if(node->fBounds[2*vdim+iter]==treebounds[2*vdim+iter]){
+
+                   // Add new Node
+                   new ((*fNodes)[GetNTNodes()]) AliTRDTKDNodeInfo(fNDim);
+
+                   AliTRDTKDNodeInfo *newnode = GetNodeInfo(GetNTNodes()-1);
+                   if(iter==0)newnode->fBounds[2*vdim+iter]=0;
+                    if(iter==1)newnode->fBounds[2*vdim+iter]=2*treebounds[2*vdim+iter];
+                   newnode->fBounds[2*vdim+!iter]=node->fBounds[2*vdim+iter];
+                   for(Int_t idim=0;idim<fNDim;idim++){
+                       if(idim==vdim)continue;
+                       newnode->fBounds[2*idim]=node->fBounds[2*idim];
+                       newnode->fBounds[2*idim+1]=node->fBounds[2*idim+1];
+                   }
+                   newnode->fVal[0]=0;
+                   newnode->fVal[1]=Float_t(1)/fNPoints;
+                   for(int idim=0; idim<fNDim; idim++){
+                       newnode->fVal[1] /= (newnode->fBounds[2*idim+1] - newnode->fBounds[2*idim]);
+                        newnode->fData[idim]=0.5*(newnode->fBounds[2*idim+1] + newnode->fBounds[2*idim]);
+                   }
+                   nnew++;
+               }
+           }
+       }
+    }
+    printf("%d Boundary Nodes Added \n",nnew);
+}
+
+//_________________________________________________________________
+AliTRDTKDInterpolator::AliTRDTKDNodeInfo::AliTRDTKDNodeInfo(Int_t ndim):
+  TObject()
+  ,fNDim(ndim)
+  ,fNBounds(2*ndim)
+  ,fNPar(1 + ndim + (ndim*(ndim+1)>>1))
+  ,fNCov(fNPar*fNPar)
+  ,fData(NULL)
+  ,fBounds(NULL)
+  ,fPar(NULL)
+  ,fCov(NULL)
+{
+  // Default constructor
+    fVal[0] = 0.; fVal[1] = 0.;
+    fData=new Float_t[fNDim];
+    fBounds=new Float_t[fNBounds];
+}
+
+//_________________________________________________________________
+AliTRDTKDInterpolator::AliTRDTKDNodeInfo::AliTRDTKDNodeInfo(const AliTRDTKDInterpolator::AliTRDTKDNodeInfo &ref):
+TObject(ref)
+   ,fNDim(ref.fNDim)
+   ,fNBounds(ref.fNBounds)
+   ,fNPar(ref.fNPar)
+   ,fNCov(ref.fNCov)
+   ,fData(NULL)
+   ,fBounds(NULL)
+   ,fPar(NULL)
+   ,fCov(NULL)
+{
+  // Copy constructor
+
+    if(ref.fData){
+       fData = new Float_t[fNDim];
+       memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
+    }
+    if(ref.fBounds){
+       fBounds = new Float_t[2*fNDim];
+       memcpy(fBounds, ref.fBounds, 2*fNDim*sizeof(Float_t));
+    }
+
+    fVal[0] = ref.fVal[0];
+    fVal[1] = ref.fVal[1];
+
+    if(ref.fPar){
+       fPar=new Double_t[fNPar];
+       memcpy(fPar, ref.fPar, fNPar*sizeof(Double_t));
+    }
+    if(ref.fCov){
+       fCov=new Double_t[fNCov];
+       memcpy(fCov, ref.fCov, fNCov*sizeof(Double_t));
+    }
+}
+
+//_________________________________________________________________
+AliTRDTKDInterpolator::AliTRDTKDNodeInfo::~AliTRDTKDNodeInfo()
+{
+  // Destructor
+  if(fData) delete [] fData;
+  if(fBounds) delete [] fBounds;
+  if(fPar) delete [] fPar;
+  if(fCov) delete [] fCov;
+}
+
+//_________________________________________________________________
+AliTRDTKDInterpolator::AliTRDTKDNodeInfo &AliTRDTKDInterpolator::AliTRDTKDNodeInfo::operator=(const AliTRDTKDInterpolator::AliTRDTKDNodeInfo &ref)
+{
+    //
+    // Assignment operator
+    //
+
+    if(&ref==this) return *this;
+    memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
+    memcpy(fBounds, ref.fBounds, 2*fNDim*sizeof(Float_t));
+
+    fVal[0] = ref.fVal[0];
+    fVal[1] = ref.fVal[1];
+
+    if(ref.fPar){
+       fPar=new Double_t[fNPar];
+       memcpy(fPar, ref.fPar, fNPar*sizeof(Double_t));
+    }
+    if(ref.fCov){
+       fCov=new Double_t[fNCov];
+       memcpy(fCov, ref.fCov, fNCov*sizeof(Double_t));
+    }
+    return *this;
+}
+
+
+//_________________________________________________________________
+void AliTRDTKDInterpolator::AliTRDTKDNodeInfo::Print(const Option_t *opt) const
+{
+  // Print the content of the node
+  printf("x [");
+  for(int idim=0; idim<fNDim; idim++) printf("%f ", fData?fData[idim]:0.);
+  printf("] f = [%e +- %e]\n", fVal[0], fVal[1]);
+
+  if(fBounds){
+      printf("range [");
+      for(int idim=0; idim<fNDim; idim++) printf("(%f %f) ", fBounds[2*idim], fBounds[2*idim+1]);
+    printf("]\n");
+  }
+  if(strcmp(opt, "a")!=0) return;
+
+  if(fPar){
+    printf("Fit parameters : \n");
+    for(int ip=0; ip<fNPar; ip++) printf("p%d[%f] ", ip, fPar[ip]);
+    printf("\n");
+  }
+  if(!fCov) return;
+  for(int ip(0), n(0); ip<fNPar; ip++){
+    for(int jp(ip); jp<fNPar; jp++) printf("c(%d %d)[%f] ", ip, jp, fCov[n++]);
+    printf("\n");
+  }
+}
+
+//_________________________________________________________________
+void AliTRDTKDInterpolator::AliTRDTKDNodeInfo::Store(TVectorD const *par, TMatrixD const *cov)
+{
+// Store the parameters and the covariance in the node
+
+    if(!fPar){fPar = new Double_t[fNPar];}
+    for(int ip=0; ip<fNPar; ip++) fPar[ip] = (*par)[ip];
+
+    if(!cov) return;
+    if(!fCov){fCov = new Double_t[fNCov];}
+    for(int ip(0), np(0); ip<fNPar; ip++)
+       for(int jp=ip; jp<fNPar; jp++) fCov[np++] = (*cov)(ip,jp);
+}
+
+//_________________________________________________________________
+Bool_t AliTRDTKDInterpolator::AliTRDTKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t &error) const
+{
+    // Recalculate the PDF for one node from the results of interpolation (parameters and covariance matrix)
+
+    result =0.; error = 1.;
+    if(!fPar) return kFALSE;
+
+    Double_t fdfdp[fNDim];
+    Int_t ipar = 0;
+    fdfdp[ipar++] = 1.;
+    for(int idim=0; idim<fNDim; idim++){
+       fdfdp[ipar++] = point[idim];
+       for(int jdim=idim; jdim<fNDim; jdim++) fdfdp[ipar++] = point[idim]*point[jdim];
+    }
+
+    // calculate estimation
+    for(int i=0; i<fNPar; i++) result += fdfdp[i]*fPar[i];
+
+    if(!fCov)return kTRUE;
+    // calculate error
+    error=0;
+
+    for(int i(0), n(0); i<fNPar; i++){
+       error += fdfdp[i]*fdfdp[i]*fCov[n++];
+       for(int j(i+1); j<fNPar; j++) error += 2.*fdfdp[i]*fdfdp[j]*fCov[n++];
+    }
+    if(error>0)error = TMath::Sqrt(error);
+    else{error=100;}
+
+    // Boundary condition
+    if(result<0){
+       result=fVal[0];
+       error=fVal[1];
+    }
+
+    return kTRUE;
+}
+
+//_____________________________________________________________________
+Bool_t AliTRDTKDInterpolator::AliTRDTKDNodeInfo::Has(const Float_t *p) const
+{
+  Int_t n(0);
+  for(int id=0; id<fNDim; id++)
+      if(p[id]>=fBounds[2*id] && p[id]<fBounds[2*id+1]) n++;
+  if(n==fNDim) return kTRUE;
+  return kFALSE;
+}
+
diff --git a/STEER/STEERBase/AliTRDTKDInterpolator.h b/STEER/STEERBase/AliTRDTKDInterpolator.h
new file mode 100644 (file)
index 0000000..e930f36
--- /dev/null
@@ -0,0 +1,104 @@
+#ifndef ROOT_ALITRDTKDINTERPOLATOR_H
+#define ROOT_ALITRDTKDINTERPOLATOR_H
+
+#ifndef ROOT_TKDTree
+#include "TKDTree.h"
+#endif
+
+#ifndef ROOT_TObject
+#include "TObject.h"
+#endif
+
+#include "TVectorD.h"
+#include "TMatrixD.h"
+#include "TH2Poly.h"
+class TClonesArray;
+
+class AliTRDTKDInterpolator : public TKDTreeIF
+{
+
+
+
+  
+public:
+    enum TRDTKDMode{
+       kInterpolation=0,
+       kMinError=1,
+       kNodeVal=2
+    };
+
+      // Bucket Object class
+    class AliTRDTKDNodeInfo : public TObject
+    {
+        friend class AliTRDTKDInterpolator;
+    public:
+
+       AliTRDTKDNodeInfo(Int_t ndim = 0);
+       AliTRDTKDNodeInfo(const AliTRDTKDInterpolator::AliTRDTKDNodeInfo &ref);
+       AliTRDTKDInterpolator::AliTRDTKDNodeInfo& operator=(const AliTRDTKDInterpolator::AliTRDTKDNodeInfo &ref);
+       virtual ~AliTRDTKDNodeInfo();
+       Bool_t        CookPDF(const Double_t *point, Double_t &result, Double_t &error) const;
+       Bool_t Has(const Float_t *p) const;
+       void          Print(const Option_t * = "") const;
+       void          Store(TVectorD const *par, TMatrixD const *cov=NULL);
+
+    private:
+       Int_t fNDim;              // Dimension of Points
+        Int_t fNBounds;           // 2* Dimension of Points
+       Int_t fNPar;            // Number of Parameters
+       Int_t fNCov;            // Size of Cov Matrix
+       Float_t   *fData;         //[fNDim] Data Point
+       Float_t *fBounds;         //[fNBounds] Boundaries
+       Float_t   fVal[2];        //measured value for node
+       Double_t  *fPar;          //[fNPar] interpolator parameters
+       Double_t  *fCov;          //[fNCov] interpolator covariance matrix
+
+       ClassDef(AliTRDTKDNodeInfo, 1)  // node info for interpolator
+    };
+
+public:
+
+    AliTRDTKDInterpolator();
+    AliTRDTKDInterpolator(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data);
+    ~AliTRDTKDInterpolator();
+
+    Bool_t        Eval(const Double_t *point, Double_t &result, Double_t &error);
+    void          Print(const Option_t *opt="") const;
+
+    TH2Poly *     Projection(Int_t xdim,Int_t ydim);
+
+    Int_t         GetNDIM() const {return fNDim;}
+    Bool_t        GetRange(Int_t idim,Float_t range[2]) const;
+  
+    void          SetNPointsInterpolation(Int_t np){fNPointsI=np;};
+    Int_t         GetNPointsInterpolation(){return fNPointsI;};
+
+    void          SetUseWeights(Bool_t k=kTRUE){fUseWeights=k;}
+    void          SetPDFMode(TRDTKDMode mod){fPDFMode=mod;}
+
+    Bool_t        Build();
+
+    void          SetUseHelperNodes(Bool_t k){fUseHelperNodes=k;}
+
+private:
+
+    Int_t         GetNodeIndex(const Float_t *p);
+    AliTRDTKDInterpolator::AliTRDTKDNodeInfo*  GetNodeInfo(Int_t inode) const;
+    Int_t         GetNTNodes() const;
+    void          BuildInterpolation();
+    void          BuildBoundaryNodes();
+    AliTRDTKDInterpolator(const AliTRDTKDInterpolator &ref);
+    AliTRDTKDInterpolator &operator=(const AliTRDTKDInterpolator &ref);
+
+    Int_t fNDataNodes;         // Number of filled nodes (total-zero nodes)
+    TClonesArray  *fNodes;     //interpolation nodes
+    UChar_t       fLambda;      // number of parameters in polynom
+    Int_t         fNPointsI;    // number of points for interpolation
+    Bool_t        fUseHelperNodes; // Build Helper nodes to ensure boundary conditions
+    Bool_t fUseWeights; // Use tricubic weights
+    TRDTKDMode    fPDFMode; // Mode for PDF calculation
+
+    ClassDef(AliTRDTKDInterpolator, 1)   // data interpolator based on KD tree
+};
+
+#endif
index 2becb10..e9906bd 100644 (file)
@@ -87,6 +87,9 @@
 #pragma link C++ class AliTRDPIDReference+;
 #pragma link C++ class AliTRDPIDParams+;
 #pragma link C++ class AliTRDPIDParams::AliTRDPIDThresholds+;
+#pragma link C++ class AliTRDPIDResponseObject+;
+#pragma link C++ class AliTRDTKDInterpolator+;
+#pragma link C++ class AliTRDTKDInterpolator::AliTRDTKDNodeInfo+;
 #pragma link C++ class AliITSPidParams+;
 #pragma link C++ class AliPIDResponse+;
 #pragma link C++ class AliITSPIDResponse+;
index 8e501ee..dcc2e59 100644 (file)
@@ -1544,7 +1544,7 @@ AliTRDPIDResponse *AliTRDcalibDB::GetPIDResponse(AliTRDPIDResponse::ETRDPIDMetho
     Bool_t hasReference = kFALSE;
     while ((obj = refs())){
       if ((ref = dynamic_cast<AliTRDPIDReference *>(obj))){
-        fPIDResponse->Load(ref);
+        //fPIDResponse->Load(ref);
         hasReference = kTRUE;
         break;
       }