First implementation of neural network PID
authorcblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 26 Sep 2007 14:01:56 +0000 (14:01 +0000)
committercblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 26 Sep 2007 14:01:56 +0000 (14:01 +0000)
21 files changed:
TRD/AliTRDReconstructor.cxx
TRD/AliTRDReconstructor.h
TRD/AliTRDcalibDB.cxx
TRD/AliTRDcalibDB.h
TRD/AliTRDgtuTrack.cxx
TRD/AliTRDpidESD.cxx
TRD/AliTRDtrack.cxx
TRD/AliTRDtrack.h
TRD/AliTRDtracker.cxx
TRD/Cal/AliTRDCalPID.cxx
TRD/Cal/AliTRDCalPID.h
TRD/Cal/AliTRDCalPIDLQ.cxx [new file with mode: 0644]
TRD/Cal/AliTRDCalPIDLQ.h [new file with mode: 0644]
TRD/Cal/AliTRDCalPIDNN.cxx [new file with mode: 0644]
TRD/Cal/AliTRDCalPIDNN.h [new file with mode: 0644]
TRD/Cal/AliTRDCalPIDRefMaker.cxx
TRD/Calib/PIDLQ/Run0_999999999_v0_s0.root
TRD/Calib/PIDNN/Run0_999999999_v0_s0.root [new file with mode: 0644]
TRD/Macros/AliTRDpidCDB.C
TRD/TRDbaseLinkDef.h
TRD/libTRDbase.pkg

index a57f0a2..df406e5 100644 (file)
@@ -63,42 +63,6 @@ void AliTRDReconstructor::ConvertDigits(AliRawReader *rawReader
 }
 
 //_____________________________________________________________________________
-void AliTRDReconstructor::Reconstruct(AliRunLoader *runLoader
-                                    , AliRawReader *rawReader) const
-{
-  //
-  // Reconstruct clusters
-  //
-
-  AliInfo("Reconstruct TRD clusters from RAW data [RunLoader, RawReader]");
-
-  AliLoader *loader = runLoader->GetLoader("TRDLoader");
-  loader->LoadRecPoints("recreate");
-
-  runLoader->CdGAFile();
-  Int_t nEvents = runLoader->GetNumberOfEvents();
-
-  rawReader->Reset();
-  rawReader->Select("TRD");
-
-  for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
-
-    if (!rawReader->NextEvent()) break;
-
-    // New (fast) cluster finder
-    AliTRDclusterizer clusterer("clusterer","TRD clusterizer");
-    clusterer.Open(runLoader->GetFileName(),iEvent);
-    clusterer.Raw2ClustersChamber(rawReader);
-
-    clusterer.WriteClusters(-1);
-
-  }
-
-  loader->UnloadRecPoints();
-
-}
-
-//_____________________________________________________________________________
 void AliTRDReconstructor::Reconstruct(AliRawReader *rawReader
                                     , TTree *clusterTree) const
 {
@@ -136,63 +100,13 @@ void AliTRDReconstructor::Reconstruct(TTree *digitsTree
 }
 
 //_____________________________________________________________________________
-void AliTRDReconstructor::Reconstruct(AliRunLoader *runLoader) const
-{
-  //
-  // Reconstruct clusters
-  //
-
-  AliInfo("Reconstruct TRD clusters [AliRunLoader]");
-  AliLoader *loader = runLoader->GetLoader("TRDLoader");
-  loader->LoadRecPoints("recreate");
-
-  runLoader->CdGAFile();
-  Int_t nEvents = runLoader->GetNumberOfEvents();
-
-  for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
-    AliTRDclusterizer clusterer("clusterer","TRD clusterizer");
-    clusterer.Open(runLoader->GetFileName(),iEvent);
-    clusterer.ReadDigits();
-    clusterer.MakeClusters();
-    clusterer.WriteClusters(-1);
-  }
-
-  loader->UnloadRecPoints();
-
-}
-
-//_____________________________________________________________________________
-AliTracker *AliTRDReconstructor::CreateTracker(AliRunLoader *runLoader) const
+AliTracker *AliTRDReconstructor::CreateTracker() const
 {
   //
   // Create a TRD tracker
   //
 
-  runLoader->CdGAFile();
-
-  return new AliTRDtracker(gFile);
-
-}
-
-//_____________________________________________________________________________
-void AliTRDReconstructor::FillESD(AliRunLoader* /*runLoader*/
-                               , AliRawReader* /*rawReader*/
-                               , AliESDEvent* /*esd*/) const
-{
-  //
-  // Fill ESD
-  //
-
-}
-
-//_____________________________________________________________________________
-void AliTRDReconstructor::FillESD(AliRawReader* /*rawReader*/
-                               , TTree* /*clusterTree*/
-                               , AliESDEvent* /*esd*/) const
-{
-  //
-  // Fill ESD
-  //
+  return new AliTRDtracker(NULL);
 
 }
 
@@ -206,13 +120,3 @@ void AliTRDReconstructor::FillESD(TTree* /*digitsTree*/
   //
 
 }
-
-//_____________________________________________________________________________
-void AliTRDReconstructor::FillESD(AliRunLoader* /*runLoader*/
-                               , AliESDEvent* /*esd*/) const
-{
-  //
-  // Fill ESD
-  //
-
-}
index 624746f..517d987 100644 (file)
@@ -26,18 +26,14 @@ class AliTRDReconstructor: public AliReconstructor {
   virtual Bool_t      HasDigitConversion() const                 { return kFALSE; };
   virtual void        ConvertDigits(AliRawReader *rawReader, TTree *digitsTree) const;
 
-  virtual Bool_t      HasLocalReconstruction() const             { return kTRUE; };
-  virtual void        Reconstruct(AliRunLoader *runLoader, AliRawReader *rawReader) const;
   virtual void        Reconstruct(AliRawReader *rawReader, TTree *clusterTree) const;
   virtual void        Reconstruct(TTree *digitsTree, TTree *clusterTree) const;
-  virtual void        Reconstruct(AliRunLoader *runLoader) const;
 
-  virtual AliTracker *CreateTracker(AliRunLoader *runLoader) const;
+  virtual AliTracker *CreateTracker() const;
 
-  virtual void        FillESD(AliRunLoader *runLoader, AliRawReader *rawReader, AliESDEvent *esd) const;
-  virtual void        FillESD(AliRawReader *rawReader, TTree *clusterTree, AliESDEvent *esd) const;
+  virtual void        FillESD(AliRawReader */*rawReader*/, TTree *clusterTree, AliESDEvent *esd) const
+  {FillESD((TTree*)NULL,clusterTree,esd);}
   virtual void        FillESD(TTree *digitsTree, TTree *clusterTree, AliESDEvent *esd) const;
-  virtual void        FillESD(AliRunLoader *runLoader, AliESDEvent *esd) const;
 
   static  void        SetSeedingOn(Bool_t seeding)               { fgkSeedingOn  = seeding; }  
   static  void        SetStreamLevel(Int_t level)                { fgStreamLevel = level;   }
index 2f3e1f0..f739204 100644 (file)
@@ -244,6 +244,8 @@ const TObject *AliTRDcalibDB::GetCachedCDBObject(Int_t id)
     case kIDFEE : 
       return CacheCDBEntry(kIDFEE               ,"TRD/Calib/FEE"); 
       break;
+    case kIDPIDNN : 
+      return CacheCDBEntry(kIDPIDNN             ,"TRD/Calib/PIDNN");
     case kIDPIDLQ : 
       return CacheCDBEntry(kIDPIDLQ             ,"TRD/Calib/PIDLQ"); 
       break;
@@ -798,14 +800,17 @@ Bool_t AliTRDcalibDB::IsChamberMasked(Int_t det)
 }
 
 //_____________________________________________________________________________
-const AliTRDCalPID *AliTRDcalibDB::GetPIDLQObject()
+const AliTRDCalPID *AliTRDcalibDB::GetPIDObject(const Int_t method)
 {
   //
   // Returns the object storing the distributions for PID with likelihood
   //
-
-  return dynamic_cast<const AliTRDCalPID *> 
-         (GetCachedCDBObject(kIDPIDLQ));
+       switch(method){
+       case 0: return dynamic_cast<const AliTRDCalPID *> 
+                 (GetCachedCDBObject(kIDPIDNN));
+       case 1: return dynamic_cast<const AliTRDCalPID *>
+                 (GetCachedCDBObject(kIDPIDLQ));
+       }
 
 }
 
