First implementation of neural network PID
[u/mrichter/AliRoot.git] / 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 }
89
90 //_________________________________________________________________________
91 Bool_t AliTRDCalPIDNN::LoadReferences(Char_t *refFile)
92 {
93   //
94   // Read the TRD Neural Networks
95   //
96
97   // Read NN Root file  
98   TFile *nnFile = TFile::Open(refFile,"READ");
99   if (!nnFile || !nnFile->IsOpen()) {
100     AliError(Form("Opening TRD histgram file %s failed",refFile));
101     return kFALSE;
102   }
103   gROOT->cd();
104
105   // Read Networks
106   for (Int_t imom = 0; imom < kNMom; imom++) {
107     for (Int_t iplane = 0; iplane < kNPlane; iplane++) {
108       TMultiLayerPerceptron *nn = (TMultiLayerPerceptron *)
109          nnFile->Get(Form("NN_Mom%d_Plane%d",imom,iplane));
110       fModel->AddAt(nn,GetModelID(imom,0,iplane));
111     }
112   }
113
114   nnFile->Close();
115   delete nnFile;
116
117   return kTRUE;
118
119 }
120
121 //_________________________________________________________________________
122 TObject *AliTRDCalPIDNN::GetModel(Int_t ip, Int_t, Int_t iplane) const
123 {
124   //
125   // Returns one selected TMultiLayerPerceptron. iType not used.
126   //
127
128   if (ip<0 || ip>= kNMom) return 0x0;
129   
130   AliInfo(Form("Retrive MultiLayerPerceptron for %5.2f GeV/c for plane %d" 
131          ,fTrackMomentum[ip]
132          ,iplane));
133   
134   return fModel->At(GetModelID(ip, 0, iplane));
135
136 }
137
138 //_________________________________________________________________________
139 Double_t AliTRDCalPIDNN::GetProbability(Int_t spec, Float_t mom, Float_t *dedx
140                                       , Float_t, Int_t iplane) const
141 {
142   //
143   // Core function of AliTRDCalPID class for calculating the
144   // likelihood for species "spec" (see AliTRDtrack::kNspecie) of a
145   // given momentum "mom", a given dE/dx (slice "dedx") yield per
146   // layer in a given layer (iplane)
147   //
148
149   if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
150
151   // find the interval in momentum and track segment length which applies for this data
152   
153   Int_t imom = 1;
154   while(imom<AliTRDCalPID::kNMom-1 && mom>fTrackMomentum[imom]) imom++;
155   Double_t LNN1, LNN2;
156   Double_t mom1 = fTrackMomentum[imom-1], mom2 = fTrackMomentum[imom];
157
158   TMultiLayerPerceptron *nn = 0x0;
159   if(!(nn = (TMultiLayerPerceptron *) fModel->At(GetModelID(imom-1, spec, iplane/*, ilength*/)))){
160     //if(!(nn = (TMultiLayerPerceptron*)fModel->At(GetModelID(imom-1, spec, iplane/*, ilength*/)))){
161     AliInfo(Form("Looking for mom(%f) plane(%d)", mom-1, iplane));
162     AliError(Form("NN id %d not found in DB.", GetModelID(imom-1, spec, iplane)));
163     return 0.;
164   }
165
166   Double_t ddedx[AliTRDtrack::kNMLPslice];
167   for(int inode=0; inode<AliTRDtrack::kNMLPslice; inode++) ddedx[inode] = (Double_t)dedx[inode];
168   LNN1 = nn->Evaluate(spec, ddedx);
169   
170   if(!(nn = (TMultiLayerPerceptron*)fModel->At(GetModelID(imom, spec, iplane/*, ilength*/)))){
171     //if(!(nn = (TMultiLayerPerceptron*)fModel->At(GetModelID(imom, spec, iplane/*, ilength*/)))){
172     AliInfo(Form("Looking for mom(%f) plane(%d)", mom, iplane));
173     AliError(Form("NN id %d not found in DB.", GetModelID(imom, spec, iplane)));
174     return LNN1;
175   }
176   
177   LNN2 = nn->Evaluate(spec, ddedx);
178   
179   // return interpolation over momentum binning
180   if(mom < fTrackMomentum[0]) return LNN1;
181   else if(mom > fTrackMomentum[AliTRDCalPID::kNMom-1]) return LNN2;
182   else return LNN1 + (LNN2 - LNN1)*(mom - mom1)/(mom2 - mom1);
183   
184 }
185
186 //_________________________________________________________________________
187 void AliTRDCalPIDNN::Init()
188 {
189   //
190   // Initialization
191   //
192
193   fModel = new TObjArray(AliTRDCalPID::kNPlane * AliTRDCalPID::kNMom);
194   fModel->SetOwner();
195   
196 }
197
198 //_________________________________________________________________________
199 Int_t AliTRDCalPIDNN::GetModelID(Int_t mom, Int_t, Int_t plane) const
200 {
201   
202   // returns the ID of the NN distribution (66 MLPs, ID from 56 to 121)
203
204   return /*AliPID::kSPECIES * AliTRDCalPID::kNMom + */
205     plane * AliTRDCalPID::kNMom + mom;  
206   
207 }