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