Implementation of 1D PID LQ method

diff --git a/OCDB/TRD/Calib/PIDLQ/Run0_999999999_v0_s0.root b/OCDB/TRD/Calib/PIDLQ/Run0_999999999_v0_s0.root
index 4632f80..3bf6460 100644 (file)
Binary files a/OCDB/TRD/Calib/PIDLQ/Run0_999999999_v0_s0.root and b/OCDB/TRD/Calib/PIDLQ/Run0_999999999_v0_s0.root differ

diff --git a/TRD/AliTRDrecoParam.cxx b/TRD/AliTRDrecoParam.cxx
index dc97789..89c7d19 100644 (file)
--- a/TRD/AliTRDrecoParam.cxx
+++ b/TRD/AliTRDrecoParam.cxx
@@ -210,3 +210,29 @@ Float_t AliTRDrecoParam::GetNClusters() const
   nclusters *= 1.+fkNClusterNoise;
   return nclusters;
 }
+
+//______________________________________________________________
+void AliTRDrecoParam::SetPIDLQslices(Int_t s)
+{
+// Setting number of slices used by the PID LQ method s={1, 2}
+// If PID NN is set this function will change to PID LQ.
+ 
+  if(IsPIDNeuralNetwork()){
+    AliWarning("PID set to NN. Changing to LQ.");
+    SetPIDNeuralNetwork(kFALSE);
+  } 
+
+  switch(s){
+  case 1: 
+    if(TESTBIT(fFlags, kLQ2D)) CLRBIT(fFlags, kLQ2D);
+    break;
+  case 2:
+    SETBIT(fFlags, kLQ2D);
+    break;
+  default:
+    AliWarning(Form("N[%d] PID LQ slices not implemented. Using default 2.", s));
+    SETBIT(fFlags, kLQ2D);
+    break;
+  }
+}
+

diff --git a/TRD/AliTRDrecoParam.h b/TRD/AliTRDrecoParam.h
index 7a8bb14..4667bdd 100644 (file)
--- a/TRD/AliTRDrecoParam.h
+++ b/TRD/AliTRDrecoParam.h
@@ -30,16 +30,17 @@ public:
     kTRDreconstructionTasks = 3
   };
   enum ETRDflags {
-    kDriftGas,
-    kVertexConstraint,
-    kTailCancelation,
-    kImproveTracklet, 
-    kLUT, 
-    kGAUS,
-    kClusterSharing,
-    kSteerPID,
-    kEightSlices,
-    kCheckTimeConsistency
+    kDriftGas
+    ,kVertexConstraint
+    ,kTailCancelation
+    ,kImproveTracklet
+    ,kLUT
+    ,kGAUS
+    ,kClusterSharing
+    ,kSteerPID
+    ,kEightSlices
+    ,kCheckTimeConsistency
+    ,kLQ2D
   };
   AliTRDrecoParam();
   AliTRDrecoParam(const AliTRDrecoParam &rec);
@@ -55,6 +56,7 @@ public:
   Double_t GetNMeanClusters() const         { return fkNMeanClusters; }
   Double_t GetNSigmaClusters() const        { return fkNSigmaClusters; }
   Double_t GetFindableClusters() const      { return fkFindable; }
+  inline Int_t    GetPIDLQslices() const;
   Double_t GetMaxTheta() const              { return fkMaxTheta; }
   Double_t GetMaxPhi() const                { return fkMaxPhi;   }
   Double_t GetPlaneQualityThreshold() const { return fkPlaneQualityThreshold; }
@@ -103,6 +105,7 @@ public:
   void     SetLUT(Bool_t b=kTRUE)                             {if(b) SETBIT(fFlags, kLUT); else CLRBIT(fFlags, kLUT);}
   void     SetGAUS(Bool_t b=kTRUE)                            {if(b) SETBIT(fFlags, kGAUS); else CLRBIT(fFlags, kGAUS);}
   void     SetPIDNeuralNetwork(Bool_t b=kTRUE)                {if(b) SETBIT(fFlags, kSteerPID); else CLRBIT(fFlags, kSteerPID);}
+  void     SetPIDLQslices(Int_t s);
   void     SetTailCancelation(Bool_t b=kTRUE)                 {if(b) SETBIT(fFlags, kTailCancelation); else CLRBIT(fFlags, kTailCancelation);}
   void     SetXenon(Bool_t b = kTRUE)                         {if(b) CLRBIT(fFlags, kDriftGas); else SETBIT(fFlags, kDriftGas);}
   void     SetVertexConstrained()                             {SETBIT(fFlags, kVertexConstraint);}
