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 "AliTRDCalPIDNN.h"
  39 #include "AliTRDcalibDB.h"
  40
  41 ClassImp(AliTRDCalPIDNN)
  42
  43 //_________________________________________________________________________
  44 AliTRDCalPIDNN::AliTRDCalPIDNN()
  45   :AliTRDCalPID("pid", "NN PID references for TRD")
  46 {
  47   //
  48   //  The Default constructor
  49   //
  50
  51   Init();
  52
  53 }
  54
  55 //_________________________________________________________________________
  56 AliTRDCalPIDNN::AliTRDCalPIDNN(const Text_t *name, const Text_t *title)
  57   :AliTRDCalPID(name,title)
  58 {
  59   //
  60   //  The main constructor
  61   //
  62
  63   Init();
  64
  65 }
  66
  67 // //_____________________________________________________________________________
  68 // AliTRDCalPIDNN::AliTRDCalPIDNN(const AliTRDCalPIDNN &c)
  69 //   :TNamed(c)
  70 //   ,fMeanChargeRatio(c.fMeanChargeRatio)
  71 //   ,fModel(0x0)
  72 // {
  73 //   //
  74 //   // Copy constructor
  75 //   //
  76 //
  77 //   if (this != &c) ((AliTRDCalPIDNN &) c).Copy(*this);
  78 //
  79 // }
  80
  81 //_________________________________________________________________________
  82 AliTRDCalPIDNN::~AliTRDCalPIDNN()
  83 {
  84   //
  85   // Destructor
  86   //
  87
  88   //if (fModel) delete fModel;
  89
  90 }
  91
  92 //_________________________________________________________________________
  93 Bool_t AliTRDCalPIDNN::LoadReferences(Char_t *refFile)
  94 {
  95   //
  96   // Read the TRD Neural Networks
  97   //
  98
  99   // Read NN Root file
 100   TFile *nnFile = TFile::Open(refFile,"READ");
 101   if (!nnFile || !nnFile->IsOpen()) {
 102     AliError(Form("Opening TRD histgram file %s failed",refFile));
 103     return kFALSE;
 104   }
 105   gROOT->cd();
 106
 107   // Read Networks
 108   for (Int_t imom = 0; imom < kNMom; imom++) {
 109     for (Int_t iplane = 0; iplane < kNPlane; iplane++) {
 110       TMultiLayerPerceptron *nn = (TMultiLayerPerceptron *)
 111          nnFile->Get(Form("NN_Mom%d_Plane%d",imom,iplane));
 112       fModel->AddAt(nn,GetModelID(imom,0,iplane));
 113     }
 114   }
 115
 116   nnFile->Close();
 117   delete nnFile;
 118
 119   return kTRUE;
 120
 121 }
 122
 123 //_________________________________________________________________________
 124 TObject *AliTRDCalPIDNN::GetModel(Int_t ip, Int_t, Int_t iplane) const
 125 {
 126   //
 127   // Returns one selected TMultiLayerPerceptron. iType not used.
 128   //
 129
 130   if (ip<0 || ip>= kNMom) return 0x0;
 131
 132   AliInfo(Form("Retrive MultiLayerPerceptron for %5.2f GeV/c for plane %d"
 133          ,fTrackMomentum[ip]
 134          ,iplane));
 135
 136   return fModel->At(GetModelID(ip, 0, iplane));
 137
 138 }
 139
 140 //_________________________________________________________________________
 141 Double_t AliTRDCalPIDNN::GetProbability(Int_t spec, Float_t mom, Float_t *dedx
 142                                       , Float_t, Int_t iplane) const
 143 {
 144   //
 145   // Core function of AliTRDCalPID class for calculating the
 146   // likelihood for species "spec" (see AliTRDtrack::kNspecie) of a
 147   // given momentum "mom", a given dE/dx (slice "dedx") yield per
 148   // layer in a given layer (iplane)
 149   //
 150
 151   if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
 152
 153   // find the interval in momentum and track segment length which applies for this data
 154
 155   Int_t imom = 1;
 156   while (imom<AliTRDCalPID::kNMom-1 && mom>fTrackMomentum[imom]) imom++;
 157   Double_t lNN1, lNN2;
 158   Double_t mom1 = fTrackMomentum[imom-1], mom2 = fTrackMomentum[imom];
 159
 160   TMultiLayerPerceptron *nn = 0x0;
 161   if(!(nn = (TMultiLayerPerceptron *) fModel->At(GetModelID(imom-1, spec, iplane/*, ilength*/)))){
 162     AliInfo(Form("Looking for mom(%f) plane(%d)", mom-1, iplane));
 163     AliError(Form("NN id %d not found in DB.", GetModelID(imom-1, spec, iplane)));
 164     return 0.;
 165   }
 166
 167   Double_t ddedx[AliTRDtrack::kNMLPslice];
 168
 169   for (int inode=0; inode<AliTRDtrack::kNMLPslice; inode++) {
 170     ddedx[inode] = (((Double_t) dedx[inode]/kMLPscale)*6)          // Bug fix! Needs new reference data or different calculation of dedx!!!!
 171                  / (AliTRDcalibDB::Instance()->GetNumberOfTimeBins()/AliTRDtrack::kNMLPslice);
 172   }
 173
 174   lNN1 = nn->Evaluate(spec, ddedx);
 175
 176   if(!(nn = (TMultiLayerPerceptron*)fModel->At(GetModelID(imom, spec, iplane/*, ilength*/)))){
 177     AliInfo(Form("Looking for mom(%f) plane(%d)", mom, iplane));
 178     AliError(Form("NN id %d not found in DB.", GetModelID(imom, spec, iplane)));
 179     return lNN1;
 180   }
 181   lNN2 = nn->Evaluate(spec, ddedx);
 182
 183   // return interpolation over momentum binning
 184   if      (mom < fTrackMomentum[0]) {
 185     return lNN1;
 186   }
 187   else if (mom > fTrackMomentum[AliTRDCalPID::kNMom-1]) {
 188     return lNN2;
 189   }
 190   else {
 191     return lNN1 + (lNN2 - lNN1)*(mom - mom1)/(mom2 - mom1);
 192   }
 193
 194 }
 195
 196 //_________________________________________________________________________
 197 void AliTRDCalPIDNN::Init()
 198 {
 199   //
 200   // Initialization
 201   //
 202
 203   fModel = new TObjArray(AliTRDCalPID::kNPlane * AliTRDCalPID::kNMom);
 204   fModel->SetOwner();
 205
 206 }
 207
 208 //_________________________________________________________________________
 209 Int_t AliTRDCalPIDNN::GetModelID(Int_t mom, Int_t, Int_t plane) const
 210 {
 211
 212   // returns the ID of the NN distribution (66 MLPs, ID from 56 to 121)
 213
 214   return /*AliPID::kSPECIES * AliTRDCalPID::kNMom + */
 215     plane * AliTRDCalPID::kNMom + mom;
 216
 217 }