index a153728..d91e82a 100644 (file)
@@ -66,7 +66,7 @@ class AliTRDcalibDB : public TObject {
   Bool_t                     IsChamberMasked(Int_t det);
 
   const AliTRDCalMonitoring *GetMonitoringObject();
-  const AliTRDCalPID        *GetPIDLQObject();
+  const AliTRDCalPID        *GetPIDObject(const Int_t method);
 
   // Related functions, these depend on calibration data
   static Float_t             GetOmegaTau(Float_t vdrift, Float_t bz);
@@ -76,7 +76,7 @@ class AliTRDcalibDB : public TObject {
  protected:
 
   // For caching see also implentation of GetCachedCDBObject in the .cxx file
-  enum { kCDBCacheSize = 15 };   // Number of cached objects
+  enum { kCDBCacheSize = 16 };   // Number of cached objects
   enum { kIDVdriftPad = 0
        , kIDVdriftChamber
        , kIDT0Pad
@@ -88,6 +88,7 @@ class AliTRDcalibDB : public TObject {
        , kIDChamberPos
        , kIDStackPos
        , kIDSuperModulePos
+       , kIDPIDNN
        , kIDPIDLQ
        , kIDMonitoringData
        , kIDChamberStatus
index abbb809..3b7335b 100644 (file)
@@ -390,7 +390,7 @@ void AliTRDgtuTrack::MakePID()
     AliError("No instance of AliTRDcalibDB.");
     return;  
   }
-  const AliTRDCalPID *pd = calibration->GetPIDLQObject();
+  const AliTRDCalPID *pd = calibration->GetPIDObject(1);
   
   AliTRDltuTracklet *trk;
   Int_t   nTracklets = GetNtracklets();
@@ -455,8 +455,8 @@ void AliTRDgtuTrack::MakePID()
                 dedx[0] = dedx[1] = q*3.; dedx[2] = 0.;
                 Float_t length = 3.7;
 
-    probEle *= pd->GetProbability(0, TMath::Abs(fPt), dedx, length);
-    probPio *= pd->GetProbability(2, TMath::Abs(fPt), dedx, length);
+    probEle *= pd->GetProbability(0, TMath::Abs(fPt), dedx, length, 0);
+    probPio *= pd->GetProbability(2, TMath::Abs(fPt), dedx, length, 0);
 
   }
 
index 87b15a7..2ffc94a 100644 (file)
@@ -130,7 +130,7 @@ Int_t AliTRDpidESD::MakePID(AliESDEvent *event)
        }
        
        // Retrieve the CDB container class with the probability distributions
-       const AliTRDCalPID *pd = calibration->GetPIDLQObject();
+       const AliTRDCalPID *pd = calibration->GetPIDObject(1);
        if (!pd) {
                AliErrorGeneral("AliTRDpidESD::MakePID()"
                        ,"No access to AliTRDCalPID");
@@ -180,7 +180,7 @@ Int_t AliTRDpidESD::MakePID(AliESDEvent *event)
 
                        // 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->GetProbability(iSpecies, mom, dedx, length, iPlan);
                                //p[iSpecies] *= pd->GetProbabilityT(iSpecies, mom, timebin);
                        }
                }
index 6a8af33..f2c518e 100644 (file)
@@ -469,7 +469,7 @@ void AliTRDtrack::CookdEdx(Double_t low, Double_t up)
 }                     
 
 //_____________________________________________________________________________
-void AliTRDtrack::CookdEdxTimBin()
+void AliTRDtrack::CookdEdxTimBin(const Int_t/* tid*/)
 {
   //
   // Set fdEdxPlane and fTimBinPlane and also get the 
@@ -544,6 +544,53 @@ void AliTRDtrack::CookdEdxTimBin()
 }
 
 //_____________________________________________________________________________
+void   AliTRDtrack::CookdEdxNN(Float_t *dedx){
+
+  // This function calcuates the input for the neural networks 
+  // which are used for the PID. 
+
+       Int_t ntb = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();//number of time bins in chamber
+       Int_t plane;   // plane of current cluster
+       Int_t tb;      // time bin of current cluster
+       Int_t slice;   // curent slice
+       AliTRDcluster *cluster = 0x0; // pointer to current cluster
+       const Int_t lMLPscale = 16000; // scaling of the MLP input to be smaller than 1
+
+       // Reset class and local contors/variables
+       for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++){
+         for (Int_t iSlice = 0; iSlice < kNMLPslice; iSlice++) {
+           *(dedx + (kNMLPslice * iPlane) + iSlice) = 0.;
+         }
+       }
+
+       // start looping over clusters attached to this track
+       for (Int_t iClus = 0; iClus < GetNumberOfClusters(); iClus++) {
+         cluster = fClusters[iClus]; //(AliTRDcluster*)tracker->GetCluster(fIndex[iClus]);
+         if(!cluster) continue;
+         
+         // Read info from current cluster
+         plane  = AliTRDgeometry::GetPlane(cluster->GetDetector());
+         if (plane < 0 || plane >= AliESDtrack::kNPlane) {
+           AliError(Form("Wrong plane %d", plane));
+           continue;
+         }
+
+         tb     = cluster->GetLocalTimeBin();
+         if(tb == 0 || tb >= ntb){
+           AliWarning(Form("time bin 0 or > %d in cluster %d", ntb, iClus));
+           AliInfo(Form("dQ/dl %f", fdQdl[iClus]));
+           continue;
+         }
+
+         slice   = tb * kNMLPslice / ntb;
+         
+         *(dedx+(kNMLPslice * plane) + slice) += fdQdl[iClus]/lMLPscale;
+       } // End of loop over cluster
+        
+}
+
+
+//_____________________________________________________________________________
 void   AliTRDtrack::SetTrackSegmentDirMom(const Int_t plane)
 {
   //
@@ -580,7 +627,7 @@ Float_t     AliTRDtrack::GetTrackLengthPlane(Int_t plane) const
 }
 
 //_____________________________________________________________________________
-Int_t AliTRDtrack::CookPID(AliESDtrack *esd)
+Bool_t AliTRDtrack::CookPID(Int_t &pidQuality)
 {
   //
   // This function calculates the PID probabilities based on TRD signals
@@ -598,65 +645,80 @@ Int_t AliTRDtrack::CookPID(AliESDtrack *esd)
   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
   if (!calibration) {
     AliError("No access to calibration data");
-    return -1;
+    return kFALSE;
   }
        
-       // Retrieve the CDB container class with the probability distributions
-       const AliTRDCalPID *pd = calibration->GetPIDLQObject();
-       if (!pd) {
-               AliError("No access to AliTRDCalPID");
-               return -1;
-       }
+  // Retrieve the CDB container class with the probability distributions
+  const AliTRDCalPID *pd = calibration->GetPIDObject(fPIDmethod == kNN ? 0 : 1);
+  if (!pd) {
+    AliError("No access to AliTRDCalPID");
+    return kFALSE;
+  }
 
+  // Calculate the input for the NN if fPIDmethod is kNN
+  Float_t ldEdxNN[AliTRDCalPID::kNPlane * kNMLPslice], *dedx = 0x0;
+  if(fPIDmethod == kNN) {
+    CookdEdxNN(&ldEdxNN[0]);
+  }
 
-       Double_t p0 = 1./AliPID::kSPECIES;
-       if(AliPID::kSPECIES != 5){
-               AliError("Probabilities array defined only for 5 values. Please modify !!");
-               return -1;
-       }
-       Double_t p[] = {p0, p0, p0, p0, p0};
-       Float_t length;  // track segment length in chamber
-       Int_t nPlanePID = 0;
-       // Skip tracks which have no TRD signal at all
-       if (fdEdx == 0.) return -1;
-       
-       for (Int_t iPlane = 0; iPlane < AliTRDgeometry::kNplan; iPlane++) {
+  // Sets the a priori probabilities
+  for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
+    fPID[ispec] = 1./AliPID::kSPECIES; 
+  }
 
-               length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())/TMath::Sqrt((1. - fSnp[iPlane]*fSnp[iPlane]) / (1. + fTgl[iPlane]*fTgl[iPlane]));
 
-               // check data
-               if((fdEdxPlane[iPlane][0] + fdEdxPlane[iPlane][1] + fdEdxPlane[iPlane][2]) <=  0.
-               || fTimBinPlane[iPlane] == -1.) continue;
+  if(AliPID::kSPECIES != 5){
+    AliError("Probabilities array defined only for 5 values. Please modify !!");
+    return kFALSE;
+  }
 
-               // 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, fMom[iPlane], fdEdxPlane[iPlane], length);
-                       //p[iSpecies] *= pd->GetProbabilityT(iSpecies, fMom[iPlane], fTimBinPlane[iPlane]);
-               }
-       }
-       if(nPlanePID == 0) return 0;
+  pidQuality = 0;
+  Float_t length;  // track segment length in chamber
 
