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