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