-       // normalize probabilities
-       Double_t probTotal = 0.;
-       for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) probTotal += p[iSpecies];
+  // Skip tracks which have no TRD signal at all
+  if (fdEdx == 0.) return kFALSE;
+       
+  for (Int_t iPlane = 0; iPlane < AliTRDgeometry::kNplan; iPlane++) {
+
+    length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())/TMath::Sqrt((1. - fSnp[iPlane]*fSnp[iPlane]) / (1. + fTgl[iPlane]*fTgl[iPlane]));
+
+    // check data
+    if((fdEdxPlane[iPlane][0] + fdEdxPlane[iPlane][1] + fdEdxPlane[iPlane][2]) <=  0.
+       || fTimBinPlane[iPlane] == -1.) continue;
+
+    // this track segment has fulfilled all requierments
+    pidQuality++;
+
+    // Get the probabilities for the different particle species
+    for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
+      switch(fPIDmethod){
+      case kLQ:
+       dedx = fdEdxPlane[iPlane];
+       break;
+      case kNN:
+       dedx = &ldEdxNN[iPlane*kNMLPslice];
+       break;
+      }
+      fPID[iSpecies] *= pd->GetProbability(iSpecies, fMom[iPlane], dedx, length, iPlane);
+    }
+  }
+  if(pidQuality == 0) return kTRUE;
 
-       if(probTotal <= 0.) {
-               AliWarning("The total probability over all species <= 0. This may be caused by some error in the reference histograms.");
-               return -1;
-       }
-       for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) p[iSpecies] /= probTotal;
+  // normalize probabilities
+  Double_t probTotal = 0.;
+  for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) probTotal += fPID[iSpecies];
 
+  if(probTotal <= 0.) {
+    AliWarning("The total probability over all species <= 0. This may be caused by some error in the reference histograms.");
+    return kFALSE;
+  }
+  for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) fPID[iSpecies] /= probTotal;
        
-       // book PID to the ESD track
-       esd->SetTRDpid(p);
-       esd->SetTRDpidQuality(nPlanePID);
-       
-       return 0;
+  return kTRUE;
 }
 
+
 //_____________________________________________________________________________
 Bool_t AliTRDtrack::PropagateTo(Double_t xk, Double_t xx0, Double_t xrho)
 {
@@ -784,7 +846,7 @@ Bool_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, Int_t index
 
 //_____________________________________________________________________________
 Int_t AliTRDtrack::UpdateMI(AliTRDcluster *c, Double_t chisq, Int_t index
-                          , Double_t h01, Int_t /*plane*/)
+                          , Double_t h01, Int_t /*plane*/, Int_t tid)
 {
   //
   // Assignes found cluster to the track and updates track information
index 486757a..629d91c 100644 (file)
@@ -31,7 +31,13 @@ class AliTRDtrack : public AliKalmanTrack {
        , kNplane    =   6
        , kNcham     =   5
        , kNsect     =  18
-       , kNslice    =   3};
+       , kNslice    =   3
+       , kNMLPslice =   8};
+  
+  enum AliTRDPIDMethod {
+         kNN = 0
+       , kLQ = 1
+  };
 
    AliTRDtrack();
    AliTRDtrack(/*const */AliTRDcluster *c, Int_t index, const Double_t xx[5], const Double_t cc[15], Double_t xr, Double_t alpha);
