New 2-dim PID methods by Alexandru
authorcblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 8 May 2007 14:16:16 +0000 (14:16 +0000)
committercblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 8 May 2007 14:16:16 +0000 (14:16 +0000)
12 files changed:
TRD/AliTRDgtuTrack.cxx
TRD/AliTRDpidESD.cxx
TRD/AliTRDpidESD.h
TRD/AliTRDtrack.h
TRD/AliTRDtracker.cxx
TRD/Cal/AliTRDCalPIDLQ.cxx
TRD/Cal/AliTRDCalPIDLQ.h
TRD/Cal/AliTRDCalPIDLQRef.cxx [new file with mode: 0644]
TRD/Cal/AliTRDCalPIDLQRef.h [new file with mode: 0644]
TRD/Calib/PIDLQ/Run0_0_v0_s4.root [new file with mode: 0644]
TRD/TRDbaseLinkDef.h
TRD/libTRDbase.pkg

index 47befff..dad3d16 100644 (file)
@@ -448,8 +448,14 @@ void AliTRDgtuTrack::MakePID()
 
     // HEAD28Mar06 new distributions: factor = 6.0
 
-    probEle *= pd->GetProbability(0,TMath::Abs(fPt),q*6.0);
-    probPio *= pd->GetProbability(2,TMath::Abs(fPt),q*6.0);
+               // A.Bercuci on 2nd May 2007
+               // Temporary modification to deal with more energy slices. The values
+               // attached to charge slices and track length are dummy
+               Double_t dedx[3];
+               dedx[0] = dedx[1] = q*3.; dedx[2] = 0.;
+               Double_t length = 3.7;
+    probEle *= pd->GetProbability(0, TMath::Abs(fPt), dedx, length);
+    probPio *= pd->GetProbability(2, TMath::Abs(fPt), dedx, length);
 
   }
 
index f5a9f55..3de702c 100644 (file)
@@ -13,6 +13,8 @@
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
+/* $Id$ */
+
 ////////////////////////////////////////////////////////////////////////////
 //                                                                        //
 // Implementation of the TRD PID class                                    //
@@ -23,8 +25,9 @@
 // momentum (mom) and particle type k                                     //
 // from the precalculated distributions.                                  //
 //                                                                        //
-// Original version:                                                      //
-// Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de>                   //
+// Authors :                                                              //
+// Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de> (Original version)//
+// Alex Bercuci (a.bercuci@gsi.de)                                        //
 //                                                                        //
 ////////////////////////////////////////////////////////////////////////////
 
@@ -35,6 +38,8 @@
 #include "AliTRDpidESD.h"
 #include "AliTRDgeometry.h"
 #include "AliTRDcalibDB.h"
+#include "AliRun.h"
+#include "AliTRDtrack.h"
 #include "Cal/AliTRDCalPIDLQ.h"
 
 ClassImp(AliTRDpidESD)
@@ -44,7 +49,8 @@ ClassImp(AliTRDpidESD)
   Int_t AliTRDpidESD::fMinPlane          = 0;
 
 //_____________________________________________________________________________
-AliTRDpidESD::AliTRDpidESD():TObject()
+AliTRDpidESD::AliTRDpidESD()
+  :TObject()
 {
   //
   // Default constructor
@@ -53,7 +59,8 @@ AliTRDpidESD::AliTRDpidESD():TObject()
 }
 
 //_____________________________________________________________________________
-AliTRDpidESD::AliTRDpidESD(const AliTRDpidESD &p):TObject(p)
+AliTRDpidESD::AliTRDpidESD(const AliTRDpidESD &p)
+  :TObject(p)
 {
   //
   // AliTRDpidESD copy constructor
@@ -94,108 +101,168 @@ Int_t AliTRDpidESD::MakePID(AliESD *event)
   //
   // This function calculates the PID probabilities based on TRD signals
   //
-  // So far this method produces probabilities based on the total charge
-  // in each layer and the position of the maximum time bin in each layer.
-  // In a final version this should also exploit the charge measurement in
-  // the different slices of a given layer.
-  //  
-
-  Double_t p[10];
-  Int_t    nSpecies  = AliPID::kSPECIES;
-  Int_t    nPlanePID = 0;
-  Double_t mom       = 0.0;
-  Double_t probTotal = 0.0;
-
-  AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
-  if (!calibration) {
-    AliErrorGeneral("AliTRDpidESD::MakePID"
-                   ,"No access to calibration data\n");
-    return -1;
-  }  
-
-  // Retrieve the CDB container class with the probability distributions
-  const AliTRDCalPIDLQ *pd = calibration->GetPIDLQObject();
-  if (!pd) {
-    AliErrorGeneral("AliTRDpidESD::MakePID"
-                   ,"No access to AliTRDCalPIDLQ\n");
-    return -1;
-  }
-
-  // Loop through all ESD tracks
-  Int_t ntrk = event->GetNumberOfTracks();
-  for (Int_t i = 0; i < ntrk; i++) {
-
-    AliESDtrack *t = event->GetTrack(i);
-
-    // Check the ESD track status
-    if (fCheckTrackStatus) {
-      if (((t->GetStatus() & AliESDtrack::kTRDout  ) == 0) &&
-         ((t->GetStatus() & AliESDtrack::kTRDrefit) == 0)) {
-        continue;
-      }
-    }
-
-    // Check for ESD kink tracks
-    if (fCheckKinkStatus) {
-      if (t->GetKinkIndex(0) != 0) {
-        continue;
-      }
-    }
-
-    // Skip tracks that have no TRD signal at all
-    if (t->GetTRDsignal() == 0) {
-      continue;
-    }
-
-    mom       = t->GetP();
-    probTotal = 0.0;
-    nPlanePID = 0;
-    for (Int_t iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
-      p[iSpecies] = 1.0;
-    }
-
-    // Check the different detector layers
-    for (Int_t iPlan = 0; iPlan < AliTRDgeometry::kNplan; iPlan++) {
-
-      // Use the total charge in a given plane
-      Double_t dedx    = t->GetTRDsignals(iPlan,-1);
-      Int_t    timebin = t->GetTRDTimBin(iPlan);
-      if ((dedx    >  0.0) && 
-          (timebin > -1.0)) {
-
-        nPlanePID++;
-
-        // Get the probabilities for the different particle species
-        for (Int_t iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
-
-         p[iSpecies] *= pd->GetProbability(iSpecies,mom,dedx);
-         p[iSpecies] *= pd->GetProbabilityT(iSpecies,mom,timebin);
-         p[iSpecies] *= 100.0; // ??????????????
-
+  // The method produces probabilities based on the charge
+  // and the position of the maximum time bin in each layer.
+  // The dE/dx information can be used as global charge or 2 to 3
+  // slices. Check AliTRDCalPIDLQ and AliTRDCalPIDLQRef for the actual
+  // implementation.
+  //
+  // Author
+  // Alex Bercuci (A.Bercuci@gsi.de) 2nd May 2007
+
+       AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
+       if (!calibration) {
+               AliErrorGeneral("AliTRDpidESD::MakePID()"
+                       ,"No access to calibration data");
+               return -1;
+       }
+       
+       // Retrieve the CDB container class with the probability distributions
+       const AliTRDCalPIDLQ *pd = calibration->GetPIDLQObject();
+       if (!pd) {
+               AliErrorGeneral("AliTRDpidESD::MakePID()"
+                       ,"No access to AliTRDCalPIDLQ");
+               return -1;
        }
 
-      }
-
-    }
-
-    for (Int_t iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
-      probTotal += p[iSpecies];
-    }
 
-    for (Int_t iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
-      if ((probTotal >       0.0) &&
-          (nPlanePID > fMinPlane)) {
-        p[iSpecies] /= probTotal;
-      }
-      else {
-        p[iSpecies] = 1.0;
-      }
-    }
+       // Loop through all ESD tracks
+       Double_t p[10];
+       AliESDtrack *t = 0x0;
+       Double_t dedx[AliTRDtrack::kNslice], dEdx;
+       Int_t    timebin;
+       Float_t mom, length, probTotal;
+       Int_t nPlanePID;
+       for (Int_t i=0; i<event->GetNumberOfTracks(); i++) {
+               t = event->GetTrack(i);
+
+               // Check track
+               if(!CheckTrack(t)) continue;
+                                               
+               // Skip tracks which have no TRD signal at all
+               if (t->GetTRDsignal() == 0.) continue;
+       
+               // Loop over detector layers
+               mom       = 0.; //t->GetP();
+               length    = 0.;
+               probTotal = 0.;
+               nPlanePID    = 0;
+               for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) p[iSpecies] = 1.;
+               for (Int_t iPlan = 0; iPlan < AliTRDgeometry::kNplan; iPlan++) {
+                       // read data for track segment
+                       for(int iSlice=0; iSlice<AliTRDtrack::kNslice; iSlice++)
+                               dedx[iSlice] = t->GetTRDsignals(iPlan, iSlice);
+                       dEdx    = t->GetTRDsignals(iPlan, -1);
+                       timebin = t->GetTRDTimBin(iPlan);
+
+                       // check data
+                       if ((dEdx <=  0.) || (timebin <= -1.)) continue;
+
+                       // retrive kinematic info for this track segment
+                       if(!GetTrackSegmentKine(t, iPlan, mom, length)) continue;
+                       
+                       // this track segment has fulfilled all requierments
+                       nPlanePID++;
+                       
+                       // Get the probabilities for the different particle species
+                       for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
+                               p[iSpecies] *= pd->GetProbability(iSpecies, mom, dedx, length);
+                               p[iSpecies] *= pd->GetProbabilityT(iSpecies, mom, timebin);
+                               probTotal   += p[iSpecies];
+                       }
+               }
+
+               // normalize probabilities
+               if(probTotal > 0.)
+                       for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
+                               if(nPlanePID > fMinPlane) p[iSpecies] /= probTotal;
+                               else p[iSpecies] = 1.0;
+
+
+               // book PID to the track
+               t->SetTRDpid(p);
+       }
+       
+       return 0;
+}
 
-    t->SetTRDpid(p);
+//_____________________________________________________________________________
+Bool_t AliTRDpidESD::CheckTrack(AliESDtrack *t)
+{
+  //
+  // Check if track is eligible for PID calculations
+  //
+       
+       // Check the ESD track status
+       if (fCheckTrackStatus) {
+               if (((t->GetStatus() & AliESDtrack::kTRDout  ) == 0) &&
+                       ((t->GetStatus() & AliESDtrack::kTRDrefit) == 0)) return kFALSE;
+       }
 
-  }
+       // Check for ESD kink tracks
+       if (fCheckKinkStatus && (t->GetKinkIndex(0) != 0)) return kFALSE;
 
-  return 0;
+       return kTRUE;
+}
 
