]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/TRD/AliTRDpidRefMakerNN.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
implement LocalInit function
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDpidRefMakerNN.cxx
CommitLineData
1ee39b3a 1/*************************************************************************
9ec5bce7 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-commercialf 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 **************************************************************************/
1ee39b3a 15
16/* $Id: AliTRDpidRefMakerNN.cxx 27496 2008-07-22 08:35:45Z cblume $ */
17
18////////////////////////////////////////////////////////////////////////////
19// //
20// Builds the reference tree for the training of neural networks //
21// //
22////////////////////////////////////////////////////////////////////////////
23
24#include "TSystem.h"
25#include "TDatime.h"
26#include "TPDGCode.h"
27#include "TH1F.h"
28#include "TH2F.h"
29#include "TFile.h"
30#include "TGraphErrors.h"
31#include "TTree.h"
32#include "TEventList.h"
33#include "TMultiLayerPerceptron.h"
34
35#include "AliPID.h"
36#include "AliESDtrack.h"
37#include "AliTrackReference.h"
38
39#include "AliTRDtrackV1.h"
40#include "AliTRDReconstructor.h"
41#include "AliTRDpidUtil.h"
42#include "AliTRDpidRefMakerNN.h"
43#include "AliTRDpidUtil.h"
44
da27d983 45#include "Cal/AliTRDCalPID.h"
46#include "Cal/AliTRDCalPIDNN.h"
1ee39b3a 47#include "info/AliTRDtrackInfo.h"
48#include "info/AliTRDv0Info.h"
49
50ClassImp(AliTRDpidRefMakerNN)
51
705f8b0a 52//________________________________________________________________________
9ec5bce7 53 AliTRDpidRefMakerNN::AliTRDpidRefMakerNN()
705f8b0a 54 :AliTRDpidRefMaker()
55 ,fNet(NULL)
56 ,fTrainMomBin(kAll)
57 ,fEpochs(1000)
58 ,fMinTrain(100)
59 ,fDate(0)
60 ,fDoTraining(0)
61 ,fContinueTraining(0)
62 ,fTrainPath(NULL)
63 ,fScale(0)
64 ,fLy(0)
65 ,fNtrkl(0)
66 ,fRef(NULL)
67{
68 //
69 // Default constructor
70 //
71 SetNameTitle("refMakerNN", "PID(NN) Reference Maker");
72}
73
74//________________________________________________________________________
75 AliTRDpidRefMakerNN::AliTRDpidRefMakerNN(const char *name)
76 :AliTRDpidRefMaker(name, "PID(NN) Reference Maker")
77 ,fNet(NULL)
78 ,fTrainMomBin(kAll)
79 ,fEpochs(1000)
80 ,fMinTrain(100)
81 ,fDate(0)
82 ,fDoTraining(0)
83 ,fContinueTraining(0)
84 ,fTrainPath(NULL)
85 ,fScale(0)
86 ,fLy(0)
87 ,fNtrkl(0)
88 ,fRef(NULL)
1ee39b3a 89{
90 //
91 // Default constructor
92 //
93
9ec5bce7 94 memset(fTrain, 0, AliTRDCalPID::kNMom*sizeof(TEventList*));
95 memset(fTest, 0, AliTRDCalPID::kNMom*sizeof(TEventList*));
96 memset(fTrainData, 0, AliTRDCalPID::kNMom*sizeof(TTree*));
1ee39b3a 97
98 SetAbundance(.67);
99 SetScaledEdx(Float_t(AliTRDCalPIDNN::kMLPscale));
100 TDatime datime;
101 fDate = datime.GetDate();
1ee39b3a 102}
103
104
105//________________________________________________________________________
106AliTRDpidRefMakerNN::~AliTRDpidRefMakerNN()
107{
108}
109
110
111//________________________________________________________________________
705f8b0a 112void AliTRDpidRefMakerNN::MakeTrainTestTrees()
1ee39b3a 113{
9ec5bce7 114 // Create output file and tree
1ee39b3a 115 // Called once
116
705f8b0a 117 fRef = new TFile("TRD.CalibPIDrefMakerNN.root", "RECREATE");
9ec5bce7 118 for(Int_t ip = 0; ip < AliTRDCalPID::kNMom; ip++){
119 fTrainData[ip] = new TTree(Form("fTrainData_%d", ip), Form("NN Reference Data for MomBin %d", ip));
120 fTrainData[ip] -> Branch("fdEdx", fdEdx, Form("fdEdx[%d]/F", AliTRDpidUtil::kNNslices));
121 fTrainData[ip] -> Branch("fPID", fPID, Form("fPID[%d]/F", AliPID::kSPECIES));
122 fTrainData[ip] -> Branch("fLy", &fLy, "fLy/I");
123 fTrainData[ip] -> Branch("fNtrkl", &fNtrkl, "fNtrkl/I");
124
125 fTrain[ip] = new TEventList(Form("fTrainMom%d", ip), Form("Training list for momentum intervall %d", ip));
126 fTest[ip] = new TEventList(Form("fTestMom%d", ip), Form("Test list for momentum intervall %d", ip));
127 }
1ee39b3a 128}
129
130
131
132//________________________________________________________________________
133Bool_t AliTRDpidRefMakerNN::PostProcess()
134{
135 // Draw result to the screen
136 // Called once at the end of the query
137
9ec5bce7 138 TFile *fCalib = TFile::Open(Form("TRD.CalibPIDrefMaker.root"));
139 if (!fCalib) {
140 AliError("Calibration file not available");
141 return kFALSE;
142 }
705f8b0a 143 fData = (TTree*)fCalib->Get("refMakerNN");
9ec5bce7 144 if (!fData) {
145 AliError("Calibration data not available");
146 return kFALSE;
147 }
705f8b0a 148 TObjArray *o = NULL;
149 if(!(o = (TObjArray*)fCalib->Get("MonitorNN"))) {
9ec5bce7 150 AliWarning("Missing monitoring container.");
151 return kFALSE;
152 }
153 fContainer = (TObjArray*)o->Clone("monitor");
154
155
156 if (!fRef) {
705f8b0a 157 AliDebug(2, "Loading file TRD.CalibPIDrefMakerNN.root");
158 LoadFile("TRD.CalibPIDrefMakerNN.root");
9ec5bce7 159 }
160 else AliDebug(2, "file available");
161
162 if (!fRef) {
705f8b0a 163 MakeTrainTestTrees();
9ec5bce7 164
165 // Convert the CaliPIDrefMaker tree to 11 (different momentum bin) trees for NN training
166
167 LinkPIDdata();
168 for(Int_t ip=0; ip < AliTRDCalPID::kNMom; ip++){
169 for(Int_t is=0; is < AliPID::kSPECIES; is++) {
170 memset(fPID, 0, AliPID::kSPECIES*sizeof(Float_t));
171 fPID[is] = 1;
172 Int_t n(0); // index of data
173 for(Int_t itrk = 0; itrk<fData->GetEntries() && n<kMaxStat; itrk++){
174 if(!(fData->GetEntry(itrk))) continue;
175 if(fPIDbin!=is) continue;
176 for(Int_t ily=fPIDdataArray->fNtracklets; ily--;){
177 fLy = ily;
178 fNtrkl = fPIDdataArray->fNtracklets;
179 if((fPIDdataArray->fData[ily].fPLbin & 0xf)!= ip) continue;
180 memset(fdEdx, 0, AliTRDpidUtil::kNNslices*sizeof(Float_t));
181 for(Int_t islice=AliTRDCalPID::kNSlicesNN; islice--;){
182 fdEdx[islice]+=fPIDdataArray->fData[ily].fdEdx[islice];
183 fdEdx[islice]/=fScale;
184 }
185 fTrainData[ip] -> Fill();
186 n++;
187 }
188 }
189 AliDebug(2, Form("%d %d %d", ip, is, n));
190 }
191 }
192
193
194 fRef -> cd();
195 for(Int_t ip = 0; ip < AliTRDCalPID::kNMom; ip++){
196 fTrainData[ip] -> Write();
197 }
198
199 }
200 else AliDebug(2, "file available");
201
202
203 // build the training and the test list for the neural networks
204 for(Int_t ip = 0; ip < AliTRDCalPID::kNMom; ip++){
205 MakeTrainingLists(ip);
206 }
1ee39b3a 207 if(!fDoTraining) return kTRUE;
208
9ec5bce7 209
210
211 // train the neural networks
1ee39b3a 212 gSystem->Exec(Form("mkdir ./Networks_%d/",fDate));
213 AliDebug(2, Form("TrainMomBin [%d] [%d]", fTrainMomBin, kAll));
214
215 // train single network for a single momentum (recommended)
216 if(!(fTrainMomBin == kAll)){
9ec5bce7 217 if(fTrain[fTrainMomBin] -> GetN() < fMinTrain){
218 AliError("Warning in AliTRDpidRefMakerNN::PostProcess : Not enough events for training available! Please check Data sample!");
1ee39b3a 219 return kFALSE;
220 }
221 MakeRefs(fTrainMomBin);
1ee39b3a 222 MonitorTraining(fTrainMomBin);
223 }
224 // train all momenta
225 else{
226 for(Int_t iMomBin = 0; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
9ec5bce7 227 if(fTrain[iMomBin] -> GetN() < fMinTrain){
228 AliError(Form("Warning in AliTRDpidRefMakerNN::PostProcess : Not enough events for training available for momentum bin [%d]! Please check Data sample!", iMomBin));
229 continue;
1ee39b3a 230 }
231 MakeRefs(fTrainMomBin);
1ee39b3a 232 MonitorTraining(iMomBin);
233 }
234 }
235
236 return kTRUE; // testing protection
237}
238
239
240//________________________________________________________________________
9ec5bce7 241void AliTRDpidRefMakerNN::MakeTrainingLists(Int_t mombin)
1ee39b3a 242{
243 //
244 // build the training lists for the neural networks
245 //
246
9ec5bce7 247 if (!fRef) {
248 LoadFile(Form("TRD.Calib%s.root", GetName()));
1ee39b3a 249 }
250
9ec5bce7 251 if (!fRef) {
252 AliError("ERROR file for building training list not available");
1ee39b3a 253 return;
254 }
255
9ec5bce7 256 AliDebug(2, "\n Making training lists! \n");
1ee39b3a 257
9ec5bce7 258 Int_t nPart[AliPID::kSPECIES];
259 memset(nPart, 0, AliPID::kSPECIES*sizeof(Int_t));
1ee39b3a 260
261 // set needed branches
9ec5bce7 262 fTrainData[mombin] -> SetBranchAddress("fdEdx", fdEdx);
263 fTrainData[mombin] -> SetBranchAddress("fPID", fPID);
264 fTrainData[mombin] -> SetBranchAddress("fLy", &fLy);
265 fTrainData[mombin] -> SetBranchAddress("fNtrkl", &fNtrkl);
266
1ee39b3a 267 // start first loop to check total number of each particle type
9ec5bce7 268 for(Int_t iEv=0; iEv < fTrainData[mombin] -> GetEntries(); iEv++){
269 fTrainData[mombin] -> GetEntry(iEv);
1ee39b3a 270
271 // use only events with goes through 6 layers TRD
9ec5bce7 272 if(fNtrkl != AliTRDgeometry::kNlayer) continue;
273
274 if(fPID[AliPID::kElectron] == 1)
275 nPart[AliPID::kElectron]++;
276 else if(fPID[AliPID::kMuon] == 1)
277 nPart[AliPID::kMuon]++;
278 else if(fPID[AliPID::kPion] == 1)
279 nPart[AliPID::kPion]++;
280 else if(fPID[AliPID::kKaon] == 1)
281 nPart[AliPID::kKaon]++;
282 else if(fPID[AliPID::kProton] == 1)
283 nPart[AliPID::kProton]++;
1ee39b3a 284 }
285
286 AliDebug(2, "Particle multiplicities:");
9ec5bce7 287 AliDebug(2, Form("Momentum[%d] Elecs[%d] Muons[%d] Pions[%d] Kaons[%d] Protons[%d]", mombin, nPart[AliPID::kElectron], nPart[AliPID::kMuon], nPart[AliPID::kPion], nPart[AliPID::kKaon], nPart[AliPID::kProton]));
288
289
1ee39b3a 290
9ec5bce7 291 // // implement counter of training and test sample size
292 Int_t iTrain = 0, iTest = 0;
1ee39b3a 293
294 // set training sample size per momentum interval to 2/3
295 // of smallest particle counter and test sample to 1/3
9ec5bce7 296 iTrain = nPart[0];
297 for(Int_t iPart = 1; iPart < AliPID::kSPECIES; iPart++){
298 // exclude muons and kaons if not availyable
299 // this is neeeded since we do not have v0 candiates
300 if((iPart == AliPID::kMuon || iPart == AliPID::kKaon) && (nPart[AliPID::kMuon] == 0 || nPart[AliPID::kKaon] == 0)) continue;
301 if(iTrain > nPart[iPart])
302 iTrain = nPart[iPart];
303 }
304 iTest = Int_t( iTrain * (1-fFreq));
305 iTrain = Int_t(iTrain * fFreq);
306 AliDebug(2, Form("Momentum[%d] Train[%d] Test[%d]", mombin, iTrain, iTest));
307
1ee39b3a 308
309
310 // reset couters
9ec5bce7 311 memset(nPart, 0, AliPID::kSPECIES*sizeof(Int_t));
1ee39b3a 312
313 // start second loop to set the event lists
9ec5bce7 314 for(Int_t iEv = 0; iEv < fTrainData[mombin] -> GetEntries(); iEv++){
315 fTrainData[mombin] -> GetEntry(iEv);
316
317 // set event list
318 for(Int_t is = 0; is < AliPID::kSPECIES; is++){
319 if(nPart[is] < iTrain && fPID[is] == 1){
320 fTrain[mombin] -> Enter(iEv);
321 nPart[is]++;
322 } else if(nPart[is] < iTest+iTrain && fPID[is] == 1){
323 fTest[mombin] -> Enter(iEv);
324 nPart[is]++;
1ee39b3a 325 } else continue;
326 }
327 }
328
329 AliDebug(2, "Particle multiplicities in both lists:");
9ec5bce7 330 AliDebug(2, Form("Momentum[%d] Elecs[%d] Muons[%d] Pions[%d] Kaons[%d] Protons[%d]", mombin, nPart[AliPID::kElectron], nPart[AliPID::kMuon], nPart[AliPID::kPion], nPart[AliPID::kKaon], nPart[AliPID::kProton]));
331
332 return;
1ee39b3a 333}
334
335
336//________________________________________________________________________
337void AliTRDpidRefMakerNN::MakeRefs(Int_t mombin)
338{
339 //
340 // train the neural networks
341 //
342
343
9ec5bce7 344 if (!fTrainData[mombin]) LoadFile(Form("TRD.Calib%s.root", GetName()));
1ee39b3a 345
9ec5bce7 346 if (!fTrainData[mombin]) {
1ee39b3a 347 AliError("Tree for training list not available");
348 return;
349 }
350
351 TDatime datime;
352 fDate = datime.GetDate();
353
354 AliDebug(2, Form("Training momentum bin %d", mombin));
355
356 // set variable to monitor the training and to save the development of the networks
357 Int_t nEpochs = fEpochs/kMoniTrain;
358 AliDebug(2, Form("Training %d times %d epochs", kMoniTrain, nEpochs));
359
360 // make directories to save the networks
361 gSystem->Exec(Form("rm -r ./Networks_%d/MomBin_%d",fDate, mombin));
362 gSystem->Exec(Form("mkdir ./Networks_%d/MomBin_%d",fDate, mombin));
363
364 // variable to check if network can load weights from previous training
9ec5bce7 365 Bool_t bFirstLoop = kTRUE;
1ee39b3a 366
367 // train networks over several loops and save them after each loop
368 for(Int_t iLoop = 0; iLoop < kMoniTrain; iLoop++){
9ec5bce7 369 fTrainData[mombin] -> SetEventList(fTrain[mombin]);
370 fTrainData[mombin] -> SetEventList(fTest[mombin]);
1ee39b3a 371
9ec5bce7 372 AliDebug(2, Form("Momentum[%d] Trainingloop[%d]", mombin, iLoop));
1ee39b3a 373
9ec5bce7 374 // check if network is already implemented
375 if(bFirstLoop == kTRUE){
376 fNet = new TMultiLayerPerceptron("fdEdx[0],fdEdx[1],fdEdx[2],fdEdx[3],fdEdx[4],fdEdx[5],fdEdx[6],fdEdx[7]:15:7:fPID[0],fPID[1],fPID[2],fPID[3],fPID[4]!",fTrainData[mombin],fTrain[mombin],fTest[mombin]);
377 fNet -> SetLearningMethod(TMultiLayerPerceptron::kStochastic); // set learning method
378 fNet -> TMultiLayerPerceptron::SetEta(0.001); // set learning speed
379 if(!fContinueTraining){
380 if(AliLog::GetDebugLevel("","AliTRDpidRefMakerNN")>=2) fNet -> Train(nEpochs,"text update=10, graph");
381 else fNet -> Train(nEpochs,"");
1ee39b3a 382 }
9ec5bce7 383 else{
384 fNet -> LoadWeights(Form("./Networks_%d/MomBin_%d/Net_%d",fTrainPath, mombin, kMoniTrain - 1));
385 if(AliLog::GetDebugLevel("","AliTRDpidRefMakerNN")>=2) fNet -> Train(nEpochs,"text update=10, graph+");
386 else fNet -> Train(nEpochs,"+");
1ee39b3a 387 }
9ec5bce7 388 bFirstLoop = kFALSE;
389 }
390 else{
391 if(AliLog::GetDebugLevel("","AliTRDpidRefMakerNN")>=2) fNet -> Train(nEpochs,"text update=10, graph+");
392 else fNet -> Train(nEpochs,"+");
393 }
1ee39b3a 394
9ec5bce7 395 // save weights for monitoring of the training
396 fNet -> DumpWeights(Form("./Networks_%d/MomBin_%d/Net_%d",fDate, mombin, iLoop));
1ee39b3a 397 } // end training loop
398}
399
400
401
402//________________________________________________________________________
403void AliTRDpidRefMakerNN::MonitorTraining(Int_t mombin)
404{
405 //
406 // train the neural networks
407 //
408
9ec5bce7 409 if (!fTrainData[mombin]) LoadFile(Form("TRD.Calib%s.root", GetName()));
410 if (!fTrainData[mombin]) {
1ee39b3a 411 AliError("Tree for training list not available");
412 return;
413 }
414
415 // init networks and set event list
416 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
9ec5bce7 417 fNet = new TMultiLayerPerceptron("fdEdx[0],fdEdx[1],fdEdx[2],fdEdx[3],fdEdx[4],fdEdx[5],fdEdx[6],fdEdx[7]:15:7:fPID[0],fPID[1],fPID[2],fPID[3],fPID[4]!",fTrainData[mombin],fTrain[mombin],fTest[mombin]);
418 fTrainData[mombin] -> SetEventList(fTrain[mombin]);
419 fTrainData[mombin] -> SetEventList(fTest[mombin]);
1ee39b3a 420 }
421
422 // implement variables for likelihoods
423 Float_t like[AliPID::kSPECIES][AliTRDgeometry::kNlayer];
424 memset(like, 0, AliPID::kSPECIES*AliTRDgeometry::kNlayer*sizeof(Float_t));
425 Float_t likeAll[AliPID::kSPECIES], totProb;
426
427 Double_t pionEffiTrain[kMoniTrain], pionEffiErrTrain[kMoniTrain];
428 Double_t pionEffiTest[kMoniTrain], pionEffiErrTest[kMoniTrain];
429 memset(pionEffiTrain, 0, kMoniTrain*sizeof(Double_t));
430 memset(pionEffiErrTrain, 0, kMoniTrain*sizeof(Double_t));
431 memset(pionEffiTest, 0, kMoniTrain*sizeof(Double_t));
432 memset(pionEffiErrTest, 0, kMoniTrain*sizeof(Double_t));
433
434 // init histos
435 const Float_t epsilon = 1/(2*(AliTRDpidUtil::kBins-1)); // get nice histos with bin center at 0 and 1
436 TH1F *hElecs, *hPions;
437 hElecs = new TH1F("hElecs","Likelihood for electrons", AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon);
438 hPions = new TH1F("hPions","Likelihood for pions", AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon);
439
9ec5bce7 440 TGraphErrors *gEffisTrain = new TGraphErrors(kMoniTrain);
441 gEffisTrain -> SetLineColor(4);
442 gEffisTrain -> SetMarkerColor(4);
443 gEffisTrain -> SetMarkerStyle(29);
444 gEffisTrain -> SetMarkerSize(1);
445
446 TGraphErrors *gEffisTest = new TGraphErrors(kMoniTrain);
447 gEffisTest -> SetLineColor(2);
448 gEffisTest -> SetMarkerColor(2);
449 gEffisTest -> SetMarkerStyle(29);
450 gEffisTest -> SetMarkerSize(1);
1ee39b3a 451
452 AliTRDpidUtil *util = new AliTRDpidUtil();
453
454 // monitor training progress
455 for(Int_t iLoop = 0; iLoop < kMoniTrain; iLoop++){
456
457 // load weights
9ec5bce7 458 fNet -> LoadWeights(Form("./Networks_%d/MomBin_%d/Net_%d",fDate, mombin, iLoop));
459 AliDebug(2, Form("./Networks_%d/MomBin_%d/Net_%d",fDate, mombin, iLoop));
1ee39b3a 460
461 // event loop training list
1ee39b3a 462
9ec5bce7 463 for(Int_t is = 0; is < AliPID::kSPECIES; is++){
1ee39b3a 464
9ec5bce7 465 if(!((is == AliPID::kElectron) || (is == AliPID::kPion))) continue;
1ee39b3a 466
9ec5bce7 467 Int_t iChamb = 0;
468 // reset particle probabilities
1ee39b3a 469 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
9ec5bce7 470 likeAll[iPart] = 1./AliPID::kSPECIES;
1ee39b3a 471 }
9ec5bce7 472 totProb = 0.;
1ee39b3a 473
9ec5bce7 474 AliDebug(2, Form("%d",fTrain[mombin] -> GetN()));
475 for(Int_t iEvent = 0; iEvent < fTrain[mombin] -> GetN(); iEvent++ ){
476 fTrainData[mombin] -> GetEntry(fTrain[mombin] -> GetEntry(iEvent));
477 // use event only if it is electron or pion
478 if(!(fPID[is] == 1.0)) continue;
479 // get the probabilities for each particle type in each chamber
480 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
481 like[iPart][iChamb] = fNet -> Result(fTrain[mombin] -> GetEntry(iEvent), iPart);
482 likeAll[iPart] *= like[iPart][iChamb];
483 }
484 //end chamber loop
485 iChamb++;
486
487 // get total probability and normalize it
488 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
489 totProb += likeAll[iPart];
490 }
491 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
492 likeAll[iPart] /= totProb;
493 }
494
495 if(iChamb == 5){
496 // fill likelihood distributions
497 if(fPID[AliPID::kElectron] == 1)
498 hElecs -> Fill(likeAll[AliPID::kElectron]);
499 if(fPID[AliPID::kPion] == 1)
500 hPions -> Fill(likeAll[AliPID::kElectron]);
501 iChamb = 0;
502 // reset particle probabilities
503 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
504 likeAll[iPart] = 1./AliPID::kSPECIES;
505 }
506 totProb = 0.;
507 }
508 } // end event loop
509 } // end species loop
510
1ee39b3a 511
512 // calculate the pion efficiency and fill the graph
513 util -> CalculatePionEffi(hElecs, hPions);
514 pionEffiTrain[iLoop] = util -> GetPionEfficiency();
515 pionEffiErrTrain[iLoop] = util -> GetError();
516
517 gEffisTrain -> SetPoint(iLoop, iLoop+1, pionEffiTrain[iLoop]);
518 gEffisTrain -> SetPointError(iLoop, 0, pionEffiErrTrain[iLoop]);
519 hElecs -> Reset();
520 hPions -> Reset();
521 AliDebug(2, Form("TrainingLoop[%d] PionEfficiency[%f +/- %f]", iLoop, pionEffiTrain[iLoop], pionEffiErrTrain[iLoop]));
522 // end training loop
523
524
525
9ec5bce7 526 // monitor validation progress
527 for(Int_t is = 0; is < AliPID::kSPECIES; is++){
1ee39b3a 528
9ec5bce7 529 if(!((is == AliPID::kElectron) || (is == AliPID::kPion))) continue;
530
531 Int_t iChamb = 0;
1ee39b3a 532 // reset particle probabilities
533 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
9ec5bce7 534 likeAll[iPart] = 1./AliPID::kSPECIES;
1ee39b3a 535 }
536 totProb = 0.;
537
9ec5bce7 538 for(Int_t iEvent = 0; iEvent < fTest[mombin] -> GetN(); iEvent++ ){
539 fTrainData[mombin] -> GetEntry(fTest[mombin] -> GetEntry(iEvent));
540 // use event only if it is electron or pion
541 if(!(fPID[is] == 1.0)) continue;
542
543 // get the probabilities for each particle type in each chamber
544 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
545 like[iPart][iChamb] = fNet -> Result(fTest[mombin] -> GetEntry(iEvent), iPart);
546 likeAll[iPart] *= like[iPart][iChamb];
547 }
548 iChamb++;
549
550 // get total probability and normalize it
551 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
552 totProb += likeAll[iPart];
553 }
554 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
555 likeAll[iPart] /= totProb;
556 }
557
558 if(iChamb == 5){
559 // fill likelihood distributions
560 if(fPID[AliPID::kElectron] == 1)
561 hElecs -> Fill(likeAll[AliPID::kElectron]);
562 if(fPID[AliPID::kPion] == 1)
563 hPions -> Fill(likeAll[AliPID::kElectron]);
564 iChamb = 0;
565 // reset particle probabilities
566 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
567 likeAll[iPart] = 1./AliPID::kSPECIES;
568 }
569 totProb = 0.;
570 }
571 } // end event loop
572 } // end species loop
573
1ee39b3a 574 // calculate the pion efficiency and fill the graph
575 util -> CalculatePionEffi(hElecs, hPions);
576 pionEffiTest[iLoop] = util -> GetPionEfficiency();
577 pionEffiErrTest[iLoop] = util -> GetError();
578
579 gEffisTest -> SetPoint(iLoop, iLoop+1, pionEffiTest[iLoop]);
580 gEffisTest -> SetPointError(iLoop, 0, pionEffiErrTest[iLoop]);
581 hElecs -> Reset();
582 hPions -> Reset();
583 AliDebug(2, Form("TestLoop[%d] PionEfficiency[%f +/- %f] \n", iLoop, pionEffiTest[iLoop], pionEffiErrTest[iLoop]));
584
9ec5bce7 585 } // end validation loop
1ee39b3a 586
587 util -> Delete();
588
589 gEffisTest -> Draw("PAL");
590 gEffisTrain -> Draw("PL");
591
592}
593
594
595//________________________________________________________________________
9ec5bce7 596Bool_t AliTRDpidRefMakerNN::LoadFile(const Char_t *InFileNN)
1ee39b3a 597{
598 //
599 // Loads the files and sets the event list
9ec5bce7 600 // for neural network training.
1ee39b3a 601 // Useable for training outside of the makeResults.C macro
602 //
603
9ec5bce7 604 fRef = TFile::Open(Form("%s", InFileNN));
605 if(!fRef) return 0;
606 for(Int_t ip = 0; ip < AliTRDCalPID::kNMom; ip++){
607 fTrainData[ip] = (TTree*)fRef -> Get(Form("fTrainData_%d", ip));
608 }
1ee39b3a 609
610 for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
9ec5bce7 611 fTrain[iMom] = new TEventList(Form("fTrainMom%d", iMom), Form("Training list for momentum intervall %d", iMom));
612 fTest[iMom] = new TEventList(Form("fTestMom%d", iMom), Form("Test list for momentum intervall %d", iMom));
1ee39b3a 613 }
9ec5bce7 614 return 1;
1ee39b3a 615}
616
617