@@ -55,6 +61,10 @@ class AliTRDtrack : public AliKalmanTrack {
            Int_t    GetClusterIndex(Int_t i) const                          { return fIndex[i];                    }
            Double_t GetC() const                                            { return AliExternalTrackParam::GetC(GetBz()); }
            Double_t GetPredictedChi2(const AliTRDcluster *c, Double_t h01) const;
+
+          inline AliTRDPIDMethod GetPIDMethod() const {return fPIDmethod;}
+          inline void      SetPIDMethod(AliTRDPIDMethod method) {fPIDmethod = method;}
+                                       // end
            Float_t  GetPIDsignals(Int_t iPlane, Int_t iSlice) const         { return fdEdxPlane[iPlane][iSlice];   }
            Int_t    GetPIDTimBin(Int_t i) const                             { return fTimBinPlane[i];              }
            Double_t GetLikelihoodElectron() const                           { return fLhElectron;                  }
@@ -94,11 +104,12 @@ class AliTRDtrack : public AliKalmanTrack {
 // A. Bercuci new functions
                                        void            SetTrackSegmentDirMom(const Int_t plane);
                                        void            CookdEdx(Double_t low = 0.05, Double_t up = 0.7);
-                                       void            CookdEdxTimBin();
-                                       Int_t           CookPID(AliESDtrack *p);
+                                       void            CookdEdxTimBin(const Int_t tid);
+                                       Bool_t          CookPID(Int_t &pidQuality); 
                                        void            SetCluster(AliTRDcluster* cl, Int_t index = -1) {fClusters[(index == -1) ? GetNumberOfClusters()-1 : index] = cl;}
        inline  AliTRDcluster* GetCluster(Int_t layer){ return (layer >= 0 && layer < GetNumberOfClusters()) ?  fClusters[layer] : 0x0;}
        inline  Float_t GetMomentumPlane(Int_t plane) const {return (plane >= 0 && plane < kNplane) ? fMom[plane] : 0.;}
+       inline  const Double_t* GetPID() const {return fPID;}
        inline  Float_t GetSnpPlane(Int_t plane) const {return (plane >= 0 && plane < kNplane) ? fSnp[plane] : 0.;}
        inline  Float_t GetTglPlane(Int_t plane) const {return (plane >= 0 && plane < kNplane) ? fTgl[plane] : 0.;}
                Float_t GetTrackLengthPlane(Int_t plane) const;
@@ -111,12 +122,13 @@ class AliTRDtrack : public AliKalmanTrack {
            Bool_t   Rotate(Double_t angle, Bool_t absolute = kFALSE);
            Float_t  StatusForTOF();
            Bool_t   Update(const AliTRDcluster *c, Double_t chi2, Int_t i, Double_t h01);
-           Int_t    UpdateMI(/*const */AliTRDcluster *c, Double_t chi2, Int_t i, Double_t h01, Int_t plane);
+           Int_t    UpdateMI(/*const */AliTRDcluster *c, Double_t chi2, Int_t i, Double_t h01, Int_t plane, Int_t tid=0);
  
  protected:
 
   AliTRDtrack      &operator=(const AliTRDtrack &t);
 
+          void     CookdEdxNN(Float_t *dedx);
            Double_t GetBz() const;
            Bool_t   Update(const AliCluster */*c*/, Double_t /*chi2*/, Int_t /*idx*/) { return 0;   }
            Double_t GetPredictedChi2(const AliCluster* /*c*/) const                   { return 0.0; }
@@ -126,6 +138,7 @@ class AliTRDtrack : public AliKalmanTrack {
            Float_t  fDE;                                //  Integrated delta energy
            Float_t  fdEdxPlane[kNplane][kNslice];       //  dE/dx from all 6 planes in 3 slices each
            Int_t    fTimBinPlane[kNplane];              //  Time bin of Max cluster from all 6 planes
+           Double_t fPID[AliPID::kSPECIES];             // PID probabilities
 
                                        // A.Bercuci
                                         Float_t  fMom[kNplane];                      //  Track momentum at chamber entrance
@@ -134,7 +147,7 @@ class AliTRDtrack : public AliKalmanTrack {
            AliTRDcluster*   fClusters[kMAXCLUSTERSPERTRACK];
            Bool_t   fClusterOwner;         // indicates the track is owner of cluster
            // end A.Bercuci
-                                        
+          AliTRDPIDMethod fPIDmethod;   // switch between different PID methods                                         
                                         Bool_t   fStopped;                           //  Track stop indication
            Int_t    fIndex[kMAXCLUSTERSPERTRACK];       //  Global indexes of clusters  
            Int_t    fIndexBackup[kMAXCLUSTERSPERTRACK]; //  Backup indexes of clusters - used in iterations
@@ -153,7 +166,7 @@ class AliTRDtrack : public AliKalmanTrack {
            Float_t  fBudget[3];                         //  Integrated material budget
    AliTRDtrack     *fBackupTrack;                       //! Backup track
 
-   ClassDef(AliTRDtrack,8)                              //  TRD reconstructed tracks
+   ClassDef(AliTRDtrack,9)                              //  TRD reconstructed tracks
 
 };                     
 
index f839b3f..636ec25 100644 (file)
@@ -567,7 +567,7 @@ Int_t AliTRDtracker::PropagateBack(AliESDEvent *event)
       Int_t foundClr = track->GetNumberOfClusters();
       if (foundClr >= foundMin) {
        track->CookdEdx(); 
-       track->CookdEdxTimBin(); // A.Bercuci 25.07.07
+       track->CookdEdxTimBin(seed->GetID()); // A.Bercuci 25.07.07
        CookLabel(track,1 - fgkLabelFraction);
        if (track->GetBackupTrack()) {
           UseClusters(track->GetBackupTrack());
@@ -779,6 +779,7 @@ Int_t AliTRDtracker::RefitInward(AliESDEvent *event)
   Float_t foundMin = fgkMinClustersInTrack * timeBins; 
   Int_t   nseed    = 0;
   Int_t   found    = 0;
+  Int_t   pidQ     = 0;
   //Int_t innerTB    = fTrSec[0]->GetInnerTimeBin();
   AliTRDtrack seed2;
   
@@ -830,14 +831,15 @@ Int_t AliTRDtracker::RefitInward(AliESDEvent *event)
     //AliTRDtrack *pt = seed2;
     //AliTRDtrack &t = *pt; 
     FollowProlongation(*pt); 
-    if (pt->GetNumberOfClusters() >= foundMin) {
+    //if (pt->GetNumberOfClusters() >= foundMin) {
       //UseClusters(&t);
       //CookLabel(pt, 1-fgkLabelFraction);
       pt->CookdEdx();
-      pt->CookdEdxTimBin();
-      pt->CookPID(seed);
+      pt->CookdEdxTimBin(seed->GetID());
+      pt->SetPIDMethod(AliTRDtrack::kLQ);  //switch between TRD PID methods
+      pt->CookPID(pidQ);
       //pt->Calibrate(); // slot for calibration
-    }
+      //}
     found++;
 
     Double_t xTPC = 250.0;
@@ -862,7 +864,7 @@ Int_t AliTRDtracker::RefitInward(AliESDEvent *event)
       delete seed2;
       if (PropagateToX(*pt2,xTPC,fgkMaxStep)) { 
         pt2->CookdEdx( ); 
-        pt2->CookdEdxTimBin(); 
+        pt2->CookdEdxTimBin(seed->GetID()); 
        seed->UpdateTrackParams(pt2,AliESDtrack::kTRDrefit);
        fHRefit->Fill(6);
 
@@ -3981,8 +3983,9 @@ Int_t AliTRDtracker::Freq(Int_t n, const Int_t *inlist
   // The size of output array has is 2*n 
   //
 
-  if (n <= 0)
+  if (n <= 0) {
     return 0;
+  }
 
   Int_t *sindexS = new Int_t[n];   // Temporary array for sorting
   Int_t *sindexF = new Int_t[2*n];   
@@ -4085,7 +4088,7 @@ AliTRDtrack *AliTRDtracker::RegisterSeed(AliTRDseed *seeds, Double_t *params)
   }
   else {
     track->CookdEdx();
-    track->CookdEdxTimBin();
+    track->CookdEdxTimBin(-1);
     CookLabel(track,0.9);
   }
 
index 31c1983..9469ac4 100644 (file)
@@ -17,8 +17,7 @@
 
 //////////////////////////////////////////////////////////////////////////
 //                                                                      //
-// Container for the distributions of dE/dx and the time bin of the     //
-// max. cluster for electrons and pions                                 //
+// Container for PID information                                        //
 //                                                                      //
 // Authors:                                                             //
 //   Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de>               //
 
 ClassImp(AliTRDCalPID)
 
-Char_t* AliTRDCalPID::fpartName[AliPID::kSPECIES] = {"electron", "muon", "pion", "kaon", "proton"};
+Char_t* AliTRDCalPID::fPartName[AliPID::kSPECIES] = { "electron", "muon", "pion", "kaon", "proton"};
+Char_t* AliTRDCalPID::fPartSymb[AliPID::kSPECIES] = { "EL", "MU", "PI", "KA", "PR"};
+Float_t AliTRDCalPID::fTrackMomentum[kNMom]       = {  0.6,  0.8,  1.0,  1.5,  2.0
+                                                   ,   3.0,  4.0,  5.0,  6.0,  8.0
+                                                   ,  10.0};
 
-Char_t* AliTRDCalPID::fpartSymb[AliPID::kSPECIES] = {"EL", "MU", "PI", "KA", "PR"};
-
-Float_t AliTRDCalPID::fTrackMomentum[kNMom] = {0.6,  0.8,  1.0,  1.5,  2.0,  3.0,  4.0,  5.0,  6.0,  8.0,  10.0};
-  
-Float_t AliTRDCalPID::fTrackSegLength[kNLength] = {3.7, 3.9, 4.2, 5.0};
-
-    
 //_________________________________________________________________________
 AliTRDCalPID::AliTRDCalPID()
   :TNamed("pid", "PID for TRD")
-  ,fMeanChargeRatio(0)
-  ,fHistdEdx(0x0)
-  ,fHistTimeBin(0x0)
+  ,fModel(0x0)
 {
   //
   //  The Default constructor
   //
 
-  Init();
-
 }
 
-//_________________________________________________________________________
+//_____________________________________________________________________________
 AliTRDCalPID::AliTRDCalPID(const Text_t *name, const Text_t *title) 
   :TNamed(name,title)
-  ,fMeanChargeRatio(0)
-  ,fHistdEdx(0x0)
-  ,fHistTimeBin(0x0)
+  ,fModel(0x0)
 {
   //
   //  The main constructor
-  //
-  
-  Init();
-
-}
-
-//_____________________________________________________________________________
-AliTRDCalPID::AliTRDCalPID(const AliTRDCalPID &c) 
-  :TNamed(c)
-  ,fMeanChargeRatio(c.fMeanChargeRatio)
-  ,fHistdEdx(0x0)
-  ,fHistTimeBin(0x0)
-{
-  //
-  // Copy constructor
-  //
+  //  
 
-  if (this != &c) ((AliTRDCalPID &) c).Copy(*this);
-  
 }
 
 //_________________________________________________________________________
@@ -102,312 +75,8 @@ AliTRDCalPID::~AliTRDCalPID()
   // Destructor
   //
   
-  CleanUp();
-
-}
-
-//_________________________________________________________________________
-void AliTRDCalPID::CleanUp()
-{
-  //
-  // Delets all newly created objects
-  //
-
-  if (fHistdEdx) {
-    delete fHistdEdx;
-    fHistdEdx = 0x0;
-  }
-  
-  if (fHistTimeBin) {
-    delete fHistTimeBin;
-    fHistTimeBin = 0x0;
-  }
-}
-
-//_____________________________________________________________________________
-AliTRDCalPID &AliTRDCalPID::operator=(const AliTRDCalPID &c)
-{
-  //
-  // Assignment operator
-  //
-
-  if (this != &c) ((AliTRDCalPID &) c).Copy(*this);
-  return *this;
-
-}
-
-//_____________________________________________________________________________
-void AliTRDCalPID::Copy(TObject &c) const
-{
-  //
-  // Copy function
-  //
-
-  AliTRDCalPID& target = (AliTRDCalPID &) c;
-  
-  target.CleanUp();
-  
-  target.fMeanChargeRatio = fMeanChargeRatio;
-
-  if (fHistdEdx) {
-    target.fHistdEdx = (TObjArray*) fHistdEdx->Clone();
-  }
-  if (fHistTimeBin) {
-    target.fHistTimeBin = (TObjArray*) fHistTimeBin->Clone();
-  }
-
-  TObject::Copy(c);
-
-}
-
-//_________________________________________________________________________
-void AliTRDCalPID::Init()
-{
-  //
-  // Initialization
-  //
-
-  fHistdEdx    = new TObjArray(AliPID::kSPECIES * kNMom/* * kNLength*/);
-  fHistdEdx->SetOwner();
-  fHistTimeBin = new TObjArray(2 * kNMom);
-  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 AliTRDCalPIDRefMaker();
-       
-  // ADC Gain normalization
-  fMeanChargeRatio = 1.0;
-}
-
-//_________________________________________________________________________
-Bool_t AliTRDCalPID::LoadLQReferences(Char_t *refFile)
-{
-  //
-  // Read the TRD dEdx histograms.
-  //
-
-       Int_t nTimeBins = 22;
-       // Get number of time bins from CDB
-       AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
-       if(!calibration){
-               AliWarning(Form("No AliTRDcalibDB available. Using %d time bins.", nTimeBins));
-       }else{
-               if(calibration->GetRun() > -1) nTimeBins = calibration->GetNumberOfTimeBins();
-               else AliWarning(Form("Run number not set. Using %d time bins.", nTimeBins));
-       }
-
-       
-  // Read histogram Root file  
-  TFile *histFile = TFile::Open(refFile, "READ");
-  if (!histFile || !histFile->IsOpen()) {
-    AliError(Form("Opening TRD histgram file %s failed", refFile));
-    return kFALSE;
+  if (fModel) {
+    delete fModel;
   }
-  gROOT->cd();
-
-  // Read histograms
-  for (Int_t iparticle = 0; iparticle < AliPID::kSPECIES; iparticle++){
-    for (Int_t imom = 0; imom < kNMom; 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() != nTimeBins) 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, nTimeBins));
-                       ht->Scale(1./ht->Integral());
-                       fHistTimeBin->AddAt(ht, ((iparticle==AliPID::kElectron)?0:1)*kNMom + 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);
-  
-  return kTRUE;
-
-}
-
-// //_________________________________________________________________________
-// Double_t  AliTRDCalPID::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  AliTRDCalPID::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();
-// 
-// }
-
-//_________________________________________________________________________
-TH1* AliTRDCalPID::GetHistogram(Int_t k, Int_t ip/*, Int_t il*/) const
-{
-  //
-  // Returns one selected dEdx histogram
-  //
-
-  if (k < 0 || k >= AliPID::kSPECIES) return 0x0;
-       if(ip<0 || ip>= kNMom ) return 0x0;
-
-       AliInfo(Form("Retrive dEdx histogram for %s of %5.2f GeV/c", fpartName[k], fTrackMomentum[ip]));
-  
-  return (TH1*)fHistdEdx->At(GetHistID(k, ip));
-
-}
-
-//_________________________________________________________________________
-TH1* AliTRDCalPID::GetHistogramT(Int_t k, Int_t ip) const
-{
-  //
-  // Returns one selected time bin max histogram
-  //
 
-  if (k < 0 || k >= AliPID::kSPECIES) return 0x0;
-       if(ip<0 || ip>= kNMom ) 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)*kNMom+ip);
-}
-
-
-
-//_________________________________________________________________________
-Double_t AliTRDCalPID::GetProbability(Int_t spec, Float_t mom, Float_t *dedx, Float_t length) 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
-  //
-
-       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<kNLength-1 && length>fTrackSegLength[ilength]){
-               ilength++;
-       }
-       Int_t imom = 1;
-  while(imom<kNMom-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();
-       Float_t x = dedx[0]+dedx[1], y = dedx[2];
-        Bool_t kX = (x < ax->GetBinUpEdge(nbinsx));
-       Bool_t kY = (y < ay->GetBinUpEdge(nbinsy));
-       if(kX)
-               if(kY) LQ1 = hist->GetBinContent( hist->FindBin(x, y)); 
-    //fEstimator->Estimate2D2(hist, x, y);
-               else LQ1 = hist->GetBinContent(ax->FindBin(x), nbinsy);
-       else
-               if(kY) LQ1 = hist->GetBinContent(nbinsx, ay->FindBin(y));
-               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(x, y)); 
-    //fEstimator->Estimate2D2(hist, x, y);
-               else LQ2 = hist->GetBinContent(ax->FindBin(x), nbinsy);
-       else
-               if(kY) LQ2 = hist->GetBinContent(nbinsx, ay->FindBin(y));
-               else LQ2 = hist->GetBinContent(nbinsx, nbinsy);
-       
-        // return interpolation over momentum binning
-        if(mom < fTrackMomentum[0]) return LQ1;
-        else if(mom > fTrackMomentum[kNMom-1]) return LQ2;
-        else return LQ1 + (LQ2 - LQ1)*(mom - mom1)/(mom2 - mom1);
-
-}
-
-//_________________________________________________________________________
-Double_t AliTRDCalPID::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 (spec) (0 for e) and (2 for pi)
-  // from the precalculated timbin distributions 
-  //
-  
-       if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
-
-  Int_t iTBin = timbin+1;
-  
-  // Everything which is not an electron counts as a pion for time bin max
-  //if(spec != AliPID::kElectron) spec = AliPID::kPion;
-  
-
-  
-       Int_t imom = 1;
-  while(imom<kNMom-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)*kNMom+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)*kNMom+imom-1));
-               return 0.;
-       }
-       Double_t LQ1 = hist->GetBinContent(iTBin);
-
-       if(!(hist = (TH1F*) fHistTimeBin->At(((spec==AliPID::kElectron)?0:1)*kNMom+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)*kNMom+imom));
-               return LQ1;
-       }
-       Double_t LQ2 = hist->GetBinContent(iTBin);
-
-       // return interpolation over momentum binning
-  return LQ1 + (LQ2 - LQ1)*(mom - mom1)/(mom2 - mom1);
 }