+//_____________________________________________________________________________
+Bool_t AliTRDpidESD::GetTrackSegmentKine(AliESDtrack *t, Int_t plan, Float_t &mom, Float_t &length)
+{
+  //
+  // Retrive momentum "mom" and track "length" in TRD chamber from plane
+  // "plan" according to information stored in AliESDtrack "t".
+  // 
+
+       if(!gAlice){
+               AliErrorGeneral("AliTRDpidESD::GetTrackSegmentKine()"
+               ,"No gAlice object to retrive TRDgeometry and Magnetic fied  - this has to be removed in the future.");
+               return kFALSE;
+       }
+       
+       // Retrieve TRD geometry -> Maybe there is a better way to do this
+       Bool_t kSelfGeom = kFALSE;
+       AliTRDgeometry *TRDgeom =0x0;
+       if(gAlice) TRDgeom = AliTRDgeometry::GetGeometry(gAlice->GetRunLoader());
+       if(!TRDgeom){
+               AliWarningGeneral("AliTRDpidESD::GetTrackSegmentKine()", "Cannot load TRD geometry from gAlice! Build a new one.\n");
+               TRDgeom = new AliTRDgeometry();
+               kSelfGeom = kTRUE;
+       }
+       const Float_t kAmHalfWidth = AliTRDgeometry::AmThick() / 2.;
+  const Float_t kDrWidth = AliTRDgeometry::DrThick();
+       
+
+       // retrive the magnetic field
+       Double_t xyz0[3] = { 0., 0., 0.}, xyz1[3];
+       Double_t b[3], alpha;
+       gAlice->Field(xyz0,b);      // b[] is in kilo Gauss
+       Float_t field = b[2] * 0.1; // Tesla
+
+               
+       // find momentum at chamber entrance and track length in chamber
+       AliExternalTrackParam *param = (plan<3) ? new AliExternalTrackParam(*t->GetInnerParam()) : new AliExternalTrackParam(*t->GetOuterParam());
+
+       param->PropagateTo(TRDgeom->GetTime0(plan)+kAmHalfWidth, field);
+       param->GetXYZ(xyz0);
+       alpha = param->GetAlpha();
+       param->PropagateTo(TRDgeom->GetTime0(plan)-kAmHalfWidth-kDrWidth, field);
+       // eliminate track segments which are crossing SM boundaries along chamber
+       if(fabs(alpha-param->GetAlpha())>.01){
+               delete param;
+               if(kSelfGeom) delete TRDgeom;
+               return kFALSE;
+       }
+       param->GetXYZ(xyz1);
+       length = sqrt(
+               (xyz1[0]-xyz0[0])*(xyz1[0]-xyz0[0])+
+               (xyz1[1]-xyz0[1])*(xyz1[1]-xyz0[1])+
+               (xyz1[2]-xyz0[2])*(xyz1[2]-xyz0[2])
+       );
+       param->GetPxPyPz(xyz1);
+       mom = sqrt(xyz1[0]*xyz1[0] + xyz1[1]*xyz1[1] + xyz1[2]*xyz1[2]);
+       delete param;
+       if(kSelfGeom) delete TRDgeom;
+
+       return kTRUE;
 }
+
index 974c6e5..0f29454 100644 (file)
@@ -16,7 +16,7 @@
 #include <TObject.h>
 
 class AliESD;
-
+class AliESDtrack;
 class AliTRDpidESD : public TObject {
 
  public:
@@ -27,7 +27,7 @@ class AliTRDpidESD : public TObject {
   AliTRDpidESD &operator=(const AliTRDpidESD &p);
 
   virtual void    Copy(TObject &p) const;
-
+  static  Bool_t  CheckTrack(AliESDtrack *t);
   static  Int_t   MakePID(AliESD *event);
 
           void    SetCheckTrackStatus(Bool_t status = kTRUE) { fCheckTrackStatus = status; };
@@ -37,7 +37,7 @@ class AliTRDpidESD : public TObject {
          Bool_t  GetCheckTrackStatus()                      { return fCheckTrackStatus;   };      
          Bool_t  GetCheckKinkStatus()                       { return fCheckKinkStatus;    };      
           Int_t   GetMinPlane()                              { return fMinPlane;           };
-
+  static  Bool_t  GetTrackSegmentKine(AliESDtrack *t, Int_t plan, Float_t &mom, Float_t &length);
  private:
 
   static  Bool_t  fCheckTrackStatus;    // Enable check on ESD track status
index 7bf255d..058b96d 100644 (file)
@@ -18,6 +18,7 @@
 
 class AliESDtrack;
 class AliTrackReference;
+class AliTRDcluster;
 
 const unsigned kMAXCLUSTERSPERTRACK = 210; 
 
@@ -30,7 +31,8 @@ class AliTRDtrack : public AliKalmanTrack {
        , kNplane    =   6
        , kNcham     =   5
        , kNsect     =  18
-       , kNslice    =   3 };
+       , kNslice    =   3
+       , kMidTimeBin=  14};
 
    AliTRDtrack();
    AliTRDtrack(const AliTRDcluster *c, Int_t index, const Double_t xx[5]
index 50b3774..9efdabc 100644 (file)
@@ -3468,57 +3468,65 @@ Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster *c)
 void AliTRDtracker::CookdEdxTimBin(AliTRDtrack &TRDtrack)
 {
   //
-  // This is setting fdEdxPlane and fTimBinPlane
-  // Sums up the charge in each plane for track TRDtrack and also get the 
-  // Time bin for Max. Cluster
-  // Prashant Shukla (shukla@physi.uni-heidelberg.de)
+  // Build the information needed for TRD PID calculations. This are:
+  // - dE/dx - total
+  // - dE/dx - slices (0-14, 15-22, NONE). The Edep per slice is
+  //    calculated as the sum of clusters charge weighted with average
+  //    value from the <PH> spectrum
+  // - Time bin for Max. Cluster
   //
+  // Authors:
+  // Prashant Shukla (shukla@physi.uni-heidelberg.de)
+  // Alex Bercuci (a.bercuci@gsi.de)
 
   Double_t  clscharge[AliESDtrack::kNPlane][AliESDtrack::kNSlice];
   Double_t  maxclscharge[AliESDtrack::kNPlane];
   Int_t     nCluster[AliESDtrack::kNPlane][AliESDtrack::kNSlice];
   Int_t     timebin[AliESDtrack::kNPlane];
 
-  // Initialization of cluster charge per plane.  
+  // Initialization of variables  
   for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
     for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
-      clscharge[iPlane][iSlice] = 0.0;
+      clscharge[iPlane][iSlice] = 0.;
       nCluster[iPlane][iSlice]  = 0;
     }
-  }
-
-  // Initialization of cluster charge per plane.  
-  for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
     timebin[iPlane]      =  -1;
     maxclscharge[iPlane] = 0.0;
   }
 
+       // Average PH spectrum used for weighting the edep. It has to
+       // come eather from caligration DB or as parameter of the class. To
+       // be discuss with software responsibles
+       const Float_t fAveragePH[] = { 5.010e-03, 2.017e-02, 6.221e-02, 6.706e-02
+                                     , 5.551e-02, 4.812e-02, 4.395e-02, 4.219e-02
+                                     , 4.172e-02, 4.156e-02, 4.162e-02, 4.204e-02
+                                     , 4.278e-02, 4.367e-02, 4.460e-02, 4.554e-02
+                                     , 4.674e-02, 4.810e-02, 5.005e-02, 5.258e-02
+                                     , 5.587e-02, 5.892e-02, 5.587e-02, 5.258e-02 };
+       
   // Loop through all clusters associated to track TRDtrack
-  Int_t nClus = TRDtrack.GetNumberOfClusters();  // from Kalmantrack
-  for (Int_t iClus = 0; iClus < nClus; iClus++) {
+  AliTRDcluster *pTRDcluster = 0x0;
+  for (Int_t iClus = 0; iClus < TRDtrack.GetNumberOfClusters(); iClus++) {
     Double_t charge = TRDtrack.GetClusterdQdl(iClus);
     Int_t    index  = TRDtrack.GetClusterIndex(iClus);
-    AliTRDcluster *pTRDcluster = (AliTRDcluster *) GetCluster(index); 
-    if (!pTRDcluster) {
-      continue;
-    }
-    Int_t    tb     = pTRDcluster->GetLocalTimeBin();
-    if (!tb) {
-      continue;
-    }
+    if (!(pTRDcluster = (AliTRDcluster*)GetCluster(index))) continue;
+    
+               // check negative time bins ?!?!?
+               Int_t    tb     = pTRDcluster->GetLocalTimeBin();
+
+               // temporary until fix in GetT0
+               if(tb<0){
+                       Warning("CookdEdxTimBin()", Form("Negative time bin value in track %d cluster %d", TRDtrack.GetLabel(), iClus));
+                       continue;
+               }
+               
     Int_t detector  = pTRDcluster->GetDetector();
     Int_t iPlane    = fGeom->GetPlane(detector);
-    if (iPlane >= AliESDtrack::kNPlane) {
-      AliError(Form("Wrong plane %d",iPlane));
-      continue;
-    }
-    Int_t iSlice    = tb * AliESDtrack::kNSlice 
-                         / AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
-    if (iSlice >= AliESDtrack::kNSlice) {
-      AliError(Form("Wrong slice %d",iSlice));
-      continue;
-    }
-    clscharge[iPlane][iSlice] = clscharge[iPlane][iSlice] + charge;
+
+               // 2 slices method
+               Int_t iSlice = (tb < AliTRDtrack::kMidTimeBin) ? 0 : 1;
+               // weighted edep with <PH>
+               clscharge[iPlane][iSlice] += charge * fAveragePH[tb];
     if (charge > maxclscharge[iPlane]) {
       maxclscharge[iPlane] = charge;
       timebin[iPlane]      = tb;
@@ -3526,31 +3534,18 @@ void AliTRDtracker::CookdEdxTimBin(AliTRDtrack &TRDtrack)
     nCluster[iPlane][iSlice]++;
   } // End of loop over cluster
 
-  // Setting the fdEdxPlane and fTimBinPlane variabales 
-  Double_t totalCharge = 0.0;
-
+  // Setting the fdEdxPlane and fTimBinPlane variables
+  Double_t totalCharge = 0.;
   for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
-    for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
-      if (nCluster[iPlane][iSlice]) {
-        clscharge[iPlane][iSlice] /= nCluster[iPlane][iSlice];
-      }
+               for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
+      //if (nCluster[iPlane][iSlice]) {
+      //  clscharge[iPlane][iSlice] /= nCluster[iPlane][iSlice];
+      //}
       TRDtrack.SetPIDsignals(clscharge[iPlane][iSlice],iPlane,iSlice);
-      totalCharge = totalCharge+clscharge[iPlane][iSlice];
+      totalCharge += clscharge[iPlane][iSlice];
     }
-    TRDtrack.SetPIDTimBin(timebin[iPlane],iPlane);     
+    TRDtrack.SetPIDTimBin(timebin[iPlane],iPlane);
   }
