Update of the TRD PID Response:
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 3 May 2011 10:14:01 +0000 (10:14 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 3 May 2011 10:14:01 +0000 (10:14 +0000)
- New reference distributions covering a momentum range from 1 to 10 GeV/c with an additional step at 1.5GeV/c. The reference in the OCDB up to now are going up to 6 GeV/c in integer momentum steps. With the change we increase the pt reach for the TRD up to 10GeV/c.

- The container type has changed: Before the references were stored in a TObjArray and associated to momentum steps and species  via the name. Momentum steps were hardcoded in the PID Response, adapted to the old TRD PID code. With the new references having a different momentum binning we have seen that this is difficult to handle. A new class AliTRDPIDReference stores the reference objects together with the reference object container. The new container can handle references for all TRD PID methods.

- To keep backward compatibility the OCDB object contains the TObjArray containing the old reference objects and the reference container with the new references. In the code currently implemented the Load function checks the name of the TObject, and if it is not valid it will be neglected. In the new code we explicitly expect the new container, and all objects of a different type will be ignored. So backward compatibility is not an issue. One file in the TRD code is affected by the change, which has to be updated the same time  commit is done.

- The TRD tender was also adapted to the change and the configuration was updated.

Concerning the OCDB entry since the TRD PID Response is now used in the reconstruction to calculate the TRD PID Probabilities we were concerned that we would cause problems there if we provide only the corresponding object for the new code while the reconstruction is still running with the old code. For this since the current OCDB entry stores the reference object in a simple TObjArray, we considered to be compatible with both versions of the code we add the new object to the old references. Since in the old code references which are not compatible with the name convention for the old objects are rejected, and in the new code we require the reference object to be of the new type, both versions should work with the new code.

Markus

12 files changed:
ANALYSIS/CMakelibTENDERSupplies.pkg
ANALYSIS/TenderSupplies/AddTaskTender.C
ANALYSIS/TenderSupplies/AliTRDTenderSupply.cxx
ANALYSIS/TenderSupplies/AliTRDTenderSupply.h
OCDB/TRD/Calib/PIDLQ1D/Run0_999999999_v1_s0.root [new file with mode: 0644]
STEER/AliTRDPIDReference.cxx [new file with mode: 0644]
STEER/AliTRDPIDReference.h [new file with mode: 0644]
STEER/AliTRDPIDResponse.cxx
STEER/AliTRDPIDResponse.h
STEER/CMakelibSTEERBase.pkg
STEER/STEERBaseLinkDef.h
TRD/AliTRDcalibDB.cxx

index a41a6d8..7672ede 100644 (file)
@@ -35,6 +35,6 @@ set ( EINCLUDE  ANALYSIS ANALYSIS/Tender STEER TOF TRD/Cal VZERO ANALYSIS/Tender
 
 if( ALICE_TARGET STREQUAL "win32gcc")
        
-               set ( PACKSOFLAGS  ${SOFLAGS} -L${ALICE_ROOT}/lib/tgt_${ALICE_TARGET} -lSTEERBase -lSTEER -lANALYSIS -lCDB -lT0base -lT0rec -lTPCbase)
+               set ( PACKSOFLAGS  ${SOFLAGS} -L${ALICE_ROOT}/lib/tgt_${ALICE_TARGET} -lSTEERBase -lSTEER -lANALYSIS -lCDB -lT0base -lT0rec -lTPCbase -lTRDbase)
 
 endif( ALICE_TARGET STREQUAL "win32gcc")
index 6f6e0d9..fb6706b 100644 (file)
@@ -43,6 +43,27 @@ AliAnalysisTask *AddTaskTender(Bool_t useV0=kFALSE){
   
   //========= Attach TRD supply ======
   AliTRDTenderSupply *trdSupply=new AliTRDTenderSupply("TRDtender");
+  trdSupply->SetLoadReferencesFromCDB();
+  // Mask Bad chambers
+  trdSupply->AddBadChamber(265);      // low drift
+  trdSupply->AddBadChamber(50);
+  trdSupply->AddBadChamber(524);      // low drift
+  trdSupply->AddBadChamber(32);       // intermediate gain
+  trdSupply->AddBadChamber(15);
+  trdSupply->AddBadChamber(231);      // low gain
+  trdSupply->AddBadChamber(273);      // intermediate gain
+  trdSupply->AddBadChamber(532);
+  trdSupply->AddBadChamber(5);        // low drift
+  trdSupply->AddBadChamber(227);
+  trdSupply->AddBadChamber(287);      // low drift
+  trdSupply->AddBadChamber(212);      // intermediate gain
+  trdSupply->AddBadChamber(228);      // low gain
+  trdSupply->AddBadChamber(52);       // low gain
+  trdSupply->AddBadChamber(169);      // low drift
+  trdSupply->AddBadChamber(236);      // low drift
+
+  trdSupply->SetPIDmethod(AliTRDTenderSupply::k1DLQpid);
+  trdSupply->SetNormalizationFactor(1./7.603);
   tender->AddSupply(trdSupply);
 
   //========= Attach PID supply ======
index 7de68ec..3d1fcbd 100644 (file)
@@ -41,6 +41,7 @@
 #include <AliESDtrack.h>
 #include <AliESDInputHandler.h>
 #include <AliAnalysisManager.h>
+#include <AliTRDPIDReference.h>
 #include <AliTRDPIDResponse.h>
 #include <AliTender.h>
 
@@ -57,16 +58,18 @@ AliTRDTenderSupply::AliTRDTenderSupply() :
   fChamberGainNew(NULL),
   fChamberVdriftOld(NULL),
   fChamberVdriftNew(NULL),
-  fFileNNref(""),
+  fRefFilename(""),
   fPIDmethod(kNNpid),
+  fNormalizationFactor(1.),
   fPthreshold(0.8),
   fNBadChambers(0),
-  fGainCorrection(kTRUE)
+  fGainCorrection(kTRUE),
+  fLoadReferencesFromCDB(kFALSE),
+  fHasReferences(kFALSE)
 {
   //
   // default ctor
   //
-  memset(fBadChamberID, 0, sizeof(Int_t) * kNChambers);
 }
 
 //_____________________________________________________
@@ -78,11 +81,14 @@ AliTRDTenderSupply::AliTRDTenderSupply(const char *name, const AliTender *tender
   fChamberGainNew(NULL),
   fChamberVdriftOld(NULL),
   fChamberVdriftNew(NULL),
-  fFileNNref(""),
+  fRefFilename(""),
   fPIDmethod(kNNpid),
+  fNormalizationFactor(1.),
   fPthreshold(0.8),
   fNBadChambers(0),
-  fGainCorrection(kTRUE)
+  fGainCorrection(kTRUE),
+  fLoadReferencesFromCDB(kFALSE),
+  fHasReferences(kFALSE)
 {
   //
   // named ctor
@@ -113,6 +119,10 @@ void AliTRDTenderSupply::Init()
     fESDpid=new AliESDpid;
     fTender->GetESDhandler()->SetESDpid(fESDpid);
   }
+  // Load References
+  if(!fLoadReferencesFromCDB) LoadReferences();
+  fESDpid->GetTRDResponse().SetGainNormalisationFactor(fNormalizationFactor);
+
   // Set Normalisation Factors
   if(mgr->GetMCtruthEventHandler()){
     // Assume MC
@@ -136,6 +146,7 @@ void AliTRDTenderSupply::ProcessEvent()
   if (fTender->RunChanged()){
     AliDebug(0, Form("AliTPCTenderSupply::ProcessEvent - Run Changed (%d)\n",fTender->GetRun()));
     if (fGainCorrection) SetChamberGain();
+    if(!fHasReferences) LoadReferences();
   }
 
   fESD = fTender->GetEvent();
@@ -162,6 +173,51 @@ void AliTRDTenderSupply::ProcessEvent()
   }
 }
 
+//_____________________________________________________
+void AliTRDTenderSupply::LoadReferences(){
+  //
+  // Load Reference from the OCDB/OADB into the PID Response
+  //
+  if(fLoadReferencesFromCDB){
+    AliDebug(1, "Loading Reference Distributions from the OCDB");
+    AliCDBEntry *en = fTender->GetCDBManager()->Get("TRD/Calib/PIDLQ1D");
+    if(!en){
+      AliError("References for 1D Likelihood Method not in OCDB");
+      return;
+    }
+    en->GetId().Print();
+    TObjArray *arr = dynamic_cast<TObjArray *>(en->GetObject());
+    if(!arr) AliError("List with the references not found");
+  
+    // Get new references
+    TIter refs(arr);
+    TObject *o = NULL;
+    AliTRDPIDReference *ref = NULL;
+    while((o = refs())){
+      if(!TString(o->IsA()->GetName()).CompareTo("AliTRDPIDReference")){
+        ref = dynamic_cast<AliTRDPIDReference *>(o);
+        break;
+      }
+    }
+    if(ref){
+      fESDpid->GetTRDResponse().Load(ref);
+      fHasReferences = kTRUE;
+      AliInfo("Reference distributions loaded into the PID Response");
+    } else {
+      AliError("References not found");
+    }
+  } else {
+    // Backward compatibility mode
+    AliInfo("Loading Reference Distributions from ROOT file");
+    if(fRefFilename.Length() != 0){
+       fESDpid->GetTRDResponse().Load(fRefFilename.Data());
+       fHasReferences = kTRUE;
+    } else{
+       AliError("No file defined");
+    }
+  }
+}
+
 //_____________________________________________________
 void AliTRDTenderSupply::SetChamberGain(){
   //
@@ -169,10 +225,7 @@ void AliTRDTenderSupply::SetChamberGain(){
   //
   
   //find previous entry from the UserInfo
-  //   TTree *tree=((TChain*)fTender->GetInputData(0))->GetTree();
-  AliAnalysisManager*mgr = AliAnalysisManager::GetAnalysisManager();
-  AliAnalysisTaskSE *task = (AliAnalysisTaskSE*)mgr->GetTasks()->First();
-  TTree *tree=((TChain*)task->GetInputData(0))->GetTree();
+  TTree *tree=((TChain*)fTender->GetInputData(0))->GetTree();
   if (!tree) {
   AliError("Tree not found in ESDhandler");
     return;
index 5eb9543..dc380ab 100644 (file)
@@ -14,8 +14,8 @@
 
 #include <AliTenderSupply.h>
 
+class AliTRDpidRecalculator;
 class AliTRDCalDet;
-class AliESDtrack;
 class AliESDEvent;
 
 class AliTRDTenderSupply: public AliTenderSupply {
@@ -30,8 +30,10 @@ public:
   AliTRDTenderSupply(const char *name, const AliTender *tender=NULL);
   virtual ~AliTRDTenderSupply();
 
-  void SetNNref(const char* file) {fFileNNref=file;}
+  void SetLoadReferencesFromCDB() { fLoadReferencesFromCDB = kTRUE; }
+  void SetRefFile(const char* file) {fRefFilename=file;}
   void SetPIDmethod(Int_t pidMethod) { fPIDmethod = pidMethod; }
+  void SetNormalizationFactor(Double_t norm) { fNormalizationFactor = norm; }
   void SetCalibLowpThreshold(Double_t pmin) { fPthreshold = pmin; };
   
   virtual void              Init();
@@ -51,6 +53,7 @@ private:
   Bool_t GetTRDchamberID(AliESDtrack * const track, Int_t *detectors);
   void SetChamberGain();
   void ApplyGainCorrection(AliESDtrack *track);
+  void LoadReferences();
   
   AliESDEvent           *fESD;       //! the ESD Event
   AliESDpid             *fESDpid;    //! ESD PID object
@@ -60,12 +63,15 @@ private:
   AliTRDCalDet *fChamberVdriftOld;   // Old drift velocity calibration
   AliTRDCalDet *fChamberVdriftNew;   // New drift velocity calibration
 
-  TString fFileNNref;                // path and name to the NNref file
+  TString fRefFilename;                // path and name to the NNref file
   Int_t fPIDmethod;                  // PID method
+  Double_t fNormalizationFactor;     // dE/dx Normalization Factor 
   Double_t fPthreshold;              // Low Momentum threshold for calibration
   Int_t fBadChamberID[kNChambers];   // List of Bad Chambers
   UInt_t fNBadChambers;              // Number of bad chambers
   Bool_t fGainCorrection;            // Apply gain correction 
+  Bool_t fLoadReferencesFromCDB;     // Load References from CDB
+  Bool_t fHasReferences;             // has references loaded
   
   AliTRDTenderSupply(const AliTRDTenderSupply&c);
   AliTRDTenderSupply& operator= (const AliTRDTenderSupply&c);
diff --git a/OCDB/TRD/Calib/PIDLQ1D/Run0_999999999_v1_s0.root b/OCDB/TRD/Calib/PIDLQ1D/Run0_999999999_v1_s0.root
new file mode 100644 (file)
index 0000000..f09dead
Binary files /dev/null and b/OCDB/TRD/Calib/PIDLQ1D/Run0_999999999_v1_s0.root differ
diff --git a/STEER/AliTRDPIDReference.cxx b/STEER/AliTRDPIDReference.cxx
new file mode 100644 (file)
index 0000000..e3e6379
--- /dev/null
@@ -0,0 +1,199 @@
+/**************************************************************************
+ * 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>
+//
+#include "TObjArray.h"
+
+#include "AliLog.h"
+
+#include "AliTRDPIDReference.h"
+
+ClassImp(AliTRDPIDReference)
+
+//____________________________________________________________
+AliTRDPIDReference::AliTRDPIDReference():
+TNamed(),
+fRefContainer(NULL),
+fMomentumBins()
+{
+       //
+       // Dummy constructor
+       //
+       SetBit(kIsOwner, kTRUE);
+}
+
+//____________________________________________________________
+AliTRDPIDReference::AliTRDPIDReference(const Char_t *name):
+                 TNamed(name, "TRD PID References"),
+                 fRefContainer(NULL),
+                 fMomentumBins()
+{
+       //
+       // Default constructor
+       //
+       SetBit(kIsOwner, kTRUE);
+}
+
+//____________________________________________________________
+AliTRDPIDReference::AliTRDPIDReference(const AliTRDPIDReference &ref):
+                 TNamed(ref),
+                 fRefContainer(ref.fRefContainer),
+                 fMomentumBins(ref.fMomentumBins)
+{
+       //
+       // Copy constructor
+       // Only copies poiters, object is not the owner of the references
+       //
+       SetBit(kIsOwner, kFALSE);
+}
+
+//____________________________________________________________
+AliTRDPIDReference &AliTRDPIDReference::operator=(const AliTRDPIDReference &ref){
+       //
+       // Assginment operator
+       // Only copies poiters, object is not the owner of the references
+       //
+       if(this != &ref){
+               TNamed::operator=(ref);
+               if(TestBit(kIsOwner) && fRefContainer) delete fRefContainer;
+               fRefContainer = ref.fRefContainer;
+               fMomentumBins = ref.fMomentumBins;
+               SetBit(kIsOwner, kFALSE);
+       }
+       return *this;
+}
+
+//____________________________________________________________
+AliTRDPIDReference::~AliTRDPIDReference(){
+       //
+       // Destructor
+       // references are deleted if the object is the owner
+       //
+       if(fRefContainer && TestBit(kIsOwner)) delete fRefContainer;
+}
+
+//____________________________________________________________
+void AliTRDPIDReference::SetNumberOfMomentumBins(Int_t nBins, Float_t *momenta){
+       //
+       // Set the momentum binning
+       //
+       if(fRefContainer) fRefContainer->Clear();
+       else{
+               fRefContainer = new TObjArray;
+               fRefContainer->SetOwner();
+       }
+       fRefContainer->Expand(nBins * AliPID::kSPECIES);
+       fMomentumBins.Set(nBins,momenta);
+}
+
+//____________________________________________________________
+void AliTRDPIDReference::AddReference(TObject *ref, AliPID::EParticleType spec, Int_t pbin){
+       //
+       // Add a new reference distribution for a given species and a givem
+       // momentum value to the reference container
+       // The reference distribution is associated with the momentum value defined for the
+       // given momentum bin
+       //
+       if(!fRefContainer){
+               AliError("Reference Container not initialized");
+               return;
+       }
+       if(pbin > fMomentumBins.GetSize()){
+               AliError("Pbin overflow");
+               return;
+       }
+       AliDebug(1, Form("Adding object with address %p to position %d", ref, spec * fMomentumBins.GetSize() + pbin));
+       fRefContainer->AddAt(ref, spec * fMomentumBins.GetSize() + pbin);
+}
+
+//____________________________________________________________
+TObject *AliTRDPIDReference::GetUpperReference(AliPID::EParticleType spec, Float_t p, Float_t &pUpper) const{
+       //
+       // Get the next reference associated with a momentum larger than the requested momentum
+       // In case no next upper reference is found, NULL is returned
+       // The momentum value the reference is associated to is stored in the reference pUpper
+       //
+       Int_t bin = -1;
+       pUpper = 20;
+       for(Int_t ip = 0; ip < fMomentumBins.GetSize(); ip++){
+               AliDebug(10, Form("Bin %d, p = %.1f", ip, fMomentumBins[ip]));
+               if(p < fMomentumBins[ip]){
+                       AliDebug(10, "Bin found");
+                       bin = ip;
+                       break;
+               }
+       }
+       AliDebug(2, Form("p = %.1f, bin = %d\n", p, bin));
+       if(bin >= 0) {
+               pUpper = fMomentumBins[bin];
+               return fRefContainer->At(spec * fMomentumBins.GetSize() + bin);
+       }
+       else return NULL;
+}
+
+//____________________________________________________________
+TObject *AliTRDPIDReference::GetLowerReference(AliPID::EParticleType spec, Float_t p, Float_t &pLower) const{
+       //
+       // Get the next reference associated with a momentum smaller than the requested momentum
+       // In case no next lower reference is found, NULL is returned
+       // The momentum value the reference is associated to is stored in the reference pLower
+       //
+       Int_t bin = -1;
+       pLower = 0;
+       for(Int_t ip = fMomentumBins.GetSize() - 1; ip >= 0; ip--){
+               AliDebug(10, Form("Bin %d, p = %.1f", ip, fMomentumBins[ip]));
+               if(p > fMomentumBins[ip]){
+                       AliDebug(10, "Bin found");
+                       bin = ip;
+                       break;
+               }
+       }
+       AliDebug(2, Form("p = %.1f, bin = %d\n", p, bin));
+       if(bin >= 0){
+               pLower = fMomentumBins[bin];
+               return fRefContainer->At(spec * fMomentumBins.GetSize() + bin);
+       }
+       else return NULL;
+}
+
+//____________________________________________________________
+void AliTRDPIDReference::Print(const Option_t*) const{
+       //
+       // Print content of the PID reference container
+       //
+       printf("Number of Momentum Bins: %d\n", GetNumberOfMomentumBins());
+       printf("=====================================\n");
+       for(Int_t ip = 0; ip < GetNumberOfMomentumBins(); ip++){
+               printf("Bin %d: p = %.1f\n", ip, fMomentumBins[ip]);
+       }
+       printf("=====================================\n");
+       if(fRefContainer){
+               printf("Content of the reference container:\n");
+               for(Int_t ip = 0; ip < GetNumberOfMomentumBins(); ip++){
+                       printf("[");
+                       for(Int_t is = 0; is < 5; is++) printf("%p|", fRefContainer->At(is * fMomentumBins.GetSize() + ip));
+                       printf("]\n");
+               }
+       }
+}
diff --git a/STEER/AliTRDPIDReference.h b/STEER/AliTRDPIDReference.h
new file mode 100644 (file)
index 0000000..0ac7636
--- /dev/null
@@ -0,0 +1,57 @@
+/**************************************************************************
+ * 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 the reference distributions for the TRD PID
+// Provides storing of the references together with the mometum steps
+// they were taken
+// More information can be found in the implementation file
+//
+#ifndef ALITRDPIDREFERENCE_H
+#define ALITRDPIDREFERENCE_H
+
+#include <TNamed.h>
+#include <TArrayF.h>
+#include "AliPID.h"
+
+class TObjArray;
+
+class AliTRDPIDReference : public TNamed{
+public:
+       AliTRDPIDReference();
+       AliTRDPIDReference(const Char_t *name);
+       AliTRDPIDReference(const AliTRDPIDReference &ref);
+       AliTRDPIDReference &operator=(const AliTRDPIDReference &ref);
+       ~AliTRDPIDReference();
+
+       void SetNumberOfMomentumBins(Int_t nBins, Float_t *momenta);
+       void AddReference(TObject *histo, AliPID::EParticleType spec, Int_t pbin);
+
+       // Derive reference
+       TObject *GetLowerReference(AliPID::EParticleType spec, Float_t p, Float_t &pLower) const;
+       TObject *GetUpperReference(AliPID::EParticleType spec, Float_t p, Float_t &pUpper) const;
+
+       Int_t GetNumberOfMomentumBins() const { return fMomentumBins.GetSize(); }
+       void Print(const Option_t *) const;
+private:
+       enum{
+               kIsOwner = BIT(14)
+       };
+       TObjArray *fRefContainer;     // Histogram container
+       TArrayF fMomentumBins;        // Momentum Bins
+
+       ClassDef(AliTRDPIDReference, 1)         // Container for TRD references
+};
+#endif
+
index 60c182c..3021cb9 100644 (file)
 
 #include "AliLog.h"
 
+#include "AliTRDPIDReference.h"
 #include "AliTRDPIDResponse.h"
 
 ClassImp(AliTRDPIDResponse)
 
-const Double_t AliTRDPIDResponse::fgkPBins[kNPBins] = {1, 2, 3, 4, 5, 6};
-
 //____________________________________________________________
 AliTRDPIDResponse::AliTRDPIDResponse():
   TObject()
-  ,fReferences(NULL)
+  ,fkPIDReference(NULL)
   ,fGainNormalisationFactor(1.)
   ,fPIDmethod(kLQ1D)
 {
   //
   // Default constructor
   //
-  for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++)
-    for(Int_t ipbin = 0; ipbin < kNPBins; ipbin++)
-      fMapRefHists[ispec][ipbin] = -1;
 }
 
 //____________________________________________________________
 AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
   TObject(ref)
-  ,fReferences(ref.fReferences)
+  ,fkPIDReference(ref.fkPIDReference)
   ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
   ,fPIDmethod(ref.fPIDmethod)
 {
   //
   // Copy constructor
-  // Flat copy of the reference histos
-  // For creating a deep copy call SetOwner
   //
-  Int_t size  = (AliPID::kSPECIES)*(kNPBins);
-  memcpy(fMapRefHists, ref.fMapRefHists, sizeof(Int_t) * size);
 }
 
 //____________________________________________________________
 AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
   //
   // Assignment operator
-  // Performs a flat copy of the reference histos
   //
   if(this == &ref) return *this;
   
   // Make copy
   TObject::operator=(ref);
-  if(IsOwner()){
-    if(fReferences){
-      fReferences->Delete();
-      delete fReferences;
-    }
-    SetBit(kIsOwner, kFALSE);
-  } else if(fReferences) delete fReferences;
-  fReferences = ref.fReferences;
-  Int_t size  = (AliPID::kSPECIES)*(kNPBins);
-  memcpy(fMapRefHists, ref.fMapRefHists, sizeof(Int_t) * size);
+  fGainNormalisationFactor = ref.fGainNormalisationFactor;
+  fkPIDReference = ref.fkPIDReference;
   fPIDmethod = ref.fPIDmethod;
   
   return *this;
@@ -102,15 +85,11 @@ AliTRDPIDResponse::~AliTRDPIDResponse(){
   //
   // Destructor
   //
-  if(IsOwner()){
-    // Destroy histos
-    if(fReferences) fReferences->Delete();
-  }
-  if(fReferences) delete fReferences;
+  if(IsOwner()) delete fkPIDReference;
 }
 
 //____________________________________________________________
-Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
+Bool_t AliTRDPIDResponse::Load(const Char_t * filename, const Char_t *refName){
   //
   // Load References into the toolkit
   //
@@ -126,101 +105,16 @@ Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
   }
   
   gROOT->cd();
-  fReferences = new TObjArray(AliPID::kSPECIES*kNPBins);
-  fReferences->SetOwner();
-  TIter iter(in->GetListOfKeys());
-  TKey *histkey = NULL;
-  TObject *tmp = NULL;
-  Int_t arrayPos = 0, pbin = 0; 
-  AliPID::EParticleType species;
-  TString histname;
-  TH1 *hnew = NULL;
-  while((histkey = dynamic_cast<TKey *>(iter()))){
-    tmp = histkey->ReadObj();
-    histname = tmp->GetName();
-    if(histname.BeginsWith("fHQel")){ // Electron histogram
-      histname.ReplaceAll("fHQel_p","");
-      species = AliPID::kElectron;
-    } else if(histname.BeginsWith("fHQmu")){ // Muon histogram
-      histname.ReplaceAll("fHQmu_p","");
-      species = AliPID::kMuon;
-    } else if(histname.BeginsWith("fHQpi")){ // Pion histogram
-      histname.ReplaceAll("fHQpi_p","");
-      species = AliPID::kPion;
-    }else if(histname.BeginsWith("fHQka")){ // Kaon histogram
-      histname.ReplaceAll("fHQka_p","");
-      species = AliPID::kKaon;
-    }else if(histname.BeginsWith("fHQpr")){ // Proton histogram
-      histname.ReplaceAll("fHQpr_p","");
-      species = AliPID::kProton;
-    } else continue;
-    pbin = histname.Atoi() - 1;
-    AliDebug(1, Form("Species %d, Histname %s, Pbin %d, Position in container %d", species, histname.Data(), pbin, arrayPos));
-    fMapRefHists[species][pbin] = arrayPos;
-    fReferences->AddAt((hnew =new TH1F(*dynamic_cast<TH1F *>(tmp))), arrayPos);
-    hnew->SetDirectory(gROOT);
-    arrayPos++;
-  } 
-  in->Close();
+  fkPIDReference = dynamic_cast<const AliTRDPIDReference *>(in->Get(refName)->Clone());
+  in->Close(); delete in;
   owd->cd();
-  AliDebug(2, Form("Successfully loaded %d Reference Histograms", arrayPos));
-  SetBit(kIsOwner, kTRUE);
-  delete in;
-  return kTRUE;
-}
-
-//____________________________________________________________
-Bool_t AliTRDPIDResponse::Load(const TObjArray *histos){
-  //
-  // Load Reference into the PID Response (from OCDB/OADB)
-  //
-  AliDebug(1, "Setting references (from the database)");
-  if(fReferences) fReferences->Clear();
-  else
-    fReferences = new TObjArray(AliPID::kSPECIES*kNPBins);
-  fReferences->SetOwner();
-  TIter iter(histos);
-  Int_t arrayPos = 0, pbin = 0; 
-  AliPID::EParticleType species;
-  TObject *tmp = NULL;
-  TString histname;
-  TH1 *hnew = NULL;
-  while((tmp = iter())){
-    histname = tmp->GetName();
-    if(histname.BeginsWith("fHQel")){ // Electron histogram
-      histname.ReplaceAll("fHQel_p","");
-      species = AliPID::kElectron;
-    } else if(histname.BeginsWith("fHQmu")){ // Muon histogram
-      histname.ReplaceAll("fHQmu_p","");
-      species = AliPID::kMuon;
-    } else if(histname.BeginsWith("fHQpi")){ // Pion histogram
-      histname.ReplaceAll("fHQpi_p","");
-      species = AliPID::kPion;
-    }else if(histname.BeginsWith("fHQka")){ // Kaon histogram
-      histname.ReplaceAll("fHQka_p","");
-      species = AliPID::kKaon;
-    }else if(histname.BeginsWith("fHQpr")){ // Proton histogram
-      histname.ReplaceAll("fHQpr_p","");
-      species = AliPID::kProton;
-    } else continue;
-    pbin = histname.Atoi() - 1;
-    AliDebug(1, Form("Species %d, Histname %s, Pbin %d, Position in container %d", species, histname.Data(), pbin, arrayPos));
-    fMapRefHists[species][pbin] = arrayPos;
-    fReferences->AddAt((hnew = new TH1F(*dynamic_cast<TH1F *>(tmp))), arrayPos);
-    hnew->SetDirectory(gROOT);
-    arrayPos++;
-  }
-  if(!arrayPos){
-    AliError("Failed loading References");
-    return kFALSE;
-  }
-  AliDebug(2, Form("Successfully loaded %d Reference Histograms", arrayPos));
   SetBit(kIsOwner, kTRUE);
+  AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDReference->GetNumberOfMomentumBins()));
   return kTRUE;
 }
 
 //____________________________________________________________
-Bool_t AliTRDPIDResponse::GetResponse(Int_t n, Double_t *dedx, Float_t *p, Double_t prob[AliPID::kSPECIES], Bool_t kNorm) const
+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
 {
   //
   // Calculate TRD likelihood values for the track based on dedx and 
@@ -240,21 +134,21 @@ Bool_t AliTRDPIDResponse::GetResponse(Int_t n, Double_t *dedx, Float_t *p, Doubl
   //   true if calculation success
   // 
 
-  if(!fReferences){
+  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.);
+  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(!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]);
+      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];
     }
@@ -287,58 +181,27 @@ Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t pl
   // Interpolation between momentum bins
   //
   AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
-  Int_t pbin = GetLowerMomentumBin(plocal);  
-  AliDebug(1, Form("Bin %d", pbin));
-  Double_t pLayer = 0.;
+  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(pbin >= 0 && pbin < kNPBins-1){
-    TH1 *refHistUpper = NULL, *refHistLower = NULL;
-    if(fMapRefHists[species][pbin] >= 0)
-      refHistLower = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][pbin]));
-    if(fMapRefHists[species][pbin+1] >= 0)
-      refHistUpper = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][pbin+1]));
-    AliDebug(1, Form("Reference Histos (Upper/Lower): [%p|%p]", refHistUpper, refHistLower));
-  
-    if (refHistLower && refHistUpper ) {
-      Double_t pLower = refHistLower->GetBinContent(refHistLower->GetXaxis()->FindBin(dEdx));
-      Double_t pUpper = refHistUpper->GetBinContent(refHistUpper->GetXaxis()->FindBin(dEdx));
+  if(refLower && refUpper){
+    Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx));
+    Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx));
   
-      pLayer = pLower + (pUpper - pLower)/(fgkPBins[pbin+1]-fgkPBins[pbin]) * (plocal - fgkPBins[pbin]); 
-    }
-  }
-  else{
-    TH1 *refHist = NULL;
-    if(pbin < 0){
-      // underflow
-      if(fMapRefHists[species][0] >= 0){
-        refHist = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][0])); 
-      }
-    } else {
-      // overflow
-      if(fMapRefHists[species][kNPBins-1] > 0)
-        refHist = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][kNPBins-1])); 
-    }
-    if (refHist)
-      pLayer = refHist->GetBinContent(refHist->GetXaxis()->FindBin(dEdx));
-  }
-  AliDebug(1, Form("Probability %f", pLayer));
-  return pLayer;
-}
-
-//____________________________________________________________
-Int_t AliTRDPIDResponse::GetLowerMomentumBin(Double_t p) const {
-  //
-  // Get the momentum bin for a given momentum value
-  //
-  Int_t bin = -1;
-  if(p > fgkPBins[kNPBins-1]) return kNPBins-1;
-  for(Int_t ibin = 0; ibin < kNPBins - 1; ibin++){
-    if(p >= fgkPBins[ibin] && p < fgkPBins[ibin+1]){
-      bin = ibin;
-      break;
-    }
+    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");
   }
-  return bin;
+  AliDebug(1, Form("Probability %f", probLayer));
+  return probLayer;
 }
 
 //____________________________________________________________
@@ -346,18 +209,18 @@ void AliTRDPIDResponse::SetOwner(){
   //
   // Make Deep Copy of the Reference Histograms
   //
-  if(!fReferences || IsOwner()) return;
-  TObjArray *ctmp = new TObjArray();
-  for(Int_t ien = 0; ien < fReferences->GetEntriesFast(); ien++)
-    ctmp->AddAt(fReferences->UncheckedAt(ien)->Clone(), ien);
-  delete fReferences;
-  fReferences = ctmp;
+  if(!fkPIDReference || IsOwner()) return;
+  const AliTRDPIDReference *tmp = fkPIDReference;
+  fkPIDReference = dynamic_cast<const AliTRDPIDReference *>(tmp->Clone());
   SetBit(kIsOwner, kTRUE);
 }
 
 //____________________________________________________________
-Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, Double_t *in, Double_t *out) const
+Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out) const
 {
+       //
+       // Recalculate dE/dx
+       //
   switch(fPIDmethod){
   case kNN: // NN 
     break;
index 37ef8fd..31b0ff1 100644 (file)
@@ -30,6 +30,7 @@
 
 class TObjArray;
 class AliVTrack;
+class AliTRDPIDReference;
 class AliTRDPIDResponse : public TObject {
   public:
     enum ETRDPIDResponseStatus {
@@ -51,11 +52,11 @@ class AliTRDPIDResponse : public TObject {
     };
     AliTRDPIDResponse();
     AliTRDPIDResponse(const AliTRDPIDResponse &ref);
-    AliTRDPIDResponse& operator=(const AliTRDPIDResponse &);
+    AliTRDPIDResponse& operator=(const AliTRDPIDResponse &ref);
     ~AliTRDPIDResponse();
     
     ETRDPIDMethod     GetPIDmethod() const { return fPIDmethod;}
-    Bool_t    GetResponse(Int_t n, Double_t *dedx, Float_t *p, Double_t prob[AliPID::kSPECIES], Bool_t kNorm=kTRUE) const;
+    Bool_t    GetResponse(Int_t n, const Double_t * const dedx, const Float_t * const p, Double_t prob[AliPID::kSPECIES], Bool_t kNorm=kTRUE) const;
     inline ETRDNslices  GetNumberOfSlices() const;
     
     Bool_t    IsOwner() const {return TestBit(kIsOwner);}
@@ -64,18 +65,14 @@ class AliTRDPIDResponse : public TObject {
     void      SetPIDmethod(ETRDPIDMethod m) {fPIDmethod=m;}
     void      SetGainNormalisationFactor(Double_t gainFactor) { fGainNormalisationFactor = gainFactor; }
 
-    Bool_t    Load(const Char_t *filename = NULL);
-    Bool_t    Load(const TObjArray *histos);
+    Bool_t    Load(const Char_t *filename = NULL, const Char_t *refName = "RefTRDLQ1D");
+    Bool_t    Load(const AliTRDPIDReference *ref) { fkPIDReference = ref; return kTRUE; }
   
   private:
-    Bool_t    CookdEdx(Int_t nSlice, Double_t *in, Double_t *out) const;
-    Int_t     GetLowerMomentumBin(Double_t p) const;
+    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;
     
-    static const Double_t fgkPBins[kNPBins];
-    TObjArray *fReferences; // Container for reference distributions
-    Int_t     fMapRefHists[AliPID::kSPECIES][kNPBins];     
-                            // Map for the position of a given historgam in the container 
+    const AliTRDPIDReference *fkPIDReference;   // PID References
     Double_t  fGainNormalisationFactor;  // Gain normalisation factor
     ETRDPIDMethod   fPIDmethod;   // PID method selector  
   
index effa58e..33da86c 100644 (file)
@@ -1,4 +1,4 @@
-set ( SRCS  AliVParticle.cxx AliVTrack.cxx AliVCluster.cxx AliVCaloCells.cxx AliVVertex.cxx AliVEvent.cxx AliMixedEvent.cxx AliVHeader.cxx AliVEventHandler.cxx AliVEventPool.cxx AliVCuts.cxx AliVVZERO.cxx AliCentrality.cxx AliEventplane.cxx AliPID.cxx AliLog.cxx AliRunTag.cxx AliLHCTag.cxx AliDetectorTag.cxx AliEventTag.cxx AliFileTag.cxx AliEventTagCuts.cxx AliRunTagCuts.cxx AliLHCTagCuts.cxx AliDetectorTagCuts.cxx AliTagCreator.cxx AliHeader.cxx AliGenEventHeader.cxx AliDetectorEventHeader.cxx AliGenPythiaEventHeader.cxx AliGenCocktailEventHeader.cxx AliGenGeVSimEventHeader.cxx AliGenHijingEventHeader.cxx AliCollisionGeometry.cxx AliGenDPMjetEventHeader.cxx AliGenHerwigEventHeader.cxx AliGenEposEventHeader.cxx AliStack.cxx AliMCEventHandler.cxx AliInputEventHandler.cxx AliTrackReference.cxx AliSysInfo.cxx AliMCEvent.cxx AliMCParticle.cxx AliMCVertex.cxx AliMagF.cxx AliMagWrapCheb.cxx AliCheb3D.cxx AliCheb3DCalc.cxx AliNeutralTrackParam.cxx AliCodeTimer.cxx AliPDG.cxx AliTimeStamp.cxx AliTriggerScalers.cxx AliTriggerScalersRecord.cxx AliExternalTrackParam.cxx AliQA.cxx AliITSPidParams.cxx AliPIDResponse.cxx AliITSPIDResponse.cxx AliTPCPIDResponse.cxx AliTOFPIDResponse.cxx AliTRDPIDResponse.cxx AliDAQ.cxx AliRefArray.cxx)
+set ( SRCS  AliVParticle.cxx AliVTrack.cxx AliEventplane.cxx AliVCluster.cxx AliVCaloCells.cxx AliVVertex.cxx AliVEvent.cxx AliMixedEvent.cxx AliVHeader.cxx AliVEventHandler.cxx AliVEventPool.cxx AliVCuts.cxx AliVVZERO.cxx AliCentrality.cxx AliPID.cxx AliLog.cxx AliRunTag.cxx AliLHCTag.cxx AliDetectorTag.cxx AliEventTag.cxx AliFileTag.cxx AliEventTagCuts.cxx AliRunTagCuts.cxx AliLHCTagCuts.cxx AliDetectorTagCuts.cxx AliTagCreator.cxx AliHeader.cxx AliGenEventHeader.cxx AliDetectorEventHeader.cxx AliGenPythiaEventHeader.cxx AliGenCocktailEventHeader.cxx AliGenGeVSimEventHeader.cxx AliGenHijingEventHeader.cxx AliCollisionGeometry.cxx AliGenDPMjetEventHeader.cxx AliGenHerwigEventHeader.cxx AliGenEposEventHeader.cxx AliStack.cxx AliMCEventHandler.cxx AliInputEventHandler.cxx AliTrackReference.cxx AliSysInfo.cxx AliMCEvent.cxx AliMCParticle.cxx AliMCVertex.cxx AliMagF.cxx AliMagWrapCheb.cxx AliCheb3D.cxx AliCheb3DCalc.cxx AliNeutralTrackParam.cxx AliCodeTimer.cxx AliPDG.cxx AliTimeStamp.cxx AliTriggerScalers.cxx AliTriggerScalersRecord.cxx AliExternalTrackParam.cxx AliQA.cxx AliITSPidParams.cxx AliTRDPIDReference.cxx AliPIDResponse.cxx AliITSPIDResponse.cxx AliTPCPIDResponse.cxx AliTOFPIDResponse.cxx AliTRDPIDResponse.cxx AliDAQ.cxx AliRefArray.cxx)
 
 string(REPLACE ".cxx" ".h" HDRS  "${SRCS}")
 
index a4acecb..fcf2acf 100644 (file)
@@ -83,6 +83,7 @@
 #pragma link C++ class  AliExternalTrackParam+;
 #pragma link C++ class AliQA+;
 
+#pragma link C++ class AliTRDPIDReference+;
 #pragma link C++ class AliITSPidParams+;
 #pragma link C++ class AliPIDResponse+;
 #pragma link C++ class AliITSPIDResponse+;
index 1f50533..c558283 100644 (file)
@@ -35,6 +35,7 @@
 #include "AliCDBEntry.h"
 #include "AliLog.h"
 
+#include "AliTRDPIDReference.h"
 #include "AliTRDcalibDB.h"
 
 #include "Cal/AliTRDCalROC.h"
@@ -1163,7 +1164,19 @@ AliTRDPIDResponse *AliTRDcalibDB::GetPIDResponse(AliTRDPIDResponse::ETRDPIDMetho
     fPIDResponse = new AliTRDPIDResponse;
     // Load Reference Histos from OCDB
     fPIDResponse->SetPIDmethod(method);
-    fPIDResponse->Load(dynamic_cast<const TObjArray *>(GetCachedCDBObject(kIDPIDLQ1D)));
+    const TObjArray *references = dynamic_cast<const TObjArray *>(GetCachedCDBObject(kIDPIDLQ1D));
+    TIter refs(references);
+    TObject *obj = NULL;
+    AliTRDPIDReference *ref = NULL;
+    Bool_t hasReference = kFALSE;
+    while((obj = refs())){
+       if((ref = dynamic_cast<AliTRDPIDReference *>(obj))){
+               fPIDResponse->Load(ref);
+               hasReference = kTRUE;
+               break;
+       }
+    }
+    if(!hasReference) AliError("Reference histograms not found in the OCDB");
   }
   return fPIDResponse;
 }