@@ -236,4 +239,12 @@ inline void AliTRDrecoParam::SetTCParams(Double_t *par)
   if(!par) return;
   memcpy(fTCParams, par, 8*sizeof(Double_t));
 }
+
+//___________________________________________________
+inline Int_t AliTRDrecoParam::GetPIDLQslices() const
+{
+  if(IsPIDNeuralNetwork()) return -1;
+  return TESTBIT(fFlags, kLQ2D) ? 2 : 1;
+}
+
 #endif

diff --git a/TRD/AliTRDseedV1.cxx b/TRD/AliTRDseedV1.cxx
index ad35037..5147d38 100644 (file)
--- a/TRD/AliTRDseedV1.cxx
+++ b/TRD/AliTRDseedV1.cxx
@@ -585,8 +585,9 @@ Bool_t AliTRDseedV1::CookPID()
   AliDebug(4, Form("p=%6.4f[GeV/c] dEdx{%7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f} l=%4.2f[cm]", GetMomentum(), fdEdx[0], fdEdx[1], fdEdx[2], fdEdx[3], fdEdx[4], fdEdx[5], fdEdx[6], fdEdx[7], length));
 
   // Sets the a priori probabilities
+  Bool_t kPIDNN(fkReconstructor->GetPIDMethod()==AliTRDpidUtil::kNN);
   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++)
-    fProb[ispec] = pd->GetProbability(ispec, GetMomentum(), &fdEdx[0], length, GetPlane());
+    fProb[ispec] = pd->GetProbability(ispec, GetMomentum(), &fdEdx[0], length, kPIDNN?GetPlane():fkReconstructor->GetRecoParam()->GetPIDLQslices());
   
   return kTRUE;
 }

diff --git a/TRD/Cal/AliTRDCalPID.h b/TRD/Cal/AliTRDCalPID.h
index dee363f..e2b73e4 100644 (file)
--- a/TRD/Cal/AliTRDCalPID.h
+++ b/TRD/Cal/AliTRDCalPID.h
@@ -21,11 +21,9 @@
 
 class AliTRDCalPID : public TNamed
 {
-
- public:
-
+public:
   enum {
-    kNMom   = 11,
+    kNMom      = 11,
     kNSlicesLQ = 2,
     kNSlicesNN = 8
   };
@@ -49,10 +47,10 @@ class AliTRDCalPID : public TNamed
   static  const Char_t  *GetPartName(Int_t i)               { return fPartName[i];   }
   static  const Char_t  *GetPartSymb(Int_t i)               { return fPartSymb[i];   }
 
-          void     SetPartName(Int_t i, const Char_t *name) { fPartName[i] = name;   }
-          void     SetPartSymb(Int_t i, const Char_t *symb) { fPartSymb[i] = symb;   }
+  void     SetPartName(Int_t i, const Char_t *name) { fPartName[i] = name;   }
+  void     SetPartSymb(Int_t i, const Char_t *symb) { fPartSymb[i] = symb;   }
 
- protected:
+protected:
   virtual void     Init() = 0;
 
   static const Char_t   *fPartName[AliPID::kSPECIES];     //! Names of particle species

diff --git a/TRD/Cal/AliTRDCalPIDLQ.cxx b/TRD/Cal/AliTRDCalPIDLQ.cxx
index 80d28c2..d6a543a 100644 (file)
--- a/TRD/Cal/AliTRDCalPIDLQ.cxx
+++ b/TRD/Cal/AliTRDCalPIDLQ.cxx
@@ -34,6 +34,7 @@
 #include "AliPID.h"
 
 #include "TKDPDF.h"
+
 #include "AliTRDCalPIDLQ.h"
 #include "AliTRDcalibDB.h"
 
@@ -65,7 +66,6 @@ AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t *name, const Text_t *title)
 
 }
 
-
 //_________________________________________________________________________
 Bool_t AliTRDCalPIDLQ::LoadReferences(Char_t *refFile)
 {
@@ -82,7 +82,7 @@ Bool_t AliTRDCalPIDLQ::LoadReferences(Char_t *refFile)
     return kFALSE;
   }
   TObjArray *pdf(NULL);
-  if (!( pdf = dynamic_cast<TObjArray*>(gFile->Get("PDF_2DLQ")))) {
+  if (!( pdf = dynamic_cast<TObjArray*>(gFile->Get("PDF_LQ")))) {
     AliError("PID data not available");
     return kFALSE;
   }
@@ -92,34 +92,44 @@ Bool_t AliTRDCalPIDLQ::LoadReferences(Char_t *refFile)
 }
 
 //_________________________________________________________________________