-
-  // Still needed ????
-  //  Int_t i;
-  //  Int_t nc=TRDtrack.GetNumberOfClusters(); 
-  //  Float_t dedx=0;
-  //  for (i=0; i<nc; i++) dedx += TRDtrack.GetClusterdQdl(i);
-  //  dedx /= nc;
-  //  for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
-  //    TRDtrack.SetPIDsignals(dedx, iPlane);
-  //    TRDtrack.SetPIDTimBin(timbin[iPlane], iPlane);
-  //  }
-
 }
 
 //_____________________________________________________________________________
index 6d8edca..44fd604 100644 (file)
 
 /* $Id$ */
 
-///////////////////////////////////////////////////////////////////////////////
-//                                                                           //
-// Container for the distributions of dE/dx and the time bin of the          //
-// max. cluster for electrons and pions                                      //
-//                                                                           //
-// Author:                                                                   //
-//   Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de>                    //
-//                                                                           //
-///////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+//                                                                      //
+// Container for the distributions of dE/dx and the time bin of the     //
+// max. cluster for electrons and pions                                 //
+//                                                                      //
+// Authors:                                                             //
+//   Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de>               //
+//   Alex Bercuci <a.bercuci@gsi.de>                                    //
+//                                                                      //
+//////////////////////////////////////////////////////////////////////////
 
 #include <TH1F.h>
+#include <TH2F.h>
 #include <TFile.h>
+#include <TTree.h>
 #include <TROOT.h>
+#include <TPDGCode.h>
+#include <TParticle.h>
+#include <TPrincipal.h>
 
 #include "AliLog.h"
+#include "AliHeader.h"
+#include "AliGenEventHeader.h"
+#include "AliRun.h"
+#include "AliRunLoader.h"
+#include "AliStack.h"
 #include "AliPID.h"
+#include "AliESD.h"
+#include "AliESDtrack.h"
 
 #include "AliTRDCalPIDLQ.h"
+#include "AliTRDCalPIDLQRef.h"
+#include "AliTRDpidESD.h"
+#include "AliTRDcalibDB.h"
+#include "AliTRDtrack.h"
+#include "AliTRDgeometry.h"
 
 ClassImp(AliTRDCalPIDLQ)
 
 Char_t* AliTRDCalPIDLQ::fpartName[AliPID::kSPECIES] = {"electron", "muon", "pion", "kaon", "proton"};
+Char_t* AliTRDCalPIDLQ::fpartSymb[AliPID::kSPECIES] = {"EL", "MU", "PI", "KA", "PR"};
+//const Int_t AliTRDCalPIDLQ::fNMom = 11;
+//const Int_t AliTRDCalPIDLQ::fNLength = 4;
     
 //_________________________________________________________________________
 AliTRDCalPIDLQ::AliTRDCalPIDLQ()
-  :TNamed()
-  ,fNMom(0)
-  ,fTrackMomentum(0)
+  :TNamed("pid", "PID for TRD")
+  //,fNMom(0)
+  ,fTrackMomentum(0x0)
+  //,fNLength(0)
+  ,fTrackSegLength(0x0)
+  ,fNTimeBins(0)
   ,fMeanChargeRatio(0)
   ,fNbins(0)
   ,fBinSize(0)
-  ,fHistdEdx(0)
-  ,fHistTimeBin(0)
+  ,fHistdEdx(0x0)
+  ,fHistTimeBin(0x0)
+  ,fEstimator(0x0)
 {
   //
   //  The Default constructor
   //
-  
+
+       Init();
 }
 
 //_________________________________________________________________________
 AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t *name, const Text_t *title) 
   :TNamed(name,title)
-  ,fNMom(0)
-  ,fTrackMomentum(0)
+  //,fNMom(0)
+  ,fTrackMomentum(0x0)
+  //,fNLength(0)
+  ,fTrackSegLength(0x0)
+  ,fNTimeBins(0)
   ,fMeanChargeRatio(0)
   ,fNbins(0)
   ,fBinSize(0)
-  ,fHistdEdx(0)
-  ,fHistTimeBin(0)
+  ,fHistdEdx(0x0)
+  ,fHistTimeBin(0x0)
+  ,fEstimator(0x0)
 {
   //
   //  The main constructor
@@ -77,32 +107,24 @@ AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t *name, const Text_t *title)
 //_____________________________________________________________________________
 AliTRDCalPIDLQ::AliTRDCalPIDLQ(const AliTRDCalPIDLQ &c) 
   :TNamed(c)
-  ,fNMom(c.fNMom)
-  ,fTrackMomentum(0)
+  //,fNMom(c.fNMom)
+  ,fTrackMomentum(0x0)
+  //,fNLength(c.fNLength)
+  ,fTrackSegLength(0x0)
+  ,fNTimeBins(c.fNTimeBins)
   ,fMeanChargeRatio(c.fMeanChargeRatio)
   ,fNbins(c.fNbins)
   ,fBinSize(c.fBinSize)
-  ,fHistdEdx(0)
-  ,fHistTimeBin(0)
+  ,fHistdEdx(0x0)
+  ,fHistTimeBin(0x0)
+  ,fEstimator(0x0)
 {
   //
   // Copy constructor
   //
 
-  AliTRDCalPIDLQ& target = (AliTRDCalPIDLQ &) c;
+  if (this != &c) ((AliTRDCalPIDLQ &) c).Copy(*this);
   
-  target.fTrackMomentum = new Double_t[fNMom];
-  for (Int_t i=0; i<fNMom; ++i) {
-    target.fTrackMomentum[i] = fTrackMomentum[i];
-  }
-  if (fHistdEdx) {
-    target.fHistdEdx = (TObjArray*) fHistdEdx->Clone();
-  }  
-
-  if (fHistTimeBin) {
-    target.fHistTimeBin = (TObjArray*) fHistTimeBin->Clone();
-  }
-
 }
 
 //_________________________________________________________________________
@@ -125,17 +147,27 @@ void AliTRDCalPIDLQ::CleanUp()
 
   if (fHistdEdx) {
     delete fHistdEdx;
-    fHistdEdx = 0;
+    fHistdEdx = 0x0;
   }
   
   if (fHistTimeBin) {
     delete fHistTimeBin;
-    fHistTimeBin = 0;
+    fHistTimeBin = 0x0;
   }
 
   if (fTrackMomentum) {
     delete[] fTrackMomentum;
-    fTrackMomentum = 0;
+    fTrackMomentum = 0x0;
+  }
+
+  if (fTrackSegLength) {
+    delete[] fTrackSegLength;
+    fTrackSegLength = 0x0;
+  }
+
+  if (fEstimator) {
+    delete fEstimator;
+    fEstimator = 0x0;
   }
 
 }
@@ -163,16 +195,23 @@ void AliTRDCalPIDLQ::Copy(TObject &c) const
   
   target.CleanUp();
   
-  target.fNMom            = fNMom;
   target.fNbins           = fNbins;
   target.fBinSize         = fBinSize;
   target.fMeanChargeRatio = fMeanChargeRatio;
-  
+  target.fNTimeBins       = fNTimeBins;
+
+  //target.fNMom            = fNMom;
   target.fTrackMomentum = new Double_t[fNMom];
   for (Int_t i=0; i<fNMom; ++i) {
     target.fTrackMomentum[i] = fTrackMomentum[i];
   }
 
+  //target.fNLength         = fNLength;
+  target.fTrackSegLength = new Double_t[fNLength];
+  for (Int_t i=0; i<fNLength; ++i) {
+    target.fTrackSegLength[i] = fTrackSegLength[i];
+  }
+
   if (fHistdEdx) {
     target.fHistdEdx = (TObjArray*) fHistdEdx->Clone();
   }
@@ -180,6 +219,8 @@ void AliTRDCalPIDLQ::Copy(TObject &c) const
     target.fHistTimeBin = (TObjArray*) fHistTimeBin->Clone();
   }
 
+  target.fEstimator = new AliTRDCalPIDLQRef(*fEstimator);
+
   TObject::Copy(c);
 
 }
@@ -191,33 +232,54 @@ void AliTRDCalPIDLQ::Init()
   // Initialization
   //
 
-  const Int_t kNMom = 11;
+  fNMom    = 11;
+  fNLength =  4;
 
-  fNMom = kNMom;
   fTrackMomentum = new Double_t[fNMom];
-  Double_t trackMomentum[kNMom] = {  0.6,  0.8,  1.0,  1.5,  2.0
-                                  ,  3.0,  4.0,  5.0,  6.0,  8.0
-                                  , 10.0 };
-  for (Int_t imom = 0; imom < kNMom; imom++) {
-    fTrackMomentum[imom] = trackMomentum[imom];
-  }  
-
-  fHistdEdx    = new TObjArray(AliPID::kSPECIES * fNMom);
+  fTrackMomentum[0] = 0.6;
+  fTrackMomentum[1] = 0.8;
+  fTrackMomentum[2] = 1.0;
+  fTrackMomentum[3] = 1.5;
+  fTrackMomentum[4] = 2.0;
+  fTrackMomentum[5] = 3.0;
+  fTrackMomentum[6] = 4.0;
+  fTrackMomentum[7] = 5.0;
+  fTrackMomentum[8] = 6.0;
+  fTrackMomentum[9] = 8.0;
+  fTrackMomentum[10] = 10.0;
+  
+  fTrackSegLength = new Double_t[fNLength];
+  fTrackSegLength[0] = 3.7;
+  fTrackSegLength[1] = 3.9;
+  fTrackSegLength[2] = 4.2;
+  fTrackSegLength[3] = 5.0;
+
+  fHistdEdx    = new TObjArray(AliPID::kSPECIES * fNMom/* * fNLength*/);
   fHistdEdx->SetOwner();
-  fHistTimeBin = new TObjArray(AliPID::kSPECIES * fNMom);
+  fHistTimeBin = new TObjArray(2 * fNMom);
   fHistTimeBin->SetOwner();  
-  
+       // Initialization of estimator at object instantiation because late
+       // initialization in function GetProbability() is not working due to
+       // constantness of this function. 
+       fEstimator = new AliTRDCalPIDLQRef();
+
+       // Number of Time bins
+       AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
+       if(!calibration){
+               AliWarning("No AliTRDcalibDB available. Using 22 time bins.");
+               fNTimeBins = 22;
+       } else fNTimeBins = calibration->GetNumberOfTimeBins();
+       
   // ADC Gain normalization
   fMeanChargeRatio = 1.0;
   
   // Number of bins and bin size
   fNbins   = 0;
   fBinSize = 0.0;
-
 }
 
 //_________________________________________________________________________