-
index 5dfdbb0..9440134 100644 (file)
@@ -7,12 +7,11 @@
 
 ///////////////////////////////////////////////////////////////////////////////
 //                                                                           //
-// 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>                                         //
+//                                                                           //
+//  Alex Bercuci <A.Bercuci@gsi.de>                                          //
+//  Alex Wilk <wilka@uni-muenster.de>                                        //
+//  Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de>                     //
 //                                                                           //
 ///////////////////////////////////////////////////////////////////////////////
 
 #include <TNamed.h>
 #endif
 
-class TH1;
-class TObjArray;
-class AliTRDCalPID : public TNamed {
-public:
-       enum {
-               kNMom = 11,
-               kNLength = 4
-       };
+class AliTRDCalPID : public TNamed
+{
+
+ public:
+
+  enum {
+    kNMom   = 11,
+    kNPlane = 6
+  };
 
   AliTRDCalPID();
   AliTRDCalPID(const Text_t *name, const Text_t *title);
-  AliTRDCalPID(const AliTRDCalPID& pd);
-  virtual        ~AliTRDCalPID();
-  AliTRDCalPID&   operator=(const AliTRDCalPID &c);
-  virtual void    Copy(TObject &c) const;
-           Bool_t   LoadLQReferences(Char_t* refFile);
-           Bool_t   LoadNNReferences(Char_t* /*refFile*/) {return kTRUE;}
-  //void         SetMeanChargeRatio(Double_t ratio)     { fMeanChargeRatio = ratio;  }
-
-  //Double_t     GetMeanChargeRatio() const             { return fMeanChargeRatio;   }
-  static   Double_t GetMomentum(Int_t ip)            { return (ip<0 || ip>=kNMom)    ? -1. : fTrackMomentum[ip];  }
-  static   Double_t GetLength(Int_t il)              { return (il<0 || il>=kNLength) ? -1. : fTrackSegLength[il]; }
-  //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;
-           Double_t GetProbability(Int_t spec, Float_t mom, Float_t *dedx, Float_t length) const;
-           Double_t GetProbabilityT(Int_t spec, Double_t mom, Int_t timbin) const;
+  virtual         ~AliTRDCalPID();
+
+  virtual Bool_t   LoadReferences(Char_t *refFile) = 0;
+  static  Double_t GetMomentum(Int_t ip)            
+                                                      { return (ip<0 || ip>=kNMom) 
+                                                        ? -1.0 
+                                                        : fTrackMomentum[ip]; 
+                                                      }
+  virtual TObject *GetModel(Int_t ip, Int_t iType, Int_t iPlane) const = 0;
+  virtual Double_t GetProbability(Int_t spec, Float_t mom, Float_t *dedx
+                                , Float_t length, Int_t plane) const = 0;
+  static  Char_t  *GetPartName(Int_t i)               { return fPartName[i]; }
+  static  Char_t  *GetPartSymb(Int_t i)               { return fPartSymb[i]; }
+
+          void     SetPartName(Int_t i, Char_t *name) { fPartName[i] = name; }
+          void     SetPartSymb(Int_t i, Char_t *symb) { fPartSymb[i] = symb; }
 
  protected:
 
-           void     Init();
-  inline  Int_t     GetHistID(Int_t part, Int_t mom/*, Int_t length=0*/) const { return part*kNMom + mom; }
-           void     CleanUp();
- public:
-  static  Char_t   *fpartName[5];     //! Names of particle species
-  static  Char_t   *fpartSymb[5];     //! Symbols of particle species
-  
+  virtual void     Init() = 0;
+  virtual Int_t    GetModelID(Int_t mom, Int_t spec, Int_t plane) const = 0;
+
+ private:
+
+  AliTRDCalPID(const AliTRDCalPID& pd);
+  AliTRDCalPID    &operator=(const AliTRDCalPID &c);
+
  protected:
-       static Float_t   fTrackMomentum[kNMom];     // Track momenta for which response functions are available
-        static Float_t   fTrackSegLength[kNLength]; // Track segment lengths for which response functions are available
-         Double_t  fMeanChargeRatio;  //  Ratio of mean charge from real Det. to prob. dist.
-         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(AliTRDCalPID, 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
 
-};
+  static  Float_t   fTrackMomentum[kNMom]; //  Track momenta for which response functions are available
+  TObjArray        *fModel;                //  Model for probability estimate
 
-#endif
+  ClassDef(AliTRDCalPID, 3)                //  Base class for TRD PID methods
 
+};
+#endif
diff --git a/TRD/Cal/AliTRDCalPIDLQ.cxx b/TRD/Cal/AliTRDCalPIDLQ.cxx
new file mode 100644 (file)
index 0000000..20c52b7
--- /dev/null
@@ -0,0 +1,231 @@
+/**************************************************************************
+ * 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$ */
+
+//////////////////////////////////////////////////////////////////////////
+//                                                                      //
+// 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 <TROOT.h>
+
+#include "AliLog.h"
+#include "AliPID.h"
+
+#include "AliTRDCalPIDLQ.h"
+#include "AliTRDcalibDB.h"
+
+ClassImp(AliTRDCalPIDLQ)
+
+Float_t AliTRDCalPIDLQ::fTrackSegLength[kNLength] = { 3.7, 3.9, 4.2, 5.0 };
+
+//_________________________________________________________________________
+AliTRDCalPIDLQ::AliTRDCalPIDLQ()
+  :AliTRDCalPID("pid", "LQ PID references for TRD")
+{
+  //
+  //  The Default constructor
+  //
+
+  Init();
+
+}
+
+//_________________________________________________________________________
+AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t *name, const Text_t *title)
+  :AliTRDCalPID(name,title)
+{
+  //
+  //  The main constructor
+  //
+  
+  Init();
+
+}
+
+//_________________________________________________________________________
+AliTRDCalPIDLQ::~AliTRDCalPIDLQ()
+{
+  //
+  // Destructor
+  //
+  
+}
+
+//_________________________________________________________________________
+Bool_t AliTRDCalPIDLQ::LoadReferences(Char_t *refFile)
+{
+  //
+  // Read the TRD dEdx histograms.
+  //
+
+       Int_t nTimeBins = 22;
+       // Get number of time bins from CDB
+       AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
+       if(!calibration){
+               AliWarning(Form("No AliTRDcalibDB available. Using %d time bins.", nTimeBins));
+       }else{
+               if(calibration->GetRun() > -1) nTimeBins = calibration->GetNumberOfTimeBins();
+               else AliWarning(Form("Run number not set. Using %d time bins.", nTimeBins));
+       }
+
+       
+  // Read histogram Root file  
+  TFile *histFile = TFile::Open(refFile, "READ");
+  if (!histFile || !histFile->IsOpen()) {
+    AliError(Form("Opening TRD histgram file %s failed", refFile));
+    return kFALSE;
+  }
+  gROOT->cd();
+
+  // Read histograms
+  for (Int_t iparticle = 0; iparticle < AliPID::kSPECIES; iparticle++){
+    for (Int_t imom = 0; imom < kNMom; imom++){
+                       TH2D* hist = (TH2D*)histFile->Get(Form("h2dEdx%s%d", fPartSymb[iparticle], imom/*, ilength*/))->Clone();
+                       hist->Scale(1./hist->Integral());
+                       fModel->AddAt(hist, GetModelID(imom, iparticle, 0));
+
+//                     if (iparticle != AliPID::kElectron && iparticle != AliPID::kPion) continue;
+// 
+//                     TH1F* ht = (TH1F*)histFile->Get(Form("h1MaxTB%s%02d", fPartSymb[iparticle], imom))->Clone();
+//                     if(ht->GetNbinsX() != nTimeBins) 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, nTimeBins));
+//                     ht->Scale(1./ht->Integral());
+//                     fHistTimeBin->AddAt(ht, ((iparticle==AliPID::kElectron)?0:1)*kNMom + 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);
+  
+  return kTRUE;
+
+
+}
+
+//_________________________________________________________________________
+TObject* AliTRDCalPIDLQ::GetModel(Int_t ip, Int_t iType, Int_t iplane) const          // iType not needed
+{
+  //
+  // Returns one selected dEdx histogram
+  //
+
+  if (iType < 0 || iType >= AliPID::kSPECIES) return 0x0;
+       if(ip<0 || ip>= kNMom ) return 0x0;
+
+       AliInfo(Form("Retrive dEdx histogram for %s of %5.2f GeV/c", fPartName[iType], fTrackMomentum[ip]));
+  
+  return fModel->At(GetModelID(ip, iType, iplane));
+}
+
+//_________________________________________________________________________
+Double_t AliTRDCalPIDLQ::GetProbability(Int_t spec, Float_t mom, Float_t *dedx, Float_t length, Int_t iplane) 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
+  //
+
+       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<kNLength-1 && length>fTrackSegLength[ilength]){
+               ilength++;
+       }
+       Int_t imom = 1;
+  while(imom<kNMom-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*)fModel->At(GetModelID(imom-1, spec, iplane)))){
+               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.", GetModelID(imom-1, spec, iplane)));
+               return 0.;
+       }
+       ax = hist->GetXaxis(); nbinsx = ax->GetNbins();
+       ay = hist->GetYaxis(); nbinsy = ay->GetNbins();
+       Float_t x = dedx[0]+dedx[1], y = dedx[2];
+        Bool_t kX = (x < ax->GetBinUpEdge(nbinsx));
+       Bool_t kY = (y < ay->GetBinUpEdge(nbinsy));
+       if(kX)
+               if(kY) LQ1 = hist->GetBinContent( hist->FindBin(x, y)); 
+    //fEstimator->Estimate2D2(hist, x, y);
+               else LQ1 = hist->GetBinContent(ax->FindBin(x), nbinsy);
+       else
+               if(kY) LQ1 = hist->GetBinContent(nbinsx, ay->FindBin(y));
+               else LQ1 = hist->GetBinContent(nbinsx, nbinsy);
+
+
+       if(!(hist = (TH2D*)fModel->At(GetModelID(imom, spec, iplane)))){
+               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.", GetModelID(imom, spec, iplane)));
+               return LQ1;
+       }
+       if(kX)
+               if(kY) LQ2 = hist->GetBinContent( hist->FindBin(x, y)); 
+    //fEstimator->Estimate2D2(hist, x, y);
+               else LQ2 = hist->GetBinContent(ax->FindBin(x), nbinsy);
+       else
+               if(kY) LQ2 = hist->GetBinContent(nbinsx, ay->FindBin(y));
+               else LQ2 = hist->GetBinContent(nbinsx, nbinsy);
+       
+        // return interpolation over momentum binning
+        if(mom < fTrackMomentum[0]) return LQ1;
+        else if(mom > fTrackMomentum[kNMom-1]) return LQ2;
+        else return LQ1 + (LQ2 - LQ1)*(mom - mom1)/(mom2 - mom1);
+
+}
+
+//_________________________________________________________________________
+void AliTRDCalPIDLQ::Init()
+{
+  //
+  // Initialization
+  //
+
+  fModel = new TObjArray(AliPID::kSPECIES  * kNMom);
+  fModel -> SetOwner();
+
+}
+
+//_________________________________________________________________________
+Int_t AliTRDCalPIDLQ::GetModelID(Int_t mom, Int_t spec, Int_t) const
+{  
+  // returns the ID of the LQ distribution (55 Histos, ID from 1 to 55)
+
+       return spec * AliTRDCalPID::kNMom + mom;
+}
diff --git a/TRD/Cal/AliTRDCalPIDLQ.h b/TRD/Cal/AliTRDCalPIDLQ.h
new file mode 100644 (file)
index 0000000..c86c783
--- /dev/null
@@ -0,0 +1,57 @@
+#ifndef ALITRDCALPIDLQ_H
+#define ALITRDCALPIDLQ_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+// PID distributions for the LQ method                                    //
+//                                                                        //
+// Author:                                                                //
+// Alex Bercuci <A.Bercuci@gsi.de>                                        //
+//                                                                        //          
+////////////////////////////////////////////////////////////////////////////
+
+#ifndef ALITRDCALPID_H
+#include "AliTRDCalPID.h"
+#endif
+
+class AliTRDCalPIDLQ : public AliTRDCalPID
+{
+
+ public:
+
+  enum {
+    kNLength = 4
+  };
+       
+  AliTRDCalPIDLQ();
+  AliTRDCalPIDLQ(const Text_t *name, const Text_t *title);
+  virtual        ~AliTRDCalPIDLQ();
+
+  Bool_t          LoadReferences(Char_t* refFile);
+  TObject*        GetModel(Int_t ip, Int_t iType, Int_t iPlane) const;
+  static Double_t GetLength(Int_t il) { return (il<0 || il>=kNLength) ? -1. : fTrackSegLength[il]; }
+         Double_t GetProbability(Int_t spec, Float_t mom, Float_t *dedx
+                               , Float_t length, Int_t plane) const;
+
+ private:
+
+  AliTRDCalPIDLQ(const AliTRDCalPIDLQ& pd);
+  AliTRDCalPIDLQ&   operator=(const AliTRDCalPIDLQ &c);
+
+  void     Init();
+  Int_t    GetModelID(Int_t mom, Int_t , Int_t) const;
+
+ protected:
+
+  static Float_t  fTrackSegLength[kNLength];  // Track segment lengths
+
+  ClassDef(AliTRDCalPIDLQ, 1)                 // LQ PID reference manager
+
+};
+#endif
+
diff --git a/TRD/Cal/AliTRDCalPIDNN.cxx b/TRD/Cal/AliTRDCalPIDNN.cxx
new file mode 100644 (file)
index 0000000..078ef81
--- /dev/null
@@ -0,0 +1,207 @@
+/**************************************************************************
+ * 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$ */
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+// Container of the distributions for the neural network                  //
+//                                                                        //
+// Author:                                                                //
+// Alex Wilk <wilka@uni-muenster.de>                                      //
+//                                                                        //
+////////////////////////////////////////////////////////////////////////////
+
+#include <TFile.h>
+#include <TROOT.h>
+#include <TObject.h>
+#include <TMultiLayerPerceptron.h>
+
+#include "AliLog.h"
+#include "AliPID.h"
+#include "AliESD.h"
+#include "AliESDtrack.h"
+#include "AliTRDtrack.h"
+
+#include "AliTRDCalPIDNN.h"
+#include "AliTRDcalibDB.h"
+
+ClassImp(AliTRDCalPIDNN)
+
+//_________________________________________________________________________
+AliTRDCalPIDNN::AliTRDCalPIDNN()
+  :AliTRDCalPID("pid", "NN PID references for TRD")
+{
+  //
+  //  The Default constructor
+  //
+
+  Init();
+
+}
+
+//_________________________________________________________________________
+AliTRDCalPIDNN::AliTRDCalPIDNN(const Text_t *name, const Text_t *title) 
+  :AliTRDCalPID(name,title)
+{
+  //
+  //  The main constructor
+  //
+  
+  Init();
+
+}
+
+// //_____________________________________________________________________________
+// AliTRDCalPIDNN::AliTRDCalPIDNN(const AliTRDCalPIDNN &c) 
+//   :TNamed(c)
+//   ,fMeanChargeRatio(c.fMeanChargeRatio)
+//   ,fModel(0x0)
+// {
+//   //
+//   // Copy constructor
+//   //
+// 
+//   if (this != &c) ((AliTRDCalPIDNN &) c).Copy(*this);
+//   
+// }
+
+//_________________________________________________________________________
+AliTRDCalPIDNN::~AliTRDCalPIDNN()
+{
+  //
+  // Destructor
+  //
+  
+}
+
+//_________________________________________________________________________
+Bool_t AliTRDCalPIDNN::LoadReferences(Char_t *refFile)
+{
+  //
+  // Read the TRD Neural Networks
+  //
+
+  // Read NN Root file  
+  TFile *nnFile = TFile::Open(refFile,"READ");
+  if (!nnFile || !nnFile->IsOpen()) {
+    AliError(Form("Opening TRD histgram file %s failed",refFile));
+    return kFALSE;
+  }
+  gROOT->cd();
+
+  // Read Networks
+  for (Int_t imom = 0; imom < kNMom; imom++) {
+    for (Int_t iplane = 0; iplane < kNPlane; iplane++) {
+      TMultiLayerPerceptron *nn = (TMultiLayerPerceptron *)
+         nnFile->Get(Form("NN_Mom%d_Plane%d",imom,iplane));
+      fModel->AddAt(nn,GetModelID(imom,0,iplane));
+    }
+  }
+
+  nnFile->Close();
+  delete nnFile;
+
+  return kTRUE;
+
+}
+
+//_________________________________________________________________________
+TObject *AliTRDCalPIDNN::GetModel(Int_t ip, Int_t, Int_t iplane) const
+{
+  //
+  // Returns one selected TMultiLayerPerceptron. iType not used.
+  //
+
+  if (ip<0 || ip>= kNMom) return 0x0;
+  
+  AliInfo(Form("Retrive MultiLayerPerceptron for %5.2f GeV/c for plane %d" 
+         ,fTrackMomentum[ip]
+         ,iplane));
+  
+  return fModel->At(GetModelID(ip, 0, iplane));
+
+}
+
+//_________________________________________________________________________
+Double_t AliTRDCalPIDNN::GetProbability(Int_t spec, Float_t mom, Float_t *dedx
+                                      , Float_t, Int_t iplane) const
+{
+  //
+  // Core function of AliTRDCalPID class for calculating the
+  // likelihood for species "spec" (see AliTRDtrack::kNspecie) of a
+  // given momentum "mom", a given dE/dx (slice "dedx") yield per
+  // layer in a given layer (iplane)
+  //
+
+  if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
+
+  // find the interval in momentum and track segment length which applies for this data
+  
+  Int_t imom = 1;
+  while(imom<AliTRDCalPID::kNMom-1 && mom>fTrackMomentum[imom]) imom++;
+  Double_t LNN1, LNN2;
+  Double_t mom1 = fTrackMomentum[imom-1], mom2 = fTrackMomentum[imom];
+
+  TMultiLayerPerceptron *nn = 0x0;
+  if(!(nn = (TMultiLayerPerceptron *) fModel->At(GetModelID(imom-1, spec, iplane/*, ilength*/)))){
+    //if(!(nn = (TMultiLayerPerceptron*)fModel->At(GetModelID(imom-1, spec, iplane/*, ilength*/)))){
+    AliInfo(Form("Looking for mom(%f) plane(%d)", mom-1, iplane));
+    AliError(Form("NN id %d not found in DB.", GetModelID(imom-1, spec, iplane)));
+    return 0.;
+  }
+
+  Double_t ddedx[AliTRDtrack::kNMLPslice];
+  for(int inode=0; inode<AliTRDtrack::kNMLPslice; inode++) ddedx[inode] = (Double_t)dedx[inode];
+  LNN1 = nn->Evaluate(spec, ddedx);
+  
+  if(!(nn = (TMultiLayerPerceptron*)fModel->At(GetModelID(imom, spec, iplane/*, ilength*/)))){
+    //if(!(nn = (TMultiLayerPerceptron*)fModel->At(GetModelID(imom, spec, iplane/*, ilength*/)))){
+    AliInfo(Form("Looking for mom(%f) plane(%d)", mom, iplane));
+    AliError(Form("NN id %d not found in DB.", GetModelID(imom, spec, iplane)));
+    return LNN1;
+  }
+  
+  LNN2 = nn->Evaluate(spec, ddedx);
+  
+  // return interpolation over momentum binning
+  if(mom < fTrackMomentum[0]) return LNN1;
+  else if(mom > fTrackMomentum[AliTRDCalPID::kNMom-1]) return LNN2;
+  else return LNN1 + (LNN2 - LNN1)*(mom - mom1)/(mom2 - mom1);
+  
+}
+
+//_________________________________________________________________________
+void AliTRDCalPIDNN::Init()
+{
+  //
+  // Initialization
+  //
+
+  fModel = new TObjArray(AliTRDCalPID::kNPlane * AliTRDCalPID::kNMom);
+  fModel->SetOwner();
+  
+}
+
+//_________________________________________________________________________
+Int_t AliTRDCalPIDNN::GetModelID(Int_t mom, Int_t, Int_t plane) const
+{
+  
+  // returns the ID of the NN distribution (66 MLPs, ID from 56 to 121)
+
+  return /*AliPID::kSPECIES * AliTRDCalPID::kNMom + */
+    plane * AliTRDCalPID::kNMom + mom;  
+  
+}
diff --git a/TRD/Cal/AliTRDCalPIDNN.h b/TRD/Cal/AliTRDCalPIDNN.h
new file mode 100644 (file)
index 0000000..523cdf8
--- /dev/null
@@ -0,0 +1,44 @@
+#ifndef ALITRDCALPIDNN_H
+#define ALITRDCALPIDNN_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+// PID distributions for the NN method                                    //
+//                                                                        //
+// Author:                                                                //
+// Alex Wilk <wilka@uni-muenster.de>                                      //
+//                                                                        //
+////////////////////////////////////////////////////////////////////////////
+
+#ifndef ALITRDCALPID_H
+#include "AliTRDCalPID.h"
+#endif
+
+class AliTRDCalPIDNN : public AliTRDCalPID
+{
+
+ public:
+
+  AliTRDCalPIDNN();
+  AliTRDCalPIDNN(const Text_t *name, const Text_t *title);
+  virtual  ~AliTRDCalPIDNN();
+  Bool_t    LoadReferences(Char_t *refFile);
+  TObject  *GetModel(Int_t ip, Int_t iType, Int_t iPlane) const;
+  Double_t  GetProbability(Int_t spec, Float_t mom, Float_t *dedx, Float_t length, Int_t plane) const;
+
+ private:
+
+  AliTRDCalPIDNN(const AliTRDCalPIDNN &pd);
+  AliTRDCalPIDNN &operator=(const AliTRDCalPIDNN &c);
+           
+  void     Init();
+  Int_t    GetModelID(Int_t mom, Int_t , Int_t) const;
+
+  ClassDef(AliTRDCalPIDNN, 1) // NN PID reference manager
+
+};
+#endif
index a5d430c..850fd8e 100644 (file)
@@ -178,7 +178,7 @@ Bool_t AliTRDCalPIDRefMaker::BuildLQReferences(Char_t *File, Char_t *dir)
        // Print statistics header
        Int_t nPart[AliPID::kSPECIES], nTotPart;
        printf("P[GeV/c] ");
