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