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