TRD/Cal/AliTRDCalPIDNN.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 ////////////////////////////////////////////////////////////////////////////
  19 //                                                                        //
  20 // Container of the distributions for the neural network                  //
  21 //                                                                        //
  22 // Author:                                                                //
  23 // Alex Wilk <wilka@uni-muenster.de>                                      //
  24 //                                                                        //
  25 ////////////////////////////////////////////////////////////////////////////
  26
  27 #include <TFile.h>
  28 #include <TROOT.h>
  29 #include <TObject.h>
  30 #include <TMultiLayerPerceptron.h>
  31
  32 #include "AliLog.h"
  33 #include "AliPID.h"
  34 #include "AliESD.h"
  35 #include "AliESDtrack.h"
  36 #include "AliTRDtrack.h"
  37
  38 #include "AliTRDgeometry.h"
  39 #include "AliTRDCalPIDNN.h"
  40 #include "AliTRDcalibDB.h"
  41
  42 ClassImp(AliTRDCalPIDNN)
  43
  44 //_________________________________________________________________________
  45 AliTRDCalPIDNN::AliTRDCalPIDNN()
  46   :AliTRDCalPID("pid", "NN PID references for TRD")
  47 {
  48   //
  49   //  The Default constructor
  50   //
  51
  52   Init();
  53
  54 }
  55
  56 //_________________________________________________________________________
  57 AliTRDCalPIDNN::AliTRDCalPIDNN(const Text_t *name, const Text_t *title)
  58   :AliTRDCalPID(name,title)
  59 {
  60   //
  61   //  The main constructor
  62   //
  63
  64   Init();
  65
  66 }
  67
  68 // //_____________________________________________________________________________
  69 // AliTRDCalPIDNN::AliTRDCalPIDNN(const AliTRDCalPIDNN &c)
  70 //   :TNamed(c)
  71 //   ,fMeanChargeRatio(c.fMeanChargeRatio)
  72 //   ,fModel(0x0)
  73 // {
  74 //   //
  75 //   // Copy constructor
  76 //   //
  77 //
  78 //   if (this != &c) ((AliTRDCalPIDNN &) c).Copy(*this);
  79 //
  80 // }
  81
  82 //_________________________________________________________________________
  83 AliTRDCalPIDNN::~AliTRDCalPIDNN()
  84 {
  85   //
  86   // Destructor
  87   //
  88
  89   //if (fModel) delete fModel;
  90
  91 }
  92
  93 //_________________________________________________________________________
  94 Bool_t AliTRDCalPIDNN::LoadReferences(Char_t *refFile)
  95 {
  96   //
  97   // Read the TRD Neural Networks
  98   //
  99
 100   // Read NN Root file
 101   TFile *nnFile = TFile::Open(refFile,"READ");
 102   if (!nnFile || !nnFile->IsOpen()) {
 103     AliError(Form("Opening TRD histgram file %s failed",refFile));
 104     return kFALSE;
 105   }
 106   gROOT->cd();
 107
 108   // Read Networks
 109   for (Int_t imom = 0; imom < kNMom; imom++) {
 110     for (Int_t iplane = 0; iplane < AliTRDgeometry::kNlayer; iplane++) {
 111       TMultiLayerPerceptron *nn = (TMultiLayerPerceptron *)
 112          nnFile->Get(Form("NN_Mom%d_Plane%d",imom,iplane));
 113       fModel->AddAt(nn,GetModelID(imom,0,iplane));
 114     }
 115   }
 116
 117   nnFile->Close();
 118   delete nnFile;
 119
 120   return kTRUE;
 121
 122 }
 123
 124 //_________________________________________________________________________
 125 TObject *AliTRDCalPIDNN::GetModel(Int_t ip, Int_t, Int_t iplane) const
 126 {
 127   //
 128   // Returns one selected TMultiLayerPerceptron. iType not used.
 129   //
 130
 131   if (ip<0 || ip>= kNMom) return 0x0;
 132
 133   AliInfo(Form("Retrive MultiLayerPerceptron for %5.2f GeV/c for plane %d"
 134          ,fTrackMomentum[ip]
 135          ,iplane));
 136
 137   return fModel->At(GetModelID(ip, 0, iplane));
 138
 139 }
 140
 141 //_________________________________________________________________________
 142 Double_t AliTRDCalPIDNN::GetProbability(Int_t spec, Float_t mom, Float_t *dedx
 143                                       , Float_t, Int_t iplane) const
 144 {
 145   //
 146   // Core function of AliTRDCalPID class for calculating the
 147   // likelihood for species "spec" (see AliTRDtrack::kNspecie) of a
 148   // given momentum "mom", a given dE/dx (slice "dedx") yield per
 149   // layer in a given layer (iplane)
 150   //
 151
 152   if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
 153
 154   // find the interval in momentum and track segment length which applies for this data
 155
 156   Int_t imom = 1;
 157   while (imom<AliTRDCalPID::kNMom-1 && mom>fTrackMomentum[imom]) imom++;
 158   Double_t lNN1, lNN2;
 159   Double_t mom1 = fTrackMomentum[imom-1], mom2 = fTrackMomentum[imom];
 160
 161   TMultiLayerPerceptron *nn = 0x0;
 162   if(!(nn = (TMultiLayerPerceptron *) fModel->At(GetModelID(imom-1, spec, iplane/*, ilength*/)))){
 163     AliInfo(Form("Looking for mom(%f) plane(%d)", mom-1, iplane));
 164     AliError(Form("NN id %d not found in DB.", GetModelID(imom-1, spec, iplane)));
 165     return 0.;
 166   }
 167
 168   Double_t ddedx[AliTRDtrack::kNMLPslice];
 169
 170   for (int inode=0; inode<AliTRDtrack::kNMLPslice; inode++) {
 171     ddedx[inode] = (((Double_t) dedx[inode]/kMLPscale)*3)          // Bug fix! Needs new reference data or different calculation of dedx!!!!
 172                  / (AliTRDcalibDB::Instance()->GetNumberOfTimeBins()/AliTRDtrack::kNMLPslice);
 173   }
 174
 175   lNN1 = nn->Evaluate(spec, ddedx);
 176
 177   if(!(nn = (TMultiLayerPerceptron*)fModel->At(GetModelID(imom, spec, iplane/*, ilength*/)))){
 178     AliInfo(Form("Looking for mom(%f) plane(%d)", mom, iplane));
 179     AliError(Form("NN id %d not found in DB.", GetModelID(imom, spec, iplane)));
 180     return lNN1;
 181   }
 182   lNN2 = nn->Evaluate(spec, ddedx);
 183
 184   // return interpolation over momentum binning
 185   if      (mom < fTrackMomentum[0]) {
 186     return lNN1;
 187   }
 188   else if (mom > fTrackMomentum[AliTRDCalPID::kNMom-1]) {
 189     return lNN2;
 190   }
 191   else {
 192     return lNN1 + (lNN2 - lNN1)*(mom - mom1)/(mom2 - mom1);
 193   }
 194
 195 }
 196
 197 //_________________________________________________________________________
 198 void AliTRDCalPIDNN::Init()
 199 {
 200   //
 201   // Initialization
 202   //
 203
 204   fModel = new TObjArray(AliTRDgeometry::kNlayer * AliTRDCalPID::kNMom);
 205   fModel->SetOwner();
 206
 207 }
 208
 209 //_________________________________________________________________________
 210 Int_t AliTRDCalPIDNN::GetModelID(Int_t mom, Int_t, Int_t plane) const
 211 {
 212
 213   // returns the ID of the NN distribution (66 MLPs, ID from 56 to 121)
 214
 215   return /*AliPID::kSPECIES * AliTRDCalPID::kNMom + */
 216     plane * AliTRDCalPID::kNMom + mom;
 217
 218 }