TOF + macros to each detector folder
[u/mrichter/AliRoot.git] / TRD / TRDrec / 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 }