-Bool_t AliTRDCalPIDLQ::ReadData(Char_t *responseFile)
+Bool_t AliTRDCalPIDLQ::ReadReferences(Char_t *responseFile)
 {
   //
   // Read the TRD dEdx histograms.
@@ -232,211 +294,462 @@ Bool_t AliTRDCalPIDLQ::ReadData(Char_t *responseFile)
   gROOT->cd();
 
   // Read histograms
-  Char_t text[200];
-  for (Int_t particle = 0; particle < AliPID::kSPECIES; ++particle)
-  {
-    Char_t* particleKey = "";
-    switch (particle)
-    {
-      case AliPID::kElectron: particleKey = "EL"; break;
-      case AliPID::kPion: particleKey = "PI"; break;
-      case AliPID::kMuon: particleKey = "MU"; break;
-      case AliPID::kKaon: particleKey = "KA"; break;
-      case AliPID::kProton: particleKey = "PR"; break;
-    }
-    
-    for (Int_t imom = 0; imom < fNMom; imom++) 
-    {
-      sprintf(text, "h1dEdx%s%01d", particleKey, imom+1);
-      TH1F* hist = (TH1F*)histFile->Get(text)->Clone();
-      hist->Scale(1.0/hist->Integral());
-      fHistdEdx->AddAt(hist, GetHistID(particle, imom));
-  
-      if (particle == AliPID::kElectron || particle == AliPID::kPion)
-      {
-        sprintf(text,"h1MaxTimBin%s%01d", particleKey, imom+1);
-        TH1F* hist = (TH1F*)histFile->Get(text)->Clone();
-        hist->Scale(1.0/hist->Integral());
-        fHistTimeBin->AddAt(hist, GetHistID(particle,imom));
-      }
-    }
+  for (Int_t iparticle = 0; iparticle < AliPID::kSPECIES; iparticle++){
+    for (Int_t imom = 0; imom < fNMom; imom++){
+                       TH2D* hist = (TH2D*)histFile->Get(Form("h2dEdx%s%02d", fpartSymb[iparticle], imom/*, ilength*/))->Clone();
+                       hist->Scale(1./hist->Integral());
+                       fHistdEdx->AddAt(hist, GetHistID(iparticle, imom));
+
+                       if (iparticle != AliPID::kElectron && iparticle != AliPID::kPion) continue;
+
+                       TH1F* ht = (TH1F*)histFile->Get(Form("h1MaxTB%s%02d", fpartSymb[iparticle], imom))->Clone();
+                       if(ht->GetNbinsX() != fNTimeBins) AliWarning(Form("The number of time bins %d defined in h1MaxTB%s%02d differs from calibration value of %d. This may lead to erroneous results.", ht->GetNbinsX(), fpartSymb[iparticle], imom, fNTimeBins));
+                       ht->Scale(1./ht->Integral());
+                       fHistTimeBin->AddAt(ht, ((iparticle==AliPID::kElectron)?0:1)*fNMom + imom);
+               }
   }
   
   histFile->Close();
   delete histFile;
   
   // Number of bins and bin size
-  TH1F* hist = (TH1F*) fHistdEdx->At(GetHistID(AliPID::kPion, 1));
-  fNbins   = hist->GetNbinsX();
-  fBinSize = hist->GetBinWidth(1);
+  //TH1F* hist = (TH1F*) fHistdEdx->At(GetHistID(AliPID::kPion, 1));
+  //fNbins   = hist->GetNbinsX();
+  //fBinSize = hist->GetBinWidth(1);
   
   return kTRUE;
 
 }
 
-//_________________________________________________________________________
-Double_t  AliTRDCalPIDLQ::GetMean(Int_t k, Int_t ip) const
-{
-  //
-  // Gets mean of de/dx dist. of e
-  //
-
-  AliInfo(Form("Mean for particle = %s and momentum = %.2f is:\n"
-              ,fpartName[k]
-              ,fTrackMomentum[ip]));
-  if (k < 0 || k > AliPID::kSPECIES) {
-    return 0;
-  }
-
-  return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->GetMean();
-
-}
+// //_________________________________________________________________________
+// Double_t  AliTRDCalPIDLQ::GetMean(Int_t k, Int_t ip) const
+// {
+//   //
+//   // Gets mean of de/dx dist. of e
+//   //
+// 
+//   AliInfo(Form("Mean for particle = %s and momentum = %.2f is:\n"
+//               ,fpartName[k]
+//               ,fTrackMomentum[ip]));
+//   if (k < 0 || k > AliPID::kSPECIES) {
+//     return 0;
+//   }
+// 
+//   return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->GetMean();
+// 
+// }
+// 
+// //_________________________________________________________________________
+// Double_t  AliTRDCalPIDLQ::GetNormalization(Int_t k, Int_t ip) const
+// {
+//   //
+//   // Gets Normalization of de/dx dist. of e
+//   //
+// 
+//   AliInfo(Form("Normalization for particle = %s and momentum = %.2f is:\n"
+//               ,fpartName[k]
+//               ,fTrackMomentum[ip]));
+//   if (k < 0 || k > AliPID::kSPECIES) {
+//     return 0;
+//   }
+//   
+//   return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->Integral();
+// 
+// }
 
 //_________________________________________________________________________
-Double_t  AliTRDCalPIDLQ::GetNormalization(Int_t k, Int_t ip) const
-{
-  //
-  // Gets Normalization of de/dx dist. of e
-  //
-
-  AliInfo(Form("Normalization for particle = %s and momentum = %.2f is:\n"
-              ,fpartName[k]
-              ,fTrackMomentum[ip]));
-  if (k < 0 || k > AliPID::kSPECIES) {
-    return 0;
-  }
-  
-  return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->Integral();
-
-}
-
-//_________________________________________________________________________
-TH1F* AliTRDCalPIDLQ::GetHistogram(Int_t k, Int_t ip) const
+TH1* AliTRDCalPIDLQ::GetHistogram(Int_t k, Int_t ip/*, Int_t il*/) const
 {
   //
   // Returns one selected dEdx histogram
   //
 
-  AliInfo(Form("Histogram for particle = %s and momentum = %.2f is:\n"
-              ,fpartName[k]
-              ,fTrackMomentum[ip]));
-  if (k < 0 || k > AliPID::kSPECIES) {
-    return 0;
-  }
+  if (k < 0 || k >= AliPID::kSPECIES) return 0x0;
+       if(ip<0 || ip>= fNMom ) return 0x0;
+
+       AliInfo(Form("Retrive dEdx histogram for %s of %5.2f GeV/c", fpartName[k], fTrackMomentum[ip]));
   
-  return (TH1F*) fHistdEdx->At(GetHistID(k,ip));
+  return (TH1*)fHistdEdx->At(GetHistID(k, ip));
 
 }
 
 //_________________________________________________________________________
-TH1F* AliTRDCalPIDLQ::GetHistogramT(Int_t k, Int_t ip) const
+TH1* AliTRDCalPIDLQ::GetHistogramT(Int_t k, Int_t ip) const
 {
   //
   // Returns one selected time bin max histogram
   //
 
-  AliInfo(Form("Histogram for particle = %s and momentum = %.2f is:\n"
-              ,fpartName[k]
-              ,fTrackMomentum[ip]));
-  if (k < 0 || k > AliPID::kSPECIES)
-    return 0;
-  
-  return (TH1F*) fHistTimeBin->At(GetHistID(k,ip));
+  if (k < 0 || k >= AliPID::kSPECIES) return 0x0;
+       if(ip<0 || ip>= fNMom ) return 0x0;
+         
+       AliInfo(Form("Retrive MaxTB histogram for %s of %5.2f GeV/c", fpartName[k], fTrackMomentum[ip]));
 
+       return (TH1*)fHistTimeBin->At(((k==AliPID::kElectron)?0:1)*fNMom+ip);
 }
 
 //_________________________________________________________________________
-Double_t AliTRDCalPIDLQ::GetProbability(Int_t k, Double_t mom, Double_t dedx1) const
+Double_t AliTRDCalPIDLQ::GetProbability(Int_t spec, Double_t mom, Double_t *dedx, Double_t length) const 
 {
   //
-  // Gets the Probability of having dedx at a given momentum (mom)
-  // and particle type k (0 for e) and (2 for pi)
-  // from the precalculated de/dx distributions 
+       // Core function of AliTRDCalPIDLQ 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
   //
-  
-  Double_t dedx   = dedx1/fMeanChargeRatio;
-  Int_t    iEnBin = ((Int_t) (dedx/fBinSize+1));
-  if(iEnBin > fNbins) iEnBin = fNbins;
 
-  if (k < 0 || k > AliPID::kSPECIES) {
-    return 1;
-  }
-  
-  TH1F* hist1 = 0;
-  TH1F* hist2 = 0;
-  Double_t mom1 = 0;
-  Double_t mom2 = 0;
-  
-  // Lower limit
-  if (mom<=fTrackMomentum[0])  {
-    hist1 = (TH1F*) fHistdEdx->At(GetHistID(k,1));
-    hist2 = (TH1F*) fHistdEdx->At(GetHistID(k,0));
-    mom1 = fTrackMomentum[1];
-    mom2 = fTrackMomentum[0];
-  }
-    
-  // Upper Limit
-  if(mom>=fTrackMomentum[fNMom-1]) {
-    hist2 = (TH1F*) fHistdEdx->At(GetHistID(k,fNMom-1));
-    hist1 = (TH1F*) fHistdEdx->At(GetHistID(k,fNMom-2));
-    mom2 = fTrackMomentum[fNMom-1];
-    mom1 = fTrackMomentum[fNMom-2];
-  }
-    
-  // In the range
-  for (Int_t ip=1; ip<fNMom; ip++) {
-    if ((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
-      hist1 = (TH1F*) fHistdEdx->At(GetHistID(k,ip));
-      hist2 = (TH1F*) fHistdEdx->At(GetHistID(k,ip-1));
-      mom1 = fTrackMomentum[ip];
-      mom2 = fTrackMomentum[ip-1];
-    }
-  }
-  
-  Double_t slop = (hist1->GetBinContent(iEnBin) - hist2->GetBinContent(iEnBin)) 
-                / (mom1 - mom2);
-  return hist2->GetBinContent(iEnBin) + slop * (mom - mom2);
+       if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
+               
+       //Double_t dedx   = dedx1/fMeanChargeRatio;
+       
+       // find the interval in momentum and track segment length which applies for this data
+       Int_t ilength = 1;
+  while(ilength<fNLength-1 && length>fTrackSegLength[ilength]){
+               ilength++;
+       }
+       Int_t imom = 1;
+  while(imom<fNMom-1 && mom>fTrackMomentum[imom]) imom++;
+       
+       Int_t nbinsx, nbinsy;
+       TAxis *ax = 0x0, *ay = 0x0;
+       Double_t LQ1, LQ2;
+       Double_t mom1 = fTrackMomentum[imom-1], mom2 = fTrackMomentum[imom];
+       TH2 *hist = 0x0;
+       if(!(hist = (TH2D*)fHistdEdx->At(GetHistID(spec, imom-1/*, ilength*/)))){
+               AliInfo(Form("Looking for spec(%d) mom(%f) Ex(%f) Ey(%f) length(%f)", spec, mom, dedx[0], dedx[1], length));
+               AliError(Form("EHistogram id %d not found in DB.", GetHistID(spec, imom-1)));
+               return 0.;
+       }
+       ax = hist->GetXaxis(); nbinsx = ax->GetNbins();
+       ay = hist->GetYaxis(); nbinsy = ay->GetNbins();
+       Bool_t kX = (dedx[0] < ax->GetBinUpEdge(nbinsx));
+       Bool_t kY = (dedx[1] < ay->GetBinUpEdge(nbinsy));
+       if(kX)
+               if(kY) LQ1 = hist->GetBinContent( hist->FindBin(dedx[0], dedx[1])); //fEstimator->Estimate2D2(hist, (Float_t&)dedx[0], (Float_t&)dedx[1]);
+               else LQ1 = hist->GetBinContent(ax->FindBin(dedx[0]), nbinsy);
+       else
+               if(kY) LQ1 = hist->GetBinContent(nbinsx, ay->FindBin(dedx[1]));
+               else LQ1 = hist->GetBinContent(nbinsx, nbinsy);
+
+
+       if(!(hist = (TH2D*)fHistdEdx->At(GetHistID(spec, imom/*, ilength*/)))){
+               AliInfo(Form("Looking for spec(%d) mom(%f) Ex(%f) Ey(%f) length(%f)", spec, mom, dedx[0], dedx[1], length));
+               AliError(Form("EHistogram id %d not found in DB.", GetHistID(spec, imom)));
+               return LQ1;
+       }
+       if(kX)
+               if(kY) LQ2 = hist->GetBinContent( hist->FindBin(dedx[0], dedx[1])); //fEstimator->Estimate2D2(hist, (Float_t&)dedx[0], (Float_t&)dedx[1]);
+               else LQ2 = hist->GetBinContent(ax->FindBin(dedx[0]), nbinsy);
+       else
+               if(kY) LQ2 = hist->GetBinContent(nbinsx, ay->FindBin(dedx[1]));
+               else LQ2 = hist->GetBinContent(nbinsx, nbinsy);
+
+       
+       // return interpolation over momentum binning
+  return LQ1 + (LQ2 - LQ1)*(mom - mom1)/(mom2 - mom1);
 
 }
 
 //_________________________________________________________________________
-Double_t AliTRDCalPIDLQ::GetProbabilityT(Int_t k, Double_t mom, Int_t timbin) const
+Double_t AliTRDCalPIDLQ::GetProbabilityT(Int_t spec, Double_t mom, Int_t timbin) const
 {
   //
   // Gets the Probability of having timbin at a given momentum (mom)
-  // and particle type k (0 for e) and (2 for pi)
+  // and particle type (spec) (0 for e) and (2 for pi)
   // from the precalculated timbin distributions 
   //
   
-  if (timbin<=0) {
-    return 0.0;
-  }
+       if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
+  if (timbin<0 || timbin >= fNTimeBins) return 0.;
 
   Int_t iTBin = timbin+1;
   
   // Everything which is not an electron counts as a pion for time bin max
-  if (k != AliPID::kElectron) {
-    k = AliPID::kPion;
-  }
+  //if(spec != AliPID::kElectron) spec = AliPID::kPion;
+  
 
-  if (mom<=fTrackMomentum[0]) {
-    return ((TH1F*) fHistTimeBin->At(GetHistID(k,0)))->GetBinContent(iTBin);
-  }
-  if (mom>=fTrackMomentum[fNMom-1]) { 
-    return ((TH1F*) fHistTimeBin->At(GetHistID(k,fNMom-1)))->GetBinContent(iTBin);
-  }
+  
+       Int_t imom = 1;
+  while(imom<fNMom-1 && mom>fTrackMomentum[imom]) imom++;
+
+       Double_t mom1 = fTrackMomentum[imom-1], mom2 = fTrackMomentum[imom];
+       TH1F *hist = 0x0;
+       if(!(hist = (TH1F*) fHistTimeBin->At(((spec==AliPID::kElectron)?0:1)*fNMom+imom-1))){
+               AliInfo(Form("Looking for spec(%d) mom(%f) timbin(%d)", spec, mom, timbin));
+               AliError(Form("THistogram id %d not found in DB.", ((spec==AliPID::kElectron)?0:1)*fNMom+imom-1));
+               return 0.;
+       }
+       Double_t LQ1 = hist->GetBinContent(iTBin);
+
+       if(!(hist = (TH1F*) fHistTimeBin->At(((spec==AliPID::kElectron)?0:1)*fNMom+imom))){
+               AliInfo(Form("Looking for spec(%d) mom(%f) timbin(%d)", spec, mom, timbin));
+               AliError(Form("THistogram id %d not found in DB.", ((spec==AliPID::kElectron)?0:1)*fNMom+imom));
+               return LQ1;
+       }
+       Double_t LQ2 = hist->GetBinContent(iTBin);
+
+       // return interpolation over momentum binning
+  return LQ1 + (LQ2 - LQ1)*(mom - mom1)/(mom2 - mom1);
+}
 
-  for (Int_t ip=1; ip<fNMom; ip++) {
-    if ((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
-      Double_t slop = (((TH1F*) fHistTimeBin->At(GetHistID(k,ip)))->GetBinContent(iTBin) 
-                     - ((TH1F*) fHistTimeBin->At(GetHistID(k,ip-1)))->GetBinContent(iTBin)) 
-                    / (fTrackMomentum[ip] - fTrackMomentum[ip-1]);
-      // Linear interpolation
-      return ((TH1F*) fHistTimeBin->At(GetHistID(k,ip-1)))->GetBinContent(iTBin) 
-              + slop * (mom - fTrackMomentum[ip-1]);
-    }
+//__________________________________________________________________
+Bool_t AliTRDCalPIDLQ::WriteReferences(Char_t *File, Char_t *dir)
+{
+       // Build, Fill and write to file the histograms used for PID.
+       // The simulations are looked in the
+       // directories with the general form Form("p%3.1f", momentum)
+       // starting from dir (default .). Here momentum belongs to the list
+       // of known momentum to PID (see fTrackMomentum).
+       // The output histograms are
+       // written to the file "File" in cwd (default
+       // TRDPIDHistograms.root). In order to build a DB entry
+       // consider running $ALICE_ROOT/Cal/AliTRDCreateDummyCDB.C
+       // 
+       // Author:
+       // Alex Bercuci (A.Bercuci@gsi.de)
+
+       const Int_t nDirs = 1;
+  Int_t partCode[AliPID::kSPECIES] =
+    {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};
+       
+       // minimal test of simulation location
+       TFile *f = new TFile(Form("%s/p%3.1f/galice.root", dir, 2.));
+       if(!f || f->IsZombie()){
+               AliError(Form("Could not access file galice in directry \"%s/p%3.1f\". Please check the location and try again.", dir, 2.));
+               return kFALSE;
+       }
+       f->Close(); delete f;
+       f = new TFile(Form("%s/p%3.1f/AliESDs.root", dir, 2.));
+       if(!f || f->IsZombie()){
+               AliError(Form("Could not access file AliESDs in directry \"%s/p%3.1f\". Please check the location and try again.", dir, 2.));
+               return kFALSE;
+       }
+       f->Close(); delete f;
+       
+
+       // Init statistics
+       Int_t nPart[AliPID::kSPECIES], nTotPart;
+       printf("P[GeV/c] ");
+       for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) printf(" %s[%%] ", fpartSymb[ispec]);
+       printf("\n-----------------------------------------------\n");
+       
+       
+
+       // Build PID reference histograms and reference object
+       const Int_t color[] = {4, 3, 2, 7, 6};
+       for (Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
+               if (ispec != AliPID::kElectron && ispec != AliPID::kPion) continue;
+       
+               h1MaxTB[(ispec>0)?1:0] = new TH1F(Form("h1%s", fpartSymb[ispec]), "", fNTimeBins, -.5, fNTimeBins-.5);
+               h1MaxTB[(ispec>0)?1:0]->SetLineColor(color[ispec]);
   }
-  
-  return -1;
+       AliTRDCalPIDLQRef ref;
+       
+       // momentum loop
+       AliRunLoader *fRunLoader = 0x0;
+       TFile *esdFile = 0x0;
+       TTree *esdTree = 0x0;
+       AliESD *esd = 0x0;
+       AliESDtrack *esdTrack = 0x0;
+       for (Int_t imom = 0; imom < fNMom; imom++) {
+       //for (Int_t imom = 4; imom < 5; imom++) {
+               ref.Reset();
+               
+               for(Int_t idir = 0; idir<nDirs; idir++){
+                       // open run loader and load gAlice, kinematics and header
+                       fRunLoader = AliRunLoader::Open(Form("%s/p%3.1f/galice.root", dir, fTrackMomentum[imom]));
+                       if (!fRunLoader) {
+                               AliError(Form("Getting run loader for momentum %3.1f failed.", fTrackMomentum[imom]));
+                               return kFALSE;
+                       }
+                       TString s; s.Form("%s/p%3.1f/", dir, fTrackMomentum[imom]);
+                       fRunLoader->SetDirName(s);
+                       fRunLoader->LoadgAlice();
+                       gAlice = fRunLoader->GetAliRun();
+                       if (!gAlice) {
+                               AliError(Form("galice object not found for momentum %3.1f.", fTrackMomentum[imom]));
+                               return kFALSE;
+                       }
+                       fRunLoader->LoadKinematics();
+                       fRunLoader->LoadHeader();
+       
+                       // open the ESD file
+                       esdFile = TFile::Open(Form("%s/p%3.1f/AliESDs.root", dir, fTrackMomentum[imom]));
+                       if (!esdFile || esdFile->IsZombie()) {
+                               AliError(Form("Opening ESD file failed for momentum", fTrackMomentum[imom]));
+                               return kFALSE;
+                       }
+                       esd = new AliESD;
+                       esdTree = (TTree*)esdFile->Get("esdTree");
+                       if (!esdTree) {
+                               AliError(Form("ESD tree not found for momentum %3.1f.", fTrackMomentum[imom]));
+                               return kFALSE;
+                       }
+                       esdTree->SetBranchAddress("ESD", &esd);
+                       nTotPart = 0;
+                       for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) nPart[ispec] = 0;
+       
+                       
+                       // Event loop
+                       for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
+                               fRunLoader->GetEvent(iEvent);
+       
+                               // read stack info
+                               AliStack* stack = gAlice->Stack();
+                               TArrayF vertex(3);
+                               fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
+                                                       
+                               // Load event summary data
+                               esdTree->GetEvent(iEvent);
+                               if (!esd) {
+                                       AliWarning(Form("ESD object not found for event %d. [@ momentum %3.1f]", iEvent, imom));
+                                       continue;
+                               }
+       
+                               for(Int_t iTrack=0; iTrack<esd->GetNumberOfTracks(); iTrack++){
+                                       esdTrack = esd->GetTrack(iTrack);
+       
+                                       if(!AliTRDpidESD::CheckTrack(esdTrack)) continue;
+                                       //if((esdTrack->GetStatus() & AliESDtrack::kITSrefit) == 0) continue;
+                                       //if(esdTrack->GetConstrainedChi2() > 1E9) continue;
+                                       //if ((esdTrack->GetStatus() & AliESDtrack::kESDpid) == 0) continue;
+                                       if (esdTrack->GetTRDsignal() == 0.) continue;
+       
+                                       // read MC info
+                                       Int_t label = esdTrack->GetLabel();
+                                       if(label<0) continue;
+                                       if (label > stack->GetNtrack()) continue;     // background
+                                       TParticle* particle = stack->Particle(label);
+                                       if(!particle){
+                                               AliWarning(Form("Retriving particle with index %d from AliStack failed. [@ event %d track %d]", label, iEvent, iTrack));
+                                               continue;
+                                       }
+                                       if(particle->Pt() < 1.E-3) continue;
+                                       //      if (TMath::Abs(particle->Eta()) > 0.3) continue;
+                                       TVector3 dVertex(particle->Vx() - vertex[0],
+                                                                               particle->Vy() - vertex[1],
+                                                                               particle->Vz() - vertex[2]);
+                                       if (dVertex.Mag() > 1.E-4){
+                                               //AliInfo(Form("Particle with index %d generated too far from vertex. Skip from analysis. Details follows. [@ event %d track %d]", label, iEvent, iTrack));
+                                               //particle->Print();
+                                               continue;
+                                       }
+                                       Int_t iGen = -1;
+                                       for (Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++)
+                                               if(TMath::Abs(particle->GetPdgCode()) == partCode[ispec]){
+                                                       iGen = ispec;
+                                                       break;
+                                               }
+                                       if(iGen<0) continue;
+       
+                                       nPart[iGen]++; nTotPart++;
+                                       
+                                       Float_t mom, length;
+                                       Double_t dedx[AliTRDtrack::kNslice], dEdx;
+                                       Int_t timebin;
+                                       for (Int_t iPlane=0; iPlane<AliTRDgeometry::kNplan; iPlane++){
+                                               // read data for track segment
+                                               for(int iSlice=0; iSlice<AliTRDtrack::kNslice; iSlice++)
+                                                       dedx[iSlice] = esdTrack->GetTRDsignals(iPlane, iSlice);
+                                               dEdx    = esdTrack->GetTRDsignals(iPlane, -1);
+                                               timebin = esdTrack->GetTRDTimBin(iPlane);
+                       
+                                               // check data
+                                               if ((dEdx <=  0.) || (timebin <= -1.)) continue;
+                       
+                                               // retrive kinematic info for this track segment
+                                               if(!AliTRDpidESD::GetTrackSegmentKine(esdTrack, iPlane, mom, length)) continue;
+                                               
+                                               // find segment length and momentum bin
+                                               Int_t jmom = 1, refMom = -1;
+                                               while(jmom<fNMom-1 && mom>fTrackMomentum[jmom]) jmom++;
+                                               if(TMath::Abs(fTrackMomentum[jmom-1] - mom) < fTrackMomentum[jmom-1] * .2) refMom = jmom-1;
+                                               else if(TMath::Abs(fTrackMomentum[jmom] - mom) < fTrackMomentum[jmom] * .2) refMom = jmom;
+                                               if(refMom<0){
+                                                       AliInfo(Form("Momentum at plane %d entrance not in momentum window. [@ event %d track %d]", iPlane, iEvent, iTrack));
+                                                       continue;
+                                               }
+                                               /*while(jleng<fNLength-1 && length>fTrackSegLength[jleng]) jleng++;*/
+                                               
+                                               // this track segment has fulfilled all requierments
+                                               //nPlanePID++;
+
+                                               if(dedx[0] > 0. && dedx[1] > 0.){
+                                                       dedx[0] = log(dedx[0]); dedx[1] = log(dedx[1]);
+                                                       ref.GetPrincipal(iGen)->AddRow(dedx);
+                                               }
+                                               h1MaxTB[(iGen>0)?1:0]->Fill(timebin);
+                                       } // end plane loop
+                               } // end track loop
+                       } // end events loop
+                       
+                       delete esd; esd = 0x0;
+                       esdFile->Close();
+                       delete esdFile; esdFile = 0x0;
+       
+                       fRunLoader->UnloadHeader();
+                       fRunLoader->UnloadKinematics();
+                       delete fRunLoader; fRunLoader = 0x0;
+               } // end directory loop
+               
+               // use data to prepare references
+               ref.Prepare2DReferences();
+               // save this dEdx references
+               ref.SaveReferences(imom, File);
+               // save MaxTB references
+               SaveMaxTimeBin(imom, File);
+
+                       
+               // print momentum statistics
+               printf("  %3.1f  ", fTrackMomentum[imom]);
+               for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) printf(" %5.2f ", 100.*nPart[ispec]/nTotPart);
+               printf("\n");
+       } // end momentum loop
+       
+       TFile *fSave = 0x0;
+       TListIter it((TList*)gROOT->GetListOfFiles());
+       while((fSave=(TFile*)it.Next()))
+               if(strcmp(File, fSave->GetName())==0) break;
+
+       fSave->cd();
+       fSave->Close();
+       delete fSave;
+
+       return kTRUE;
+}
 