-       for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) printf(" %s[%%] ", AliTRDCalPID::fpartSymb[ispec]);
+       for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) printf(" %s[%%] ", AliTRDCalPID::GetPartSymb(ispec));
        printf("\n-----------------------------------------------\n");
        
        Float_t trackMomentum[AliTRDCalPID::kNMom];
@@ -552,7 +552,7 @@ void  AliTRDCalPIDRefMaker::Prepare2D()
        for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
                // check PCA data
                if(!fPrinc[ispec]){
-                       AliError(Form("No data defined for %s.", AliTRDCalPID::fpartName[ispec]));
+                       AliError(Form("No data defined for %s.", AliTRDCalPID::GetPartName(ispec)));
                        return;
                }
                // build reference histograms
@@ -562,7 +562,7 @@ void  AliTRDCalPIDRefMaker::Prepare2D()
                xmin = ymin = 0.;
                xmax = 8000.; ymax = 6000.;
                if(!fH2dEdx[ispec]){
-                       fH2dEdx[ispec] = new  TH2D(Form("h2%s", AliTRDCalPID::fpartSymb[ispec]), "", nbinsx, xmin, xmax, nbinsy, ymin, ymax);
+                       fH2dEdx[ispec] = new  TH2D(Form("h2%s", AliTRDCalPID::GetPartSymb(ispec)), "", nbinsx, xmin, xmax, nbinsy, ymin, ymax);
                        fH2dEdx[ispec]->SetLineColor(color[ispec]);
                }
        }
