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 "AliTRDReconstructor.h"
  41 #include "AliTRDpidUtil.h"
  42 #include "AliTRDpidRefMakerNN.h"
  43 #include "AliTRDpidUtil.h"
  44
  45 #include "Cal/AliTRDCalPID.h"
  46 #include "Cal/AliTRDCalPIDNN.h"
  47 #include "info/AliTRDtrackInfo.h"
  48 #include "info/AliTRDv0Info.h"
  49
  50 ClassImp(AliTRDpidRefMakerNN)
  51
  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)
  67 {
  68   //
  69   // Default constructor
  70   //
  71
  72   memset(fTrain, 0, AliTRDCalPID::kNMom*sizeof(TEventList*));
  73   memset(fTest, 0, AliTRDCalPID::kNMom*sizeof(TEventList*));
  74   memset(fTrainData, 0, AliTRDCalPID::kNMom*sizeof(TTree*));
  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 //________________________________________________________________________
  87 AliTRDpidRefMakerNN::~AliTRDpidRefMakerNN()
  88 {
  89 }
  90
  91
  92 //________________________________________________________________________
  93 void AliTRDpidRefMakerNN::CreateOutputObjects()
  94 {
  95   // Create output file and tree
  96   // Called once
  97
  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   }
 109
 110 }
 111
 112
 113
 114 //________________________________________________________________________
 115 Bool_t AliTRDpidRefMakerNN::PostProcess()
 116 {
 117   // Draw result to the screen
 118   // Called once at the end of the query
 119
 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   }
 189   if(!fDoTraining) return kTRUE;
 190
 191
 192
 193   // train the neural networks
 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)){
 199     if(fTrain[fTrainMomBin] -> GetN() < fMinTrain){
 200       AliError("Warning in AliTRDpidRefMakerNN::PostProcess : Not enough events for training available! Please check Data sample!");
 201       return kFALSE;
 202     }
 203     MakeRefs(fTrainMomBin);
 204     MonitorTraining(fTrainMomBin);
 205   }
 206   // train all momenta
 207   else{
 208     for(Int_t iMomBin = 0; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
 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;
 212       }
 213       MakeRefs(fTrainMomBin);
 214       MonitorTraining(iMomBin);
 215     }
 216   }
 217
 218   return kTRUE; // testing protection
 219 }
 220
 221
 222 //________________________________________________________________________
 223 void AliTRDpidRefMakerNN::MakeTrainingLists(Int_t mombin)
 224 {
 225   //
 226   // build the training lists for the neural networks
 227   //
 228
 229   if (!fRef) {
 230     LoadFile(Form("TRD.Calib%s.root", GetName()));
 231   }
 232
 233   if (!fRef) {
 234     AliError("ERROR file for building training list not available");
 235     return;
 236   }
 237
 238   AliDebug(2, "\n Making training lists! \n");
 239
 240   Int_t nPart[AliPID::kSPECIES];
 241   memset(nPart, 0, AliPID::kSPECIES*sizeof(Int_t));
 242
 243   // set needed branches
 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
 249   // start first loop to check total number of each particle type
 250   for(Int_t iEv=0; iEv < fTrainData[mombin] -> GetEntries(); iEv++){
 251     fTrainData[mombin] -> GetEntry(iEv);
 252
 253     // use only events with goes through 6 layers TRD
 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]++;
 266   }
 267
 268   AliDebug(2, "Particle multiplicities:");
 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
 272
 273   //   // implement counter of training and test sample size
 274   Int_t iTrain = 0, iTest = 0;
 275
 276   // set training sample size per momentum interval to 2/3
 277   // of smallest particle counter and test sample to 1/3
 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
 290
 291
 292   // reset couters
 293   memset(nPart, 0, AliPID::kSPECIES*sizeof(Int_t));
 294
 295   // start second loop to set the event lists
 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]++;
 307       } else continue;
 308     }
 309   }
 310
 311   AliDebug(2, "Particle multiplicities in both lists:");
 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;
 315 }
 316
 317
 318 //________________________________________________________________________
 319 void AliTRDpidRefMakerNN::MakeRefs(Int_t mombin)
 320 {
 321   //
 322   // train the neural networks
 323   //
 324
 325
 326   if (!fTrainData[mombin]) LoadFile(Form("TRD.Calib%s.root", GetName()));
 327
 328   if (!fTrainData[mombin]) {
 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
 347   Bool_t bFirstLoop = kTRUE;
 348
 349   // train networks over several loops and save them after each loop
 350   for(Int_t iLoop = 0; iLoop < kMoniTrain; iLoop++){
 351     fTrainData[mombin] -> SetEventList(fTrain[mombin]);
 352     fTrainData[mombin] -> SetEventList(fTest[mombin]);
 353
 354     AliDebug(2, Form("Momentum[%d] Trainingloop[%d]", mombin, iLoop));
 355
 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,"");
 364       }
 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,"+");
 369       }
 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     }
 376
 377     // save weights for monitoring of the training
 378     fNet -> DumpWeights(Form("./Networks_%d/MomBin_%d/Net_%d",fDate, mombin, iLoop));
 379   }   // end training loop
 380 }
 381
 382
 383
 384 //________________________________________________________________________
 385 void AliTRDpidRefMakerNN::MonitorTraining(Int_t mombin)
 386 {
 387   //
 388   // train the neural networks
 389   //
 390
 391   if (!fTrainData[mombin]) LoadFile(Form("TRD.Calib%s.root", GetName()));
 392   if (!fTrainData[mombin]) {
 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++){
 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]);
 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
 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);
 433
 434   AliTRDpidUtil *util = new AliTRDpidUtil();
 435
 436   // monitor training progress
 437   for(Int_t iLoop = 0; iLoop < kMoniTrain; iLoop++){
 438
 439     // load weights
 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));
 442
 443     // event loop training list
 444
 445     for(Int_t is = 0; is < AliPID::kSPECIES; is++){
 446
 447       if(!((is == AliPID::kElectron) || (is == AliPID::kPion))) continue;
 448
 449       Int_t iChamb = 0;
 450       // reset particle probabilities
 451       for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
 452         likeAll[iPart] = 1./AliPID::kSPECIES;
 453       }
 454       totProb = 0.;
 455
 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
 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
 508     // monitor validation progress
 509     for(Int_t is = 0; is < AliPID::kSPECIES; is++){
 510
 511       if(!((is == AliPID::kElectron) || (is == AliPID::kPion))) continue;
 512
 513       Int_t iChamb = 0;
 514       // reset particle probabilities
 515       for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
 516         likeAll[iPart] = 1./AliPID::kSPECIES;
 517       }
 518       totProb = 0.;
 519
 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
 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
 567   } //   end validation loop
 568
 569   util -> Delete();
 570
 571   gEffisTest -> Draw("PAL");
 572   gEffisTrain -> Draw("PL");
 573
 574 }
 575
 576
 577 //________________________________________________________________________
 578 Bool_t AliTRDpidRefMakerNN::LoadFile(const Char_t *InFileNN)
 579 {
 580   //
 581   // Loads the files and sets the event list
 582   // for neural network training.
 583   // Useable for training outside of the makeResults.C macro
 584   //
 585
 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   }
 591
 592   for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
 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));
 595   }
 596   return 1;
 597 }
 598
 599