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