+//__________________________________________________________________
+void   AliTRDCalPIDLQ::SaveMaxTimeBin(const Int_t mom, const char *fn)
+{
+  //
+  // Save the histograms
+  //
+
+       TFile *fSave = 0x0;
+       TListIter it((TList*)gROOT->GetListOfFiles());
+       TDirectory *pwd = gDirectory;
+       Bool_t kFOUND = kFALSE;
+       while((fSave=(TFile*)it.Next()))
+               if(strcmp(fn, fSave->GetName())==0){
+                       kFOUND = kTRUE;
+                       break;
+               }
+       if(!kFOUND) fSave = new TFile(fn, "RECREATE");
+       fSave->cd();
+
+       TH1 *h;
+       h = (TH1F*)h1MaxTB[0]->Clone(Form("h1MaxTBEL%02d", mom));
+       h->SetTitle(Form("Maximum Time Bin distribution for electrons @ %4.1f GeV", fTrackMomentum[mom]));
+       h->GetXaxis()->SetTitle("time [100 ns]");
+       h->GetYaxis()->SetTitle("Probability");
+       h->Write();
+
+       h = (TH1F*)h1MaxTB[1]->Clone(Form("h1MaxTBPI%02d", mom));
+       h->SetTitle(Form("Maximum Time Bin distribution for pions @ %4.1f GeV", fTrackMomentum[mom]));
+       h->GetXaxis()->SetTitle("time [100 ns]");
+       h->GetYaxis()->SetTitle("Probability");
+       h->Write();
+       
+       pwd->cd();
 }
