TRD module
[u/mrichter/AliRoot.git] / TRD / TRDrec / AliTRDCalPIDNN.cxx
CommitLineData
44dbae42 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>
44dbae42 29#include <TMultiLayerPerceptron.h>
30
44dbae42 31#include "AliPID.h"
4d6aee34 32#include "AliLog.h"
44dbae42 33
5d6dc395 34#include "AliTRDgeometry.h"
44dbae42 35#include "AliTRDCalPIDNN.h"
44dbae42 36
37ClassImp(AliTRDCalPIDNN)
38
39//_________________________________________________________________________
40AliTRDCalPIDNN::AliTRDCalPIDNN()
41 :AliTRDCalPID("pid", "NN PID references for TRD")
42{
43 //
44 // The Default constructor
45 //
46
47 Init();
48
49}
50
51//_________________________________________________________________________
52AliTRDCalPIDNN::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//_________________________________________________________________________
78AliTRDCalPIDNN::~AliTRDCalPIDNN()
79{
80 //
81 // Destructor
82 //
e5f336c9 83
84 //if (fModel) delete fModel;
44dbae42 85
86}
87
88//_________________________________________________________________________
89Bool_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++) {
5d6dc395 105 for (Int_t iplane = 0; iplane < AliTRDgeometry::kNlayer; iplane++) {
44dbae42 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//_________________________________________________________________________
120TObject *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"
2e32a5ae 129 ,fgTrackMomentum[ip]
44dbae42 130 ,iplane));
131
132 return fModel->At(GetModelID(ip, 0, iplane));
133
134}
135
136//_________________________________________________________________________
2e32a5ae 137Double_t AliTRDCalPIDNN::GetProbability(Int_t spec, Float_t mom
138 , const Float_t * const dedx
44dbae42 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;
2e32a5ae 153 while (imom<AliTRDCalPID::kNMom-1 && mom>fgTrackMomentum[imom]) imom++;
e5f336c9 154 Double_t lNN1, lNN2;
2e32a5ae 155 Double_t mom1 = fgTrackMomentum[imom-1], mom2 = fgTrackMomentum[imom];
44dbae42 156
157 TMultiLayerPerceptron *nn = 0x0;
158 if(!(nn = (TMultiLayerPerceptron *) fModel->At(GetModelID(imom-1, spec, iplane/*, ilength*/)))){
44dbae42 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
4d6aee34 164 Double_t ddedx[AliTRDCalPID::kNSlicesNN];
28efdace 165
4d6aee34 166 for (int inode=0; inode<AliTRDCalPID::kNSlicesNN; inode++) {
e3d45279 167 ddedx[inode] = (Double_t) dedx[inode]/kMLPscale;
e5f336c9 168 }
28efdace 169
e5f336c9 170 lNN1 = nn->Evaluate(spec, ddedx);
44dbae42 171
172 if(!(nn = (TMultiLayerPerceptron*)fModel->At(GetModelID(imom, spec, iplane/*, ilength*/)))){
44dbae42 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)));
e5f336c9 175 return lNN1;
44dbae42 176 }
e5f336c9 177 lNN2 = nn->Evaluate(spec, ddedx);
44dbae42 178
179 // return interpolation over momentum binning
2e32a5ae 180 if (mom < fgTrackMomentum[0]) {
e5f336c9 181 return lNN1;
182 }
2e32a5ae 183 else if (mom > fgTrackMomentum[AliTRDCalPID::kNMom-1]) {
e5f336c9 184 return lNN2;
185 }
186 else {
187 return lNN1 + (lNN2 - lNN1)*(mom - mom1)/(mom2 - mom1);
188 }
44dbae42 189
190}
191
192//_________________________________________________________________________
193void AliTRDCalPIDNN::Init()
194{
195 //
196 // Initialization
197 //
198
5d6dc395 199 fModel = new TObjArray(AliTRDgeometry::kNlayer * AliTRDCalPID::kNMom);
44dbae42 200 fModel->SetOwner();
201
202}
203
204//_________________________________________________________________________
9ded305e 205Int_t AliTRDCalPIDNN::GetModelID(Int_t mom, Int_t /*ii*/, Int_t plane)
44dbae42 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}