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