@@ -876,8 +876,8 @@ void  AliTRDCalPIDRefMaker::SaveReferences(const Int_t mom, const char *fn)
        // save dE/dx references
        TH2 *h2 = 0x0;
        for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
-               h2 = (TH2D*)fH2dEdx[ispec]->Clone(Form("h2dEdx%s%d", AliTRDCalPID::fpartSymb[ispec], mom));
-               h2->SetTitle(Form("2D dEdx for particle %s @ %d", AliTRDCalPID::fpartName[ispec], mom));
+               h2 = (TH2D*)fH2dEdx[ispec]->Clone(Form("h2dEdx%s%d", AliTRDCalPID::GetPartSymb(ispec), mom));
+               h2->SetTitle(Form("2D dEdx for particle %s @ %d", AliTRDCalPID::GetPartName(ispec), mom));
                h2->GetXaxis()->SetTitle("dE/dx_{TRD}^{amplif} [au]");
                h2->GetYaxis()->SetTitle("dE/dx_{TRD}^{drift} [au]");
                h2->GetZaxis()->SetTitle("Entries");
index 056a4a2..4caf0ee 100644 (file)
Binary files a/TRD/Calib/PIDLQ/Run0_999999999_v0_s0.root and b/TRD/Calib/PIDLQ/Run0_999999999_v0_s0.root differ
diff --git a/TRD/Calib/PIDNN/Run0_999999999_v0_s0.root b/TRD/Calib/PIDNN/Run0_999999999_v0_s0.root
new file mode 100644 (file)
index 0000000..a9c3ed3
Binary files /dev/null and b/TRD/Calib/PIDNN/Run0_999999999_v0_s0.root differ
index 9dbb6b8..bb251c4 100644 (file)
@@ -24,7 +24,7 @@ void makePIDRefs(const char *dir = ".", const char *file="Refs.root")
 }
 
 //___________________________________________________________________