+
index a1c9296..e47b291 100644 (file)
 //                                                                           //
 ///////////////////////////////////////////////////////////////////////////////
 
+#ifndef ROOT_TNamed
 #include <TNamed.h>
+#endif
 
-class TH1F;
+class TH1;
 class TObjArray;
+class AliTRDCalPIDLQRef;
 
 class AliTRDCalPIDLQ : public TNamed {
 
+ friend class AliTRDCalPIDLQRef;
+
  public:
 
-  AliTRDCalPIDLQ(); 
+  AliTRDCalPIDLQ();
   AliTRDCalPIDLQ(const Text_t *name, const Text_t *title);
   AliTRDCalPIDLQ(const AliTRDCalPIDLQ& pd);
-  virtual          ~AliTRDCalPIDLQ();              
-  AliTRDCalPIDLQ    &operator=(const AliTRDCalPIDLQ &c);
+  virtual        ~AliTRDCalPIDLQ();
+  AliTRDCalPIDLQ &operator=(const AliTRDCalPIDLQ &c);
+
+  virtual void Copy(TObject &c) const;
 
-  virtual void       Copy(TObject &c) const;
+  Bool_t       ReadReferences(Char_t *responseFile);
+  void         SetMeanChargeRatio(Double_t ratio)     { fMeanChargeRatio = ratio;  }
 
-          Bool_t     ReadData(Char_t *responseFile);     
-          void       SetMeanChargeRatio(Double_t ratio)     { fMeanChargeRatio = ratio;  }  
+  Double_t     GetMeanChargeRatio() const             { return fMeanChargeRatio;   }
+  Double_t     GetMomentum(Int_t ip) const            { return (ip<0 || ip>=fNMom)    ? -1. : fTrackMomentum[ip];  }
+  Double_t     GetLength(Int_t il) const              { return (il<0 || il>=fNLength) ? -1. : fTrackSegLength[il]; }
+  //Double_t   GetMean(Int_t iType, Int_t ip) const;
+  //Double_t   GetNormalization(Int_t iType, Int_t ip) const;
 
-          Double_t   GetMeanChargeRatio() const             { return fMeanChargeRatio;   } 
-          Double_t   GetMomentum(Int_t ip) const            { return fTrackMomentum[ip]; }
-          Double_t   GetMean(Int_t iType, Int_t ip) const;        
-          Double_t   GetNormalization(Int_t iType, Int_t ip) const;
+  TH1*         GetHistogram(Int_t iType, Int_t ip/*, Int_t il*/) const;
+  TH1*         GetHistogramT(Int_t iType, Int_t ip) const;
 
-          TH1F      *GetHistogram(Int_t iType, Int_t ip) const;
-          TH1F      *GetHistogramT(Int_t iType, Int_t ip) const;
+  Double_t     GetProbability(Int_t spec, Double_t mom, Double_t *dedx, Double_t length) const;
+  Double_t     GetProbabilityT(Int_t spec, Double_t mom, Int_t timbin) const;
+  //Int_t      GetNbins() const                       { return fNbins;   }
+  //Double_t   GetBinSize() const                     { return fBinSize; }
+  Bool_t       WriteReferences(Char_t *File="TRDPIDHistograms.root", Char_t *dir =".");
+
+ protected:
 
-          Double_t   GetProbability(Int_t iType, Double_t mom, Double_t dedx) const;
-          Double_t   GetProbabilityT(Int_t iType, Double_t mom, Int_t timbin) const;
-          Int_t      GetNbins() const                       { return fNbins;   }        
-          Double_t   GetBinSize() const                     { return fBinSize; } 
+  void         Init();
+  void         CleanUp();
 
+  inline Int_t GetHistID(Int_t part, Int_t mom/*, Int_t length=0*/) const { return part*fNMom + mom; }
+  void        SaveMaxTimeBin(const Int_t mom, const char *fn);
+                               
  protected:
 
-          void       Init();      
-          void       CleanUp();   
-  inline  Int_t      GetHistID(Int_t part, Int_t mom) const { return part*fNMom + mom; }
-    
-  static  Char_t    *fpartName[5];                 //! Names of particle species
-    
-          Int_t      fNMom;                        //  Number of momenta  
-          Double_t  *fTrackMomentum;               //[fNMom] Track momenta for which response functions are available
-          Double_t   fMeanChargeRatio;             //  Ratio of mean charge from real Det. to prob. dist.
-
-          Int_t      fNbins;                       //  Number of energy bins
-          Double_t   fBinSize;                     //  Size of energy bin
-    
-          TObjArray *fHistdEdx;                    //  Prob. of dEdx for 5 particles and for several momenta
-          TObjArray *fHistTimeBin;                 //  Prob. of max time bin for 5 particles and for several momenta
-    
-  ClassDef(AliTRDCalPIDLQ, 1)                      //  The TRD PID response container class
+  static  Char_t    *fpartName[5];     //! Names of particle species
+  static  Char_t    *fpartSymb[5];     //! Symbols of particle species  
+  Int_t              fNMom;            //  Number of momenta
+  Int_t              fNLength;         //  Number of track segment length intervals
+
+  Double_t          *fTrackMomentum;   //[fNMom]    Track momenta for which response functions are available
+  Double_t          *fTrackSegLength;  //[fNLength] Track segment lengths for which response functions are available
+  Int_t              fNTimeBins;       //  Number of time bins
+
+  Double_t           fMeanChargeRatio; //  Ratio of mean charge from real Det. to prob. dist.
+  Int_t              fNbins;           //  Number of energy bins
+  Double_t           fBinSize;         //  Size of energy bin
+
+  TObjArray         *fHistdEdx;        //  Prob. of dEdx for 5 particles and for several momenta
+  TObjArray         *fHistTimeBin;     //  Prob. of max time bin for 5 particles and for several momenta
+
+ private:
+
+  TH1              *h1MaxTB[2];      //!
+  AliTRDCalPIDLQRef *fEstimator;      //!
+       
+  ClassDef(AliTRDCalPIDLQ, 2)         //   The TRD PID response container class
 
 };
 
