new task introduced by Alex Wilk (calculation of reference data for NN
[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   //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 }