-void generatePIDDB(const char *file = "Refs.root")
+void generatePIDDB(const char *fileNN = "NN.root", const char *fileLQ = "LQ.root")
 {
 // Write TRD PID DB using the reference data from file "file"
 
@@ -34,24 +34,32 @@ void generatePIDDB(const char *file = "Refs.root")
        man->SetRun(0);
 
        AliCDBStorage *gStorLoc = man->GetStorage("local://$ALICE_ROOT");
-  if (!gStorLoc) return;
+       if (!gStorLoc) return;
 
   
-       AliTRDCalPID *pid = new AliTRDCalPID("pid", "TRD PID object");
-       pid->LoadLQReferences(file);
-
-       AliCDBMetaData *md= new AliCDBMetaData();
-  md->SetObjectClassName("AliTRDCalPIDLQ");
-  md->SetResponsible("Alex Bercuci");
-  md->SetBeamPeriod(1);
-  md->SetAliRootVersion("v4-06-HEAD"); //root version
-  md->SetComment("2D PID for TRD");
-
-       gStorLoc->Put(pid, AliCDBId("TRD/Calib/PIDLQ", 0, 0), md);
+       AliTRDCalPID *pidLQ = new AliTRDCalPIDLQ("pidLQ", "LQ TRD PID object");
+       pidLQ->LoadReferences(fileLQ);
+       AliCDBMetaData *md= new AliCDBMetaData();
+       md->SetObjectClassName("AliTRDCalPIDLQ");
+       md->SetResponsible("Alex Bercuci");
+       md->SetBeamPeriod(1);
+       md->SetAliRootVersion("v4-06-HEAD"); //root version
+       md->SetComment("2D PID for TRD");
+       gStorLoc->Put(pidLQ, AliCDBId("TRD/Calib/PIDLQ", 0, 999999999, 0, 1), md);
+       
+       AliTRDCalPID *pidNN = new AliTRDCalPIDNN("pidNN", "NN TRD PID object");
+       pidNN->LoadReferences(fileNN);
+       md->SetObjectClassName("AliTRDCalPIDNN");
+       md->SetResponsible("Alex Wilk");
+       md->SetBeamPeriod(1);
+       md->SetAliRootVersion("v4-06-HEAD"); //root version
+       md->SetComment("NN PID for TRD");
+       
+       gStorLoc->Put(pidNN, AliCDBId("TRD/Calib/PIDNN", 0, 999999999, 0, 1), md);
 }
 
 //___________________________________________________________________
-AliTRDCalPID* getPIDObject()
+AliTRDCalPID* getPIDObject(const char *method="NN")
 {
 // Returns PIDLQ object.
 // In order to browse histos do:
@@ -64,7 +72,7 @@ AliTRDCalPID* getPIDObject()
        CDBManager->SetDefaultStorage("local://$ALICE_ROOT");
        CDBManager->SetRun(0);
 
-       AliCDBEntry *wrap = CDBManager->Get("TRD/Calib/PIDLQ", 0);
+       AliCDBEntry *wrap = CDBManager->Get(Form("TRD/Calib/PID%s", method), 0);
        AliTRDCalPID *pid = dynamic_cast<const AliTRDCalPID *>wrap->GetObject();
        AliCDBMetaData *meta = wrap->GetMetaData();
        if(!pid){
index 5eb48ce..42c4023 100644 (file)
@@ -40,6 +40,8 @@
 #pragma link C++ class  AliTRDCalDet+;
 #pragma link C++ class  AliTRDCalFEE+;
 #pragma link C++ class  AliTRDCalPID+;
+#pragma link C++ class  AliTRDCalPIDLQ+;
+#pragma link C++ class  AliTRDCalPIDNN+;
 #pragma link C++ class  AliTRDCalPIDRefMaker+;
 #pragma link C++ class  AliTRDCalMonitoring+;
 
index 5b6ab29..b17db3c 100644 (file)
@@ -22,6 +22,8 @@ SRCS= AliTRDarrayI.cxx \
       Cal/AliTRDCalDet.cxx \
       Cal/AliTRDCalFEE.cxx \
       Cal/AliTRDCalPID.cxx \
+      Cal/AliTRDCalPIDLQ.cxx \
+      Cal/AliTRDCalPIDNN.cxx \
       Cal/AliTRDCalPIDRefMaker.cxx \
       Cal/AliTRDCalMonitoring.cxx \
       Cal/AliTRDCalChamberStatus.cxx \