diff --git a/TRD/Cal/AliTRDCalPIDLQRef.cxx b/TRD/Cal/AliTRDCalPIDLQRef.cxx
new file mode 100644 (file)
index 0000000..3ad3de2
--- /dev/null
@@ -0,0 +1,429 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  TRD calibration class for 2-dim PID reference histograms                 //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include <TROOT.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TH2D.h>
+#include <TH2I.h>
+#include <TH3D.h>
+#include <TPrincipal.h>
+#include <TLinearFitter.h>
+#include <TVectorT.h>
+#include <TCanvas.h>
+#include <TEllipse.h>
+#include <TMarker.h>
+
+#include "AliLog.h"
+#include "AliPID.h"
+#include "AliESD.h"
+
+#include "AliTRDCalPIDLQ.h"
+#include "AliTRDCalPIDLQRef.h"
+
+ClassImp(AliTRDCalPIDLQRef)
+
+//__________________________________________________________________
+AliTRDCalPIDLQRef::AliTRDCalPIDLQRef() 
+  :TObject()
+  ,fFitter2D2(0x0)
+  ,fFitter2D1(0x0)
+{
+  //
+  // AliTRDCalPIDLQRef default constructor
+  //
+
+       // histogram settings
+       for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
+               h2dEdx[ispec] = 0x0;
+               fPrinc[ispec] = 0x0;
+       }
+       // build parabolic 2D fitter
+       fFitter2D2 = new TLinearFitter(6, "1++x++y++x*x++y*y++x*y");
+       fFitter2D1 = new TLinearFitter(3, "1++x++y");
+
+}
+
+//__________________________________________________________________
+AliTRDCalPIDLQRef::AliTRDCalPIDLQRef(const AliTRDCalPIDLQRef &ref) 
+  :TObject()
+  ,fFitter2D2(0x0)
+  ,fFitter2D1(0x0)
+{
+  //
+  // AliTRDCalPIDLQRef copy constructor
+  // 
+
+       for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
+               if(ref.h2dEdx[ispec]){
+                       h2dEdx[ispec] = new  TH2D((TH2D&)*(ref.h2dEdx[ispec]));
+               } else h2dEdx[ispec] = 0x0;
+               fPrinc[ispec] = 0x0;
+       }
+       
+       // use the argument constructor because copy constructor is not
+       // working as for the ROOT version 5.15/07
+       fFitter2D2 = new TLinearFitter(6, "1++x++y++x*x++y*y++x*y");
+       fFitter2D1 = new TLinearFitter(3, "1++x++y");
+
+}
+
+//__________________________________________________________________
+AliTRDCalPIDLQRef::~AliTRDCalPIDLQRef()
+{
+  //
+  // AliTRDCalPIDQRef destructor
+  //
+
+       for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
+               if(h2dEdx[ispec]) delete h2dEdx[ispec];
+               if(fPrinc[ispec]) delete fPrinc[ispec]; 
+       }       
+       delete fFitter2D1;
+       delete fFitter2D2;
+
+}
+
+//__________________________________________________________________
+AliTRDCalPIDLQRef& AliTRDCalPIDLQRef::operator=(const AliTRDCalPIDLQRef &ref)
+{
+  //
+  // Assignment operator
+  //
+
+       for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
+               if(ref.h2dEdx[ispec]) (ref.h2dEdx[ispec])->Copy(*h2dEdx[ispec]);
+               fPrinc[ispec] = 0x0;
+       }
+       return *this;
+
+}
+
+//__________________________________________________________________
+Double_t AliTRDCalPIDLQRef::Estimate2D2(TH2 *h, Float_t &x, Float_t &y)
+{
+  //
+  // Linear interpolation of data point with a parabolic expresion using
+  // the logarithm of the bin content in a close neighbourhood. It is
+  // assumed that the bin content of h is in number of events !
+  //
+  // Observation:
+  // This function have to be used when the statistics of each bin is
+  // sufficient and all bins are populated. For cases of sparse data
+  // please refere to Estimate2D1().
+  //
+  // Author : Alex Bercuci (A.Bercuci@gsi.de)
+  //
+
+       if(!h){
+               AliError("No histogram defined.");
+               return 0.;
+       }
+
+       TAxis *ax = h->GetXaxis(), *ay = h->GetYaxis();
+       Int_t binx   = ax->FindBin(x);
+       Int_t biny   = ay->FindBin(y);
+       Int_t nbinsx = ax->GetNbins();
+       Int_t nbinsy = ay->GetNbins();
+       Double_t p[2];
+       Double_t entries;
+               
+       fFitter2D2->ClearPoints();
+       Int_t npoints=0;
+       Int_t binx0, binx1, biny0, biny1;
+       for(int bin=0; bin<5; bin++){
+               binx0 = TMath::Max(1, binx-bin);
+               binx1 = TMath::Min(nbinsx, binx+bin);
+               for(int ibin=binx0; ibin<=binx1; ibin++){
+                       biny0 = TMath::Max(1, biny-bin);
+                       biny1 = TMath::Min(nbinsy, biny+bin);
+                       for(int jbin=biny0; jbin<=biny1; jbin++){
+                               if(ibin != binx0 && ibin != binx1 && jbin != biny0 && jbin != biny1) continue;
+                               if((entries = h->GetBinContent(ibin, jbin)) == 0.) continue;
+                               p[0] = ax->GetBinCenter(ibin);
+                               p[1] = ay->GetBinCenter(jbin);
+                               fFitter2D2->AddPoint(p, log(entries), 1./sqrt(entries));
+                               npoints++;
+                       }
+               }
+               if(npoints>=25) break;
+       }
+       if(fFitter2D2->Eval() == 1){
+               printf("<I2> x = %9.4f y = %9.4f\n", x, y);
+               printf("\tbinx %d biny %d\n", binx, biny);
+               printf("\tpoints %d\n", npoints);
+
+               return 0.;
+       }
+       TVectorD vec(6);
+       fFitter2D2->GetParameters(vec);
+       Double_t result = vec[0] + x*vec[1] + y*vec[2] + x*x*vec[3] + y*y*vec[4] + x*y*vec[5];
+       return exp(result);
+
+}
+
+//__________________________________________________________________
+Double_t AliTRDCalPIDLQRef::Estimate2D1(TH2 *h, Float_t &x, Float_t &y
+                                      , Float_t &dCT, Float_t &rmin
+                                      , Float_t &rmax)
+{
+  //
+  // Linear interpolation of data point with a plane using
+  // the logarithm of the bin content in the area defined by the
+  // d(cos(phi)) and dr=(rmin, rmax). It is assumed that the bin content
+  // of h is number of events !
+  //
+
+       if(!h){
+               AliError("No histogram defined.");
+               return 0.;
+       }
+
+       TAxis *ax = h->GetXaxis(), *ay = h->GetYaxis();
+//     Int_t binx   = ax->FindBin(x);
+//     Int_t biny   = ay->FindBin(y);
+       Int_t nbinsx = ax->GetNbins();
+       Int_t nbinsy = ay->GetNbins();
+       Double_t p[2];
+       Double_t entries;
+       Double_t rxy = sqrt(x*x + y*y), rpxy;
+               
+       fFitter2D1->ClearPoints();
+       Int_t npoints=0;
+       for(int ibin=1; ibin<=nbinsx; ibin++){
+               for(int jbin=1; jbin<=nbinsy; jbin++){
+                       if((entries = h->GetBinContent(ibin, jbin)) == 0.) continue;
+                       p[0] = ax->GetBinCenter(ibin);
+                       p[1] = ay->GetBinCenter(jbin);
+                       rpxy = sqrt(p[0]*p[0] + p[1]*p[1]);
+                       if((x*p[0] + y*p[1])/rxy/rpxy < dCT) continue;
+                       if(rpxy<rmin || rpxy > rmax) continue;
+                       
+                       fFitter2D1->AddPoint(p, log(entries), 1./sqrt(entries));
+                       npoints++;
+               }
+       }
+       if(npoints<15) return 0.;
+       if(fFitter2D1->Eval() == 1){
+               printf("<O2> x = %9.4f y = %9.4f\n", x, y);
+               printf("\tpoints %d\n", npoints);
+               return 0.;
+       }
+       TVectorD vec(3);
+       fFitter2D1->GetParameters(vec);
+       Double_t result = vec[0] + x*vec[1] + y*vec[2];
+       return exp(result);
+}
+
+//__________________________________________________________________
+// Double_t    AliTRDCalPIDLQRef::Estimate3D2(TH3 *h, Float_t &x, Float_t &y, Float_t &z)
+// {
+//     // Author Alex Bercuci (A.Bercuci@gsi.de)
+//     return 0.;
+// }
+
+//__________________________________________________________________
+void  AliTRDCalPIDLQRef::Prepare2DReferences()
+{
+  //
+  // Prepares the 2-dimensional reference histograms
+  //
+       
+       // histogram settings
+       Float_t xmin, xmax, ymin, ymax;
+       Int_t nbinsx, nbinsy;
+       const Int_t color[] = {4, 3, 2, 7, 6};
+       for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
+               // build reference histograms
+               nbinsx = nbinsy = 500;
+               xmin = ymin = 0.;
+               xmax = ymax = 500.;
+               if(!h2dEdx[ispec]){
+                       h2dEdx[ispec] = new  TH2D(Form("h2%s", AliTRDCalPIDLQ::fpartSymb[ispec]), "", nbinsx, xmin, xmax, nbinsy, ymin, ymax);
+                       h2dEdx[ispec]->SetLineColor(color[ispec]);
+               }
+               // build PCA
+               if(fPrinc[ispec]) fPrinc[ispec] = new TPrincipal(2, "ND");
+       }
+
+       // build transformed rotated histograms
+       nbinsx = nbinsy = 100;
+       xmin = ymin = -6.;
+       xmax = 10.; ymax = 6.;
+       TH2I *hProj = 0x0;
+       if(!(hProj = (TH2I*)gDirectory->Get("hProj"))) hProj = new  TH2I("hProj", "", nbinsx, xmin, xmax, nbinsy, ymin, ymax);
+       else hProj->Reset();
+       
+       // build transformed smoothed histogram
+       nbinsx = nbinsy = 200;
+       xmin = -1.5; xmax = 7.5;
+       ymin = -3.; ymax = 7.;
+       TH2D *hSmooth = 0x0;
+       if(!(hSmooth = (TH2D*)gDirectory->Get("hSmooth"))) hSmooth = new TH2D("hSmooth", "", nbinsx, xmin, xmax, nbinsy, ymin, ymax);
+       else hSmooth->Reset();
+       
+
+       printf("Doing interpolation and invertion ... "); fflush(stdout);
+       Bool_t kDebugPlot = kTRUE;
+       //Bool_t kDebugPrint = kFALSE;
+
+       TCanvas *c = 0x0;
+       TEllipse *ellipse = 0x0;
+       TMarker *mark = 0x0;
+       if(kDebugPlot){
+               c=new TCanvas("c2", "Interpolation 2D", 10, 10, 500, 500);
+               ellipse = new TEllipse();
+               ellipse->SetFillStyle(0); ellipse->SetLineColor(2);
+               mark = new TMarker();
+               mark->SetMarkerColor(2); mark->SetMarkerSize(2); mark->SetMarkerStyle(2);
+       }
+       
+       TAxis *ax, *ay;
+       Double_t xy[2], lxy[2];
+       Double_t estimate, position;
+       const TVectorD *eValues;
+       Float_t x0, y0, rx, ry, rc, rmin, rmax, dr, dCT;
+       for(int ispec=0; ispec<5; ispec++){
+               hProj->Reset(); hSmooth->Reset();
+               // calculate covariance ellipse
+               fPrinc[ispec]->MakePrincipals();
+               eValues  = fPrinc[ispec]->GetEigenValues();
+               x0  = 0.;
+               y0  = 0.;
+               rx  = 3.5*sqrt((*eValues)[0]);
+               ry  = 3.5*sqrt((*eValues)[1]);
+
+               // rotate to principal axis
+               Int_t irow = 0;
+               const Double_t *xx;
+               while((xx = fPrinc[ispec]->GetRow(irow++))){
+                       fPrinc[ispec]->X2P(xx, lxy);
+                       hProj->Fill(lxy[0], lxy[1]);
+               }
+
+               // debug plot
+               if(kDebugPlot){
+                       hProj->Draw();
+                       ellipse->DrawEllipse(x0, y0, rx, ry, 0., 360., 0.);
+                       mark->DrawMarker(x0, y0);
+                       gPad->Modified(); gPad->Update();
+               }
+               
+                               
+               // do interpolation
+               ax=hSmooth->GetXaxis();
+               ay=hSmooth->GetYaxis();
+               for(int ibin=1; ibin<=ax->GetNbins(); ibin++){
+                       xy[0] = ax->GetBinCenter(ibin); 
+                       for(int jbin=1; jbin<=ay->GetNbins(); jbin++){
+                               xy[1] = ay->GetBinCenter(jbin);
+
+                               // rotate to PCA
+                               fPrinc[ispec]->X2P(xy, lxy);
+
+                               // calculate border of covariance ellipse
+                               position = lxy[0]*lxy[0]/rx/rx+lxy[1]*lxy[1]/ry/ry;
+
+                               // interpolation inside the covariance ellipse
+                               if(position < 1.) estimate = Estimate2D2(hProj, (Float_t&)lxy[0], (Float_t&)lxy[1]);
+                               else{ // interpolation outside the covariance ellipse
+                                       dCT  = .9977;
+                                       dr = .5;
+                                       rc = sqrt(lxy[0]*lxy[0]+lxy[1]*lxy[1]);
+                                       rmin = rc - dr*sqrt(rc);
+                                       rmax = rc + dr/rc;
+                                       while((estimate = Estimate2D1(hProj, (Float_t&)lxy[0], (Float_t&)lxy[1], dCT, rmin, rmax)) < 1.E-8){
+                                               rmin -= dr;
+                                               dCT -= .001; 
+                                       }
+                               }
+                               hSmooth->SetBinContent(ibin, jbin, estimate);
+                       }
+               }
+               
+               // debug plot
+               if(kDebugPlot){
+                       hSmooth->Draw("lego2 fb"); gPad->SetLogz();
+                       gPad->Modified(); gPad->Update();
+               }
+                               
+               // doing invertion
+               ax=h2dEdx[ispec]->GetXaxis();
+               ay=h2dEdx[ispec]->GetYaxis();
+               for(int ibin=1; ibin<=ax->GetNbins(); ibin++){
+                       xy[0] = ax->GetBinCenter(ibin); 
+                       for(int jbin=1; jbin<=ay->GetNbins(); jbin++){
+                               xy[1] = ay->GetBinCenter(jbin);
+                               estimate = (ispec ==0 && (ibin==1 || jbin==1)) ? 1.E-10 : hSmooth->GetBinContent(hSmooth->FindBin(log(xy[0]), log(xy[1])))/xy[0]/xy[1];
+                               h2dEdx[ispec]->SetBinContent(ibin, jbin, (estimate>0.) ? estimate : 1.E-10);
+                       }
+               }
+               h2dEdx[ispec]->Scale(1./h2dEdx[ispec]->Integral());
+       }
+       printf("Done \n");
+}
+
+//__________________________________________________________________
+void   AliTRDCalPIDLQRef::Reset()
+{
+  //
+  // Resets the reference histograms
+  //
+
+       for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
+               h2dEdx[ispec]->Reset();
+               fPrinc[ispec]->Clear();
+       }       
+}
+
+//__________________________________________________________________
+void  AliTRDCalPIDLQRef::SaveReferences(const Int_t mom, const char *fn)
+{
+  //
+  // Save the reference histograms
+  //
+
+       TFile *fSave = 0x0;
+       TListIter it((TList*)gROOT->GetListOfFiles());
+       Bool_t kFOUND = kFALSE;
+       TDirectory *pwd = gDirectory;
+       while((fSave=(TFile*)it.Next()))
+               if(strcmp(fn, fSave->GetName())==0){
+                       kFOUND = kTRUE;
+                       break;
+               }
+       if(!kFOUND) fSave = new TFile(fn, "RECREATE");
+       fSave->cd();
+
+       TH2 *h;
+       for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
+               h = (TH2D*)h2dEdx[ispec]->Clone(Form("h2dEdx%s%d", AliTRDCalPIDLQ::fpartSymb[ispec], mom));
+               h->SetTitle(Form("2D dEdx for particle %s @ %d", AliTRDCalPIDLQ::fpartName[ispec], mom));
+               h->GetXaxis()->SetTitle("dE/dx_{TRD}^{amplif} [au]");
+               h->GetYaxis()->SetTitle("dE/dx_{TRD}^{drift} [au]");
+               h->GetZaxis()->SetTitle("Entries");
+               h->Write();
+       }
+               
+       pwd->cd();
+}
diff --git a/TRD/Cal/AliTRDCalPIDLQRef.h b/TRD/Cal/AliTRDCalPIDLQRef.h
new file mode 100644 (file)
index 0000000..02c47f5
--- /dev/null
@@ -0,0 +1,53 @@
+#ifndef ALITRDCALPIDLQREF_H
+#define ALITRDCALPIDLQREF_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  TRD calibration class for 2-dim reference histograms                     //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef ROOT_TObject
+#include "TObject.h"
+#endif
+
+class TH2;
+class TH3;
+class TPrincipal;
+class TLinearFitter;
+
+class AliTRDCalPIDLQRef : public TObject {
+
+ public:
+
+  AliTRDCalPIDLQRef();
+  AliTRDCalPIDLQRef(const AliTRDCalPIDLQRef &ref);
+  ~AliTRDCalPIDLQRef();
+  AliTRDCalPIDLQRef& operator=(const AliTRDCalPIDLQRef &ref);
+
+  TPrincipal*  GetPrincipal(const Int_t spec){return (spec>=0 && spec<5) ? fPrinc[spec] : 0x0;}
+  Double_t     Estimate2D2(TH2 *h, Float_t &x, Float_t &y);
+  Double_t     Estimate2D1(TH2 *h, Float_t &x, Float_t &y, Float_t &dCT, Float_t &rmin, Float_t &rmax);
+  //Double_t     Estimate3D2(TH3 *h, Float_t &x, Float_t &y, Float_t &z);
+  void         Prepare2DReferences();
+  void         Reset();
+  void         SaveReferences(const Int_t mom, const char *fn);
+       
+ private:
+
+  TPrincipal   *fPrinc[5];      // Used for principal component analysis
+  TLinearFitter        *fFitter2D2;     // Working object for linear fitter
+  TLinearFitter        *fFitter2D1;     // Working object for linear fitter
+
+  TH2          *h2dEdx[5];      // Data holders
+
+  ClassDef(AliTRDCalPIDLQRef, 1) // Reference histograms builder
+
+};
+
+#endif
+
diff --git a/TRD/Calib/PIDLQ/Run0_0_v0_s4.root b/TRD/Calib/PIDLQ/Run0_0_v0_s4.root
new file mode 100644 (file)
index 0000000..d6022f4
Binary files /dev/null and b/TRD/Calib/PIDLQ/Run0_0_v0_s4.root differ
index b6ed649..1fa71f7 100644 (file)
@@ -35,6 +35,7 @@
 #pragma link C++ class  AliTRDCalDet+;
 #pragma link C++ class  AliTRDCalGlobals+;
 #pragma link C++ class  AliTRDCalPIDLQ+;
+#pragma link C++ class  AliTRDCalPIDLQRef+;
 #pragma link C++ class  AliTRDCalMonitoring+;
 
 #pragma link C++ class  AliTRDCalChamberStatus+;
index bcfa0a1..0750863 100644 (file)
@@ -19,6 +19,7 @@ SRCS= AliTRDarrayI.cxx \
       Cal/AliTRDCalDet.cxx \
       Cal/AliTRDCalGlobals.cxx \
       Cal/AliTRDCalPIDLQ.cxx \
+      Cal/AliTRDCalPIDLQRef.cxx \
       Cal/AliTRDCalMonitoring.cxx \
       Cal/AliTRDCalChamberStatus.cxx \
       Cal/AliTRDCalPadStatus.cxx \