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