-TObject* AliTRDCalPIDLQ::GetModel(Int_t ip, Int_t iType, Int_t iplane) const          // iType not needed
+TObject* AliTRDCalPIDLQ::GetModel(Int_t ip, Int_t is, Int_t slices) const
 {
   //
   // Returns one selected dEdx histogram
   //
 
-  if (iType < 0 || iType >= AliPID::kSPECIES) return NULL;
+  if (is < 0 || is >= AliPID::kSPECIES) return NULL;
   if(ip<0 || ip>= kNMom ) return NULL;
 
-  AliDebug(2, Form("Retrive dEdx distribution for %s @ p=%5.2f [GeV/c].", AliPID::ParticleShortName(iType), fgTrackMomentum[ip]));
-  return fModel->At(GetModelID(ip, iType, iplane));
+  AliDebug(2, Form("Retrive dEdx distribution for %s @ p=%5.2f [GeV/c].", AliPID::ParticleShortName(is), fgTrackMomentum[ip]));
+  return fModel->At(GetModelID(ip, is, slices));
+}
+
+//_________________________________________________________________________
+Int_t AliTRDCalPIDLQ::GetNrefs()
+{
+// Returns the number of PDF distribution 
+  return AliTRDCalPID::kNMom*AliPID::kSPECIES*2;
 }
 
 //_________________________________________________________________________
 Double_t AliTRDCalPIDLQ::GetProbability(Int_t spec, Float_t mom
                                       , const Float_t * const dedx
-                                      , Float_t length, Int_t iplane) const
+                                      , Float_t length
+                                      , Int_t slices) const
 {
-  //
-  // Core function of AliTRDCalPID class for calculating the
-  // likelihood for species "spec" (see AliTRDtrack::kNspecie) of a
-  // given momentum "mom" and a given dE/dx (slice "dedx") yield per
-  // layer
-  //
+//
+// Core function of AliTRDCalPID class for calculating the
+// likelihood for species "spec" (see AliTRDtrack::kNspecie) of a
+// given momentum "mom" and a given dE/dx (slice "dedx") yield per
+// layer
+//
 
   if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
-  Double_t x[AliTRDCalPID::kNSlicesLQ];
-  if(!CookdEdx(dedx, x)) return 0.;
+
+  Bool_t k2D(slices==2);
+  Double_t x[]={0., 0.};
+  if(!CookdEdx(dedx, x, k2D)) return 0.;
     
   // find the interval in momentum and track segment length which applies for this data
 /*  Int_t ilength = 1;
@@ -132,10 +142,9 @@ Double_t AliTRDCalPIDLQ::GetProbability(Int_t spec, Float_t mom
   Double_t p[2], e[2], r;
   TKDPDF *pdf(NULL);
 
-  AliDebug(2, Form("Looking for %s@p=%6.4f[GeV/c] dEdx={%7.2f %7.2f}[a.u.] l=%4.2f[cm].", AliPID::ParticleShortName(spec), mom, dedx[0], dedx[1], length));
-  if(!(pdf = dynamic_cast<TKDPDF*>(fModel->At(GetModelID(imom-1, spec, iplane))))) {
-    AliError(Form("Ref data @ idx[%d] not found in DB.", GetModelID(imom-1, spec, iplane)));
-    fModel->ls();
+  AliDebug(2, Form("Looking %cD for %s@p=%6.4f[GeV/c] log(dEdx)={%7.2f %7.2f}[a.u.] l=%4.2f[cm].", k2D?'2':'1', AliPID::ParticleShortName(spec), mom, x[0], x[1], length));
+  if(!(pdf = dynamic_cast<TKDPDF*>(fModel->At(GetModelID(imom-1, spec, slices))))) {
+    AliError(Form("%cD Ref data @ idx[%d] not found in DB.", k2D?'2':'1', GetModelID(imom-1, spec, slices)));
     return 0.;
   }
   if(!pdf->GetSize()) pdf->Bootstrap();
@@ -143,9 +152,8 @@ Double_t AliTRDCalPIDLQ::GetProbability(Int_t spec, Float_t mom
   p[0]=TMath::Abs(r); // conversion from interpolation to PDF
   AliDebug(2, Form("LQ=%6.3f+-%5.2f%% @ %4.1f[GeV/c]", p[0], 1.E2*e[0]/p[0], fgTrackMomentum[imom-1]));
 
-  if(!(pdf = dynamic_cast<TKDPDF*>(fModel->At(GetModelID(imom, spec, iplane))))){
-    AliError(Form("Ref data @ idx[%d] not found in DB.", GetModelID(imom, spec, iplane)));
-    fModel->ls();
+  if(!(pdf = dynamic_cast<TKDPDF*>(fModel->At(GetModelID(imom, spec, slices))))){
+    AliError(Form("%cD Ref data @ idx[%d] not found in DB.", k2D?'2':'1', GetModelID(imom, spec, slices)));
     return p[0];
   }
   if(!pdf->GetSize()) pdf->Bootstrap();
@@ -168,6 +176,6 @@ void AliTRDCalPIDLQ::Init()
 //
 // PID LQ list initialization
 //
-  fModel = new TObjArray(AliPID::kSPECIES  * kNMom);
-  fModel -> SetOwner();
+  fModel = new TObjArray(AliPID::kSPECIES  * kNMom * 2);
+  fModel->SetOwner();
 }

diff --git a/TRD/Cal/AliTRDCalPIDLQ.h b/TRD/Cal/AliTRDCalPIDLQ.h
index c2b9424..9fbd898 100644 (file)
--- a/TRD/Cal/AliTRDCalPIDLQ.h
+++ b/TRD/Cal/AliTRDCalPIDLQ.h
@@ -27,22 +27,21 @@ class AliTRDCalPIDLQ : public AliTRDCalPID
 public:
   enum ETRDCalPIDLQ {
     kNLength = 4 // No of bins for tracklet length discretization 
-   ,kNN2LQtransition = 4 // index of NN slices where first LQ slice ends 
+   ,kNN2LQtransition = 4 // index of NN slices where first LQ slice ends
   };
-       
   AliTRDCalPIDLQ();
   AliTRDCalPIDLQ(const Text_t *name, const Text_t *title);
   virtual        ~AliTRDCalPIDLQ(){};
 
-  inline static Bool_t  CookdEdx(const Float_t * const dedx, Double_t *x);
-  TObject*              GetModel(Int_t ip, Int_t iType, Int_t iPlane) const;
+  inline static Bool_t  CookdEdx(const Float_t * const dedx, Double_t *x, Bool_t k2D=kTRUE);
+  TObject*              GetModel(Int_t ip, Int_t iType, Int_t slices) const;
   static Double_t       GetLength(Int_t il) { return (il<0 || il>=kNLength) ? -1. : fgTrackSegLength[il]; }
-  inline static Int_t   GetModelID(Int_t mom, Int_t spec, Int_t ii=-1);
+  inline static Int_t   GetModelID(Int_t mom, Int_t spec, Int_t slices=1);
+  static Int_t          GetNrefs();
   Double_t              GetProbability(Int_t spec, Float_t mom
                                , const Float_t * const dedx
-                               , Float_t length, Int_t plane) const;
+                               , Float_t length, Int_t slices) const;
   Bool_t                LoadReferences(Char_t* refFile);
-
 protected:
   static Float_t  fgTrackSegLength[kNLength]; // Track segment lengths
 
@@ -55,33 +54,36 @@ private:
 };
 
 //_________________________________________________________________________
-inline Int_t AliTRDCalPIDLQ::GetModelID(Int_t mom, Int_t spec, Int_t /*ii*/)
-{  
-// Returns the ID of the PDF distribution 
-// 5 species * 11 momentum ranges
-
-  return spec * AliTRDCalPID::kNMom + mom;
-}
-
-//_________________________________________________________________________
-inline Bool_t  AliTRDCalPIDLQ::CookdEdx(const Float_t * const dedx, Double_t *x)
+inline Bool_t  AliTRDCalPIDLQ::CookdEdx(const Float_t * const dedx, Double_t *x, Bool_t k2D)
 {
 // Convert NN dEdx slices to the representation used by LQ
 
   x[0]=0.;x[1]=0.;
   for(Int_t islice=AliTRDCalPID::kNSlicesNN; islice--;){
-    Int_t jslice = islice>kNN2LQtransition;
+    Int_t jslice(0);
+    if(k2D) jslice = islice>kNN2LQtransition;
     x[jslice]+=dedx[islice];
   }
   
   // check data integrity
   if(x[0]<1.e-30) return kFALSE;
-  if(x[1]<1.e-30) return kFALSE;
+  if(k2D && x[1]<1.e-30) return kFALSE;
 
   x[0] = TMath::Log(x[0]);
-  x[1] = TMath::Log(x[1]);
+  x[1] = k2D ? TMath::Log(x[1]) : 0.;
   return kTRUE;
 }
 
+//_________________________________________________________________________
+inline Int_t AliTRDCalPIDLQ::GetModelID(Int_t mom, Int_t spec, Int_t ii)
+{
+// Returns the ID of the PDF distribution 
+// 5 species * 11 momentum ranges * 2 slices
+
+  Int_t n(GetNrefs()>>1);
+  return (ii-1)*n + (spec * AliTRDCalPID::kNMom + mom);
+}
+
+
 #endif