e THnSparse structure to store MC residuals
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDpidRefMakerNN.cxx
1 /*************************************************************************
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  **************************************************************************/
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 "AliTRDpidUtil.h"
41 #include "AliTRDpidRefMakerNN.h"
42 #include "AliTRDpidUtil.h"
43
44 #include "Cal/AliTRDCalPID.h"
45 #include "Cal/AliTRDCalPIDNN.h"
46 #include "info/AliTRDtrackInfo.h"
47 #include "info/AliTRDv0Info.h"
48 #include "info/AliTRDpidInfo.h"
49
50 ClassImp(AliTRDpidRefMakerNN)
51
52 //________________________________________________________________________
53   AliTRDpidRefMakerNN::AliTRDpidRefMakerNN() 
54   :AliTRDpidRefMaker()
55   ,fNet(NULL)
56   ,fTrainMomBin(kAll)
57   ,fEpochs(1000)
58   ,fMinTrain(100)
59   ,fDate(0)
60   ,fDoTraining(0)
61   ,fContinueTraining(0)
62   ,fTrainPath(0)
63   ,fScale(0)
64   ,fLy(0)
65   ,fNtrkl(0)
66   ,fRef(NULL)
67 {
68   //
69   // Default constructor
70   //
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();
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)
93   ,fTrainPath(0)
94   ,fScale(0)
95   ,fLy(0)
96   ,fNtrkl(0)
97   ,fRef(NULL)
98 {
99   //
100   // Default constructor
101   //
102
103   memset(fTrain, 0, AliTRDCalPID::kNMom*sizeof(TEventList*));
104   memset(fTest, 0, AliTRDCalPID::kNMom*sizeof(TEventList*));
105   memset(fTrainData, 0, AliTRDCalPID::kNMom*sizeof(TTree*));
106
107   SetAbundance(.67);
108   SetScaledEdx(Float_t(AliTRDCalPIDNN::kMLPscale));
109   TDatime datime;
110   fDate = datime.GetDate();
111 }
112
113
114 //________________________________________________________________________
115 AliTRDpidRefMakerNN::~AliTRDpidRefMakerNN() 
116 {
117 }
118
119
120 //________________________________________________________________________
121 void AliTRDpidRefMakerNN::MakeTrainTestTrees()
122 {
123   // Create output file and tree
124   // Called once
125
126   fRef = new TFile("TRD.CalibPIDrefMakerNN.root", "RECREATE");
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   }
137 }
138
139
140
141 //________________________________________________________________________
142 Bool_t AliTRDpidRefMakerNN::PostProcess()
143 {
144   // Draw result to the screen
145   // Called once at the end of the query
146
147   TFile *fCalib = TFile::Open(Form("AnalysisResults.root"));
148   if (!fCalib) {
149     AliError("Calibration file not available");
150     return kFALSE;
151   }
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");
158   if (!fData) {
159     AliError("Calibration data not available");
160     return kFALSE;
161   }
162   TObjArray *o = NULL;
163   if(!(o = (TObjArray*)dCalib->Get("MonitorNN"))) {
164     AliWarning("Missing monitoring container.");
165     return kFALSE;
166   }
167   fContainer = (TObjArray*)o->Clone("monitor");
168
169
170   
171   if (!fRef) {
172     AliDebug(2, "Loading file TRD.CalibPIDrefMakerNN.root");
173     LoadFile("TRD.CalibPIDrefMakerNN.root");
174   }
175   else AliDebug(2, "file available");
176
177   if (!fRef) {
178     MakeTrainingSample();
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   }
187   if(!fDoTraining) return kTRUE;
188
189
190
191   // train the neural networks
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)){
197     if(fTrain[fTrainMomBin] -> GetN() < fMinTrain){
198       AliError("Warning in AliTRDpidRefMakerNN::PostProcess : Not enough events for training available! Please check Data sample!");
199       return kFALSE;
200     }
201     MakeRefs(fTrainMomBin);
202     MonitorTraining(fTrainMomBin);
203   }
204   // train all momenta
205   else{
206     for(Int_t iMomBin = 0; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
207       if(fTrain[iMomBin] -> GetN() < fMinTrain){
208         AliError(Form("Warning in AliTRDpidRefMakerNN::PostProcess : Not enough events for training available for momentum bin [%d]! Please check Data sample!", iMomBin));
209         continue;
210       }
211       MakeRefs(iMomBin);
212       MonitorTraining(iMomBin);
213     }
214   }
215
216   return kTRUE; // testing protection
217 }
218
219
220 //________________________________________________________________________
221 Bool_t AliTRDpidRefMakerNN::MakeTrainingSample() 
222 {
223   // convert AnalysisResults.root to training file
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
290 //________________________________________________________________________
291 void AliTRDpidRefMakerNN::MakeTrainingLists(Int_t mombin) 
292 {
293   //
294   // build the training lists for the neural networks
295   //
296
297   if (!fRef) {
298     LoadFile(Form("TRD.Calib%s.root", GetName()));
299   }
300
301   if (!fRef) {
302     AliError("ERROR file for building training list not available");
303     return;
304   }
305
306   AliDebug(2, "\n Making training lists! \n");
307
308   Int_t nPart[AliPID::kSPECIES];
309   memset(nPart, 0, AliPID::kSPECIES*sizeof(Int_t));
310
311   // set needed branches
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   
317   // start first loop to check total number of each particle type
318   for(Int_t iEv=0; iEv < fTrainData[mombin] -> GetEntries(); iEv++){
319     fTrainData[mombin] -> GetEntry(iEv);
320
321     // use only events with goes through 6 layers TRD
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]++;
334   }
335
336   AliDebug(2, "Particle multiplicities:");
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   
340
341   //   // implement counter of training and test sample size
342   Int_t iTrain = 0, iTest = 0;
343
344   // set training sample size per momentum interval to 2/3 
345   // of smallest particle counter and test sample to 1/3
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
358
359
360   // reset couters
361   memset(nPart, 0, AliPID::kSPECIES*sizeof(Int_t));
362
363   // start second loop to set the event lists
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]++;
375       } else continue;
376     }
377   }
378   
379   AliDebug(2, "Particle multiplicities in both lists:");
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;
383 }
384
385
386 //________________________________________________________________________
387 void AliTRDpidRefMakerNN::MakeRefs(Int_t mombin) 
388 {
389   //
390   // train the neural networks
391   //
392   
393   
394   if (!fTrainData[mombin]) LoadFile(Form("TRD.CalibPIDrefMakerNN.root"));
395
396   if (!fTrainData[mombin]) {
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
415   Bool_t bFirstLoop = kTRUE;
416
417   // train networks over several loops and save them after each loop
418   for(Int_t iLoop = 0; iLoop < kMoniTrain; iLoop++){
419     fTrainData[mombin] -> SetEventList(fTrain[mombin]);
420     fTrainData[mombin] -> SetEventList(fTest[mombin]);
421       
422     AliDebug(2, Form("Momentum[%d] Trainingloop[%d]", mombin, iLoop));
423       
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,"");
432       }
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,"+");                   
437       }
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     }
444       
445     // save weights for monitoring of the training
446     fNet -> DumpWeights(Form("./Networks_%d/MomBin_%d/Net_%d",fDate, mombin, iLoop));
447   }   // end training loop
448 }
449
450
451
452 //________________________________________________________________________
453 void AliTRDpidRefMakerNN::MonitorTraining(Int_t mombin) 
454 {
455   //
456   // train the neural networks
457   //
458   
459   if (!fTrainData[mombin]) LoadFile(Form("TRD.CalibPIDrefMakerNN.root"));
460   if (!fTrainData[mombin]) {
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++){
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]);
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
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);
501
502   AliTRDpidUtil *util = new AliTRDpidUtil();
503   
504   // monitor training progress
505   for(Int_t iLoop = 0; iLoop < kMoniTrain; iLoop++){
506
507     // load weights
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));
510
511     // event loop training list
512
513     for(Int_t is = 0; is < AliPID::kSPECIES; is++){
514
515       if(!((is == AliPID::kElectron) || (is == AliPID::kPion))) continue;
516
517       Int_t iChamb = 0;
518       // reset particle probabilities
519       for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
520         likeAll[iPart] = 1./AliPID::kSPECIES;
521       }
522       totProb = 0.;
523
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
528         if(fPID[is] <1.e-5) continue;
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     
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
576     // monitor validation progress
577     for(Int_t is = 0; is < AliPID::kSPECIES; is++){
578
579       if(!((is == AliPID::kElectron) || (is == AliPID::kPion))) continue;
580
581       Int_t iChamb = 0;
582       // reset particle probabilities
583       for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
584         likeAll[iPart] = 1./AliPID::kSPECIES;
585       }
586       totProb = 0.;
587
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
591         if(TMath::Abs(fPID[is]- 1.)<1.e-5) continue;
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     
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     
635   } //   end validation loop
636
637   util -> Delete();
638
639   gEffisTest -> Draw("PAL");
640   gEffisTrain -> Draw("PL");
641
642 }
643
644
645 //________________________________________________________________________
646 Bool_t AliTRDpidRefMakerNN::LoadFile(const Char_t *InFileNN) 
647 {
648   //
649   // Loads the files and sets the event list
650   // for neural network training.
651   // Useable for training outside of the makeResults.C macro
652   //
653
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   }
659
660   for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
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));
663   }
664   return 1;
665 }
666
667