TRD/qaRec/AliTRDpidRefMaker.cxx

   1 #include "TSystem.h"
   2 #include "TDatime.h"
   3 #include "TPDGCode.h"
   4 #include "TH1F.h"
   5 #include "TFile.h"
   6 #include "TTree.h"
   7 #include "TTreeStream.h"
   8 #include "TEventList.h"
   9 #include "TMultiLayerPerceptron.h"
  10
  11 #include "AliPID.h"
  12 #include "AliESDEvent.h"
  13 #include "AliESDInputHandler.h"
  14 #include "AliTrackReference.h"
  15
  16 #include "AliAnalysisTask.h"
  17
  18 #include "AliTRDtrackV1.h"
  19 #include "AliTRDReconstructor.h"
  20 #include "../Cal/AliTRDCalPID.h"
  21 #include "../Cal/AliTRDCalPIDNN.h"
  22
  23 #include "AliTRDpidRefMaker.h"
  24 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
  25
  26 #define SETFLG(n,f) ((n) |= f)
  27 #define CLRFLG(n,f) ((n) &= ~f)
  28
  29 // builds the reference tree for the training of neural networks
  30
  31
  32 ClassImp(AliTRDpidRefMaker)
  33
  34 //________________________________________________________________________
  35 AliTRDpidRefMaker::AliTRDpidRefMaker()
  36   :AliTRDrecoTask("PidRefMaker", "PID(NN) Reference Maker")
  37   ,fReconstructor(0x0)
  38   ,fNN(0x0)
  39   ,fLQ(0x0)
  40   ,fLayer(0xff)
  41   ,fTrainMomBin(-1)
  42   ,fEpochs(0)
  43   ,fMinTrain(0)
  44   ,fDate(0)
  45   ,fMom(0.)
  46 {
  47   //
  48   // Default constructor
  49   //
  50
  51   fReconstructor = new AliTRDReconstructor();
  52   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
  53   memset(fv0pid, 0, AliPID::kSPECIES*sizeof(Float_t));
  54   memset(fdEdx, 0, 10*sizeof(Float_t));
  55
  56   const Int_t nnSize = AliTRDCalPID::kNMom * AliTRDCalPID::kNPlane;
  57   memset(fTrain, 0, nnSize*sizeof(TEventList*));
  58   memset(fTest, 0, nnSize*sizeof(TEventList*));
  59   memset(fNet, 0, AliTRDCalPID::kNPlane*sizeof(TMultiLayerPerceptron*));
  60
  61   DefineOutput(1, TTree::Class());
  62   DefineOutput(2, TTree::Class());
  63 }
  64
  65
  66 //________________________________________________________________________
  67 AliTRDpidRefMaker::AliTRDpidRefMaker(const Char_t *InFileNN, const Char_t *InFileLQ)
  68   :AliTRDrecoTask("PidRefMaker", "PID(NN) Reference Maker")
  69   ,fReconstructor(0x0)
  70   ,fNN(0x0)
  71   ,fLQ(0x0)
  72   ,fLayer(0xff)
  73   ,fTrainMomBin(-1)
  74   ,fEpochs(0)
  75   ,fMinTrain(0)
  76   ,fDate(0)
  77   ,fMom(0.)
  78 {
  79   //
  80   // constructor for the makeResults macro
  81   //
  82
  83   fReconstructor = new AliTRDReconstructor();
  84   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
  85   memset(fv0pid, 0, AliPID::kSPECIES*sizeof(Float_t));
  86   memset(fdEdx, 0, 10*sizeof(Float_t));
  87
  88   TFile *fInFileNN;
  89   fInFileNN = new TFile(InFileNN, "READ");
  90   fNN = (TTree*)fInFileNN -> Get("NN");
  91
  92   TFile *fInFileLQ;
  93   fInFileLQ = new TFile(InFileLQ, "READ");
  94   fLQ = (TTree*)fInFileLQ -> Get("LQ");
  95
  96   for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
  97     for(Int_t ily = 0; ily < AliTRDCalPID::kNPlane; ily++){
  98       fTrain[iMom][ily] = new TEventList(Form("fTrainMom%d_%d", iMom, ily), Form("Training list for momentum intervall %d and plane %d", iMom, ily));
  99       fTest[iMom][ily] = new TEventList(Form("fTestMom%d_%d", iMom, ily), Form("Test list for momentum intervall %d and plane %d", iMom, ily));
 100     }
 101   }
 102
 103   TDatime datime;
 104   fDate = datime.GetDate();
 105
 106   SetDebugLevel(2);
 107   SetEpochs(1000);
 108   SetTrainMomBin(kAll);
 109   SetMinTrain(1000);
 110
 111   if(fDebugLevel>=2) Printf("Epochs[%d] MomBin[%d] MinTrain[%d]", fEpochs, fTrainMomBin, fMinTrain);
 112
 113 //   DefineOutput(1, TTree::Class());
 114 //   DefineOutput(2, TTree::Class());
 115 }
 116
 117
 118 //________________________________________________________________________
 119 AliTRDpidRefMaker::~AliTRDpidRefMaker()
 120 {
 121   if(fReconstructor) delete fReconstructor;
 122 }
 123
 124
 125 //________________________________________________________________________
 126 void AliTRDpidRefMaker::CreateOutputObjects()
 127 {
 128   // Create histograms
 129   // Called once
 130
 131   OpenFile(0, "RECREATE");
 132   fContainer = new TObjArray();
 133   fContainer->AddAt(new TH1F("hPDG","hPDG",AliPID::kSPECIES,-0.5,5.5),0);
 134
 135   // open reference TTree for NN
 136   OpenFile(1, "RECREATE");
 137   fNN = new TTree("NN", "Reference data for NN");
 138   fNN->Branch("fLayer", &fLayer, "fLayer/I");
 139   fNN->Branch("fMom", &fMom, "fMom/F");
 140   fNN->Branch("fv0pid", fv0pid, Form("fv0pid[%d]/F", AliPID::kSPECIES));
 141   fNN->Branch("fdEdx", fdEdx, Form("fdEdx[%d]/F", AliTRDReconstructor::kNNslices));
 142
 143   // open reference TTree for LQ
 144   OpenFile(2, "RECREATE");
 145   fLQ = new TTree("LQ", "Reference data for LQ");
 146   fLQ->Branch("fLayer", &fLayer, "fLayer/I");
 147   fLQ->Branch("fMom", &fMom, "fMom/F");
 148   fLQ->Branch("fv0pid", fv0pid, Form("fv0pid[%d]/F", AliPID::kSPECIES));
 149   fLQ->Branch("fdEdx", fdEdx, Form("fdEdx[%d]/F", AliTRDReconstructor::kLQslices));
 150 }
 151
 152
 153 //________________________________________________________________________
 154 void AliTRDpidRefMaker::Exec(Option_t *)
 155 {
 156   // Main loop
 157   // Called for each event
 158
 159   Int_t labelsacc[10000];
 160   memset(labelsacc, 0, sizeof(Int_t) * 10000);
 161
 162   Float_t mom;
 163   ULong_t status;
 164   Int_t nTRD = 0;
 165
 166   AliTRDtrackInfo     *track = 0x0;
 167   AliTRDtrackV1    *TRDtrack = 0x0;
 168   AliTrackReference     *ref = 0x0;
 169   AliExternalTrackParam *esd = 0x0;
 170
 171   AliTRDseedV1 *TRDtracklet = 0x0;
 172
 173   //AliTRDcluster *TRDcluster = 0x0;
 174
 175   for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
 176
 177     // reset the pid information
 178     for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++)
 179       fv0pid[iPart] = 0.;
 180
 181     track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
 182     if(!track->HasESDtrack()) continue;
 183     status = track->GetStatus();
 184     if(!(status&AliESDtrack::kTPCout)) continue;
 185
 186     if(!(TRDtrack = track->GetTRDtrack())) continue;
 187     //&&(track->GetNumberOfClustersRefit()
 188
 189     // use only tracks that hit 6 chambers
 190     if(!(TRDtrack->GetNumberOfTracklets() == AliTRDCalPID::kNPlane)) continue;
 191
 192     ref = track->GetTrackRef(0);
 193     esd = track->GetOuterParam();
 194     mom = ref ? ref->P(): esd->P();
 195     fMom = mom;
 196
 197
 198     labelsacc[nTRD] = track->GetLabel();
 199     nTRD++;
 200
 201     // if no monte carlo data available -> use V0 information
 202     if(!HasMCdata()){
 203       GetV0info(TRDtrack,fv0pid);
 204     }
 205     // else use the MC info
 206     else{
 207       switch(track -> GetPDG()){
 208       case kElectron:
 209       case kPositron:
 210         fv0pid[AliPID::kElectron] = 1.;
 211         break;
 212       case kMuonPlus:
 213       case kMuonMinus:
 214         fv0pid[AliPID::kMuon] = 1.;
 215         break;
 216       case kPiPlus:
 217       case kPiMinus:
 218         fv0pid[AliPID::kPion] = 1.;
 219         break;
 220       case kKPlus:
 221       case kKMinus:
 222         fv0pid[AliPID::kKaon] = 1.;
 223         break;
 224       case kProton:
 225       case kProtonBar:
 226         fv0pid[AliPID::kProton] = 1.;
 227         break;
 228       }
 229     }
 230
 231     // set reconstructor
 232     Float_t *dedx;
 233     TRDtrack -> SetReconstructor(fReconstructor);
 234
 235     // fill the dE/dx information for NN
 236     fReconstructor -> SetOption("nn");
 237     for(Int_t ily = 0; ily < AliTRDCalPID::kNPlane; ily++){
 238       if(!(TRDtracklet = TRDtrack -> GetTracklet(ily))) continue;
 239       TRDtracklet->CookdEdx(AliTRDReconstructor::kNNslices);
 240       dedx = TRDtracklet->GetdEdx();
 241       for(Int_t iSlice = 0; iSlice < AliTRDReconstructor::kNNslices; iSlice++)
 242         dedx[iSlice] = dedx[iSlice]/AliTRDCalPIDNN::kMLPscale;
 243       memcpy(fdEdx, dedx, AliTRDReconstructor::kNNslices*sizeof(Float_t));
 244       if(fDebugLevel>=2) Printf("LayerNN : %d", ily);
 245       fLayer = ily;
 246       fNN->Fill();
 247     }
 248
 249
 250     // fill the dE/dx information for LQ
 251     fReconstructor -> SetOption("!nn");
 252     for(Int_t ily = 0; ily < AliTRDCalPID::kNPlane; ily++){
 253       if(!(TRDtracklet = TRDtrack -> GetTracklet(ily))) continue;
 254       TRDtracklet->CookdEdx(AliTRDReconstructor::kLQslices);
 255       dedx = TRDtracklet->GetdEdx();
 256       memcpy(fdEdx, dedx, AliTRDReconstructor::kLQslices*sizeof(Float_t));
 257       if(fDebugLevel>=2) Printf("LayerLQ : %d", ily);
 258       fLayer = ily;
 259       fLQ->Fill();
 260     }
 261
 262
 263     for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
 264       if(fDebugLevel>=4) Printf("PDG is %d %f", iPart, fv0pid[iPart]);
 265     }
 266   }
 267
 268   PostData(0, fContainer);
 269   PostData(1, fNN);
 270   PostData(2, fLQ);
 271 }
 272
 273
 274 //________________________________________________________
 275 void AliTRDpidRefMaker::GetRefFigure(Int_t /*ifig*/, Int_t &/*first*/, Int_t &/*last*/, Option_t */*opt*/)
 276 {
 277
 278 }
 279
 280
 281 //________________________________________________________________________
 282 Bool_t AliTRDpidRefMaker::PostProcess()
 283 {
 284   // Draw result to the screen
 285   // Called once at the end of the query
 286
 287   // build the training andthe test list for the neural networks
 288   MakeTrainingLists();
 289
 290
 291   // train the neural networks and build the refrence histos for 2-dim LQ
 292   gSystem->Exec(Form("mkdir ./Networks_%d/",fDate));
 293   if(fDebugLevel>=2) Printf("TrainMomBin [%d] [%d]", fTrainMomBin, kAll);
 294
 295   // train single network for a single momentum (recommended)
 296   if(!(fTrainMomBin == kAll)){
 297     if(fTrain[fTrainMomBin][0] -> GetN() < fMinTrain){
 298       if(fDebugLevel>=2) Printf("Warning in AliTRDpidRefMaker::PostProcess : Not enough events for training available! Please check Data sample!");
 299       return kFALSE;
 300     }
 301     TrainNetworks(fTrainMomBin);
 302     BuildLQRefs(fTrainMomBin);
 303   }
 304   // train all momenta
 305   else{
 306     for(Int_t iMomBin = 0; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
 307       if(fTrain[iMomBin][0] -> GetN() < fMinTrain){
 308         if(fDebugLevel>=2) Printf("Warning in AliTRDpidRefMaker::PostProcess : Not enough events for training available for momentum bin [%d]! Please check Data sample!", iMomBin);
 309         continue;
 310       }
 311       TrainNetworks(iMomBin);
 312       BuildLQRefs(fTrainMomBin);
 313     }
 314   }
 315
 316   // monitor training progress
 317   for(Int_t iMomBin = 0; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
 318     MonitorTraining(iMomBin);
 319   }
 320
 321   return kTRUE; // testing protection
 322 }
 323
 324
 325 //________________________________________________________________________
 326 void AliTRDpidRefMaker::Terminate(Option_t *)
 327 {
 328   // Draw result to the screen
 329   // Called once at the end of the query
 330
 331   fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
 332   if (!fContainer) {
 333     Printf("ERROR: list not available");
 334     return;
 335   }
 336 }
 337
 338
 339 //________________________________________________________________________
 340 void AliTRDpidRefMaker::GetV0info(AliTRDtrackV1 *TRDtrack, Float_t *v0pid)
 341 {
 342   // !!!! PREMILMINARY FUNCTION !!!!
 343   //
 344   // this is the place for the V0 procedure
 345   // as long as there is no one implemented,
 346   // just the probabilities
 347   // of the TRDtrack are used!
 348
 349   TRDtrack -> SetReconstructor(fReconstructor);
 350   fReconstructor -> SetOption("nn");
 351   TRDtrack -> CookPID();
 352   for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
 353     v0pid[iPart] = TRDtrack -> GetPID(iPart);
 354     if(fDebugLevel>=4) Printf("PDG is (in V0info) %d %f", iPart, v0pid[iPart]);
 355   }
 356 }
 357
 358
 359 //________________________________________________________________________
 360 void AliTRDpidRefMaker::MakeTrainingLists()
 361 {
 362   //
 363   // build the training lists for the neural networks
 364   //
 365
 366   SetDebugLevel(2);
 367   if(fDebugLevel>=2) Printf("\n Making training lists! \n");
 368
 369   Int_t nPart[AliPID::kSPECIES][AliTRDCalPID::kNMom];
 370   memset(nPart, 0, AliPID::kSPECIES*AliTRDCalPID::kNMom*sizeof(Int_t));
 371
 372   // set needed branches
 373   fNN -> SetBranchAddress("fv0pid", &fv0pid);
 374   fNN -> SetBranchAddress("fMom", &fMom);
 375   fNN -> SetBranchAddress("fLayer", &fLayer);
 376
 377   // start first loop to check total number of each particle type
 378   for(Int_t iEv=0; iEv < fNN -> GetEntries(); iEv++){
 379     fNN -> GetEntry(iEv);
 380
 381     // use only events with goes through 6 layers TRD
 382     if(!fLayer == 0)
 383       continue;
 384
 385     // set the 11 momentum bins
 386     Int_t iMomBin = -1;
 387     if(fMom < .7) iMomBin = 0;
 388     else if(fMom < .9) iMomBin = 1;
 389     else if(fMom < 1.25) iMomBin = 2;
 390     else if(fMom < 1.75) iMomBin = 3;
 391     else if(fMom < 2.5) iMomBin = 4;
 392     else if(fMom < 3.5) iMomBin = 5;
 393     else if(fMom < 4.5) iMomBin = 6;
 394     else if(fMom < 5.5) iMomBin = 7;
 395     else if(fMom < 7.0) iMomBin = 8;
 396     else if(fMom < 9.0) iMomBin = 9;
 397     else  iMomBin = 10;
 398
 399     // check PID information and count particle types per momentum interval
 400     if(fv0pid[0] == 1)
 401       nPart[0][iMomBin]++;
 402     else if(fv0pid[1] == 1)
 403       nPart[1][iMomBin]++;
 404     else if(fv0pid[2] == 1)
 405       nPart[2][iMomBin]++;
 406     else if(fv0pid[3] == 1)
 407       nPart[3][iMomBin]++;
 408     else if(fv0pid[4] == 1)
 409       nPart[4][iMomBin]++;
 410   }
 411
 412   if(fDebugLevel>=2){
 413     Printf("Particle multiplicities:");
 414     for(Int_t iMomBin = 0; iMomBin <AliTRDCalPID::kNMom; iMomBin++)
 415       Printf("Momentum[%d]  Elecs[%d] Muons[%d] Pions[%d] Kaons[%d] Protons[%d]", iMomBin, nPart[0][iMomBin], nPart[1][iMomBin], nPart[2][iMomBin], nPart[3][iMomBin], nPart[4][iMomBin]);
 416     Printf("\n");
 417   }
 418
 419   // implement counter of training and test sample size
 420   Int_t iTrain[AliTRDCalPID::kNMom], iTest[AliTRDCalPID::kNMom];
 421   memset(iTrain, 0, AliTRDCalPID::kNMom*sizeof(Int_t));
 422   memset(iTest, 0, AliTRDCalPID::kNMom*sizeof(Int_t));
 423
 424   // set training sample size per momentum interval to 2/3
 425   // of smallest particle counter and test sample to 1/3
 426   for(Int_t iMomBin = 0; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
 427     iTrain[iMomBin] = nPart[0][iMomBin];
 428     for(Int_t iPart = 1; iPart < AliPID::kSPECIES; iPart++){
 429       if(iTrain[iMomBin] > nPart[iPart][iMomBin])
 430         iTrain[iMomBin] = nPart[iPart][iMomBin];
 431     }
 432     iTrain[iMomBin] = Int_t(iTrain[iMomBin] * .66);
 433     iTest[iMomBin] = Int_t( iTrain[iMomBin] * .5);
 434     if(fDebugLevel>=2) Printf("Momentum[%d]  Train[%d] Test[%d]", iMomBin, iTrain[iMomBin], iTest[iMomBin]);
 435   }
 436   if(fDebugLevel>=2) Printf("\n");
 437
 438
 439   // reset couters
 440   memset(nPart, 0, AliPID::kSPECIES*AliTRDCalPID::kNMom*sizeof(Int_t));
 441
 442   // start second loop to set the event lists
 443   for(Int_t iEv = 0; iEv < fNN -> GetEntries(); iEv++){
 444     fNN -> GetEntry(iEv);
 445
 446     // use only events with goes through 6 layers TRD
 447     if(!fLayer == 0)
 448       continue;
 449
 450     // set the 11 momentum bins
 451     Int_t iMomBin = -1;
 452     if(fMom < .7) iMomBin = 0;
 453     else if(fMom < .9) iMomBin = 1;
 454     else if(fMom < 1.25) iMomBin = 2;
 455     else if(fMom < 1.75) iMomBin = 3;
 456     else if(fMom < 2.5) iMomBin = 4;
 457     else if(fMom < 3.5) iMomBin = 5;
 458     else if(fMom < 4.5) iMomBin = 6;
 459     else if(fMom < 5.5) iMomBin = 7;
 460     else if(fMom < 7.0) iMomBin = 8;
 461     else if(fMom < 9.0) iMomBin = 9;
 462     else  iMomBin = 10;
 463
 464     // set electrons
 465     if(fv0pid[0] == 1){
 466       if(nPart[0][iMomBin] < iTrain[iMomBin]){
 467         for(Int_t ily = 0; ily < AliTRDCalPID::kNPlane; ily++)
 468           fTrain[iMomBin][ily] -> Enter(iEv + ily);
 469         nPart[0][iMomBin]++;
 470       }
 471       else if(nPart[0][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){
 472         for(Int_t ily = 0; ily < AliTRDCalPID::kNPlane; ily++)
 473           fTest[iMomBin][ily] -> Enter(iEv + ily);
 474         nPart[0][iMomBin]++;
 475       }
 476       else
 477         continue;
 478     }
 479     // set muons
 480     else if(fv0pid[1] == 1){
 481       if(nPart[1][iMomBin] < iTrain[iMomBin]){
 482         for(Int_t ily = 0; ily < AliTRDCalPID::kNPlane; ily++)
 483           fTrain[iMomBin][ily] -> Enter(iEv + ily);
 484         nPart[1][iMomBin]++;
 485       }
 486       else if(nPart[1][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){
 487         for(Int_t ily = 0; ily < AliTRDCalPID::kNPlane; ily++)
 488           fTest[iMomBin][ily] -> Enter(iEv + ily);
 489         nPart[1][iMomBin]++;
 490       }
 491       else
 492         continue;
 493     }
 494     // set pions
 495     else if(fv0pid[2] == 1){
 496       if(nPart[2][iMomBin] < iTrain[iMomBin]){
 497         for(Int_t ily = 0; ily < AliTRDCalPID::kNPlane; ily++)
 498           fTrain[iMomBin][ily] -> Enter(iEv + ily);
 499         nPart[2][iMomBin]++;
 500       }
 501       else if(nPart[2][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){
 502         for(Int_t ily = 0; ily < AliTRDCalPID::kNPlane; ily++)
 503           fTest[iMomBin][ily] -> Enter(iEv + ily);
 504         nPart[2][iMomBin]++;
 505       }
 506       else
 507         continue;
 508     }
 509     // set kaons
 510     else if(fv0pid[3] == 1){
 511       if(nPart[3][iMomBin] < iTrain[iMomBin]){
 512         for(Int_t ily = 0; ily < AliTRDCalPID::kNPlane; ily++)
 513           fTrain[iMomBin][ily] -> Enter(iEv + ily);
 514         nPart[3][iMomBin]++;
 515       }
 516       else if(nPart[3][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){
 517         for(Int_t ily = 0; ily < AliTRDCalPID::kNPlane; ily++)
 518           fTest[iMomBin][ily] -> Enter(iEv + ily);
 519         nPart[3][iMomBin]++;
 520       }
 521       else
 522         continue;
 523     }
 524     // set protons
 525     else if(fv0pid[4] == 1){
 526       if(nPart[4][iMomBin] < iTrain[iMomBin]){
 527         for(Int_t ily = 0; ily < AliTRDCalPID::kNPlane; ily++)
 528           fTrain[iMomBin][ily] -> Enter(iEv + ily);
 529         nPart[4][iMomBin]++;
 530       }
 531       else if(nPart[4][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){
 532         for(Int_t ily = 0; ily < AliTRDCalPID::kNPlane; ily++)
 533           fTest[iMomBin][ily] -> Enter(iEv + ily);
 534         nPart[4][iMomBin]++;
 535       }
 536       else
 537         continue;
 538     }
 539   }
 540
 541   if(fDebugLevel>=2){
 542     Printf("Particle multiplicities in both lists:");
 543     for(Int_t iMomBin = 0; iMomBin <AliTRDCalPID::kNMom; iMomBin++)
 544       Printf("Momentum[%d]  Elecs[%d] Muons[%d] Pions[%d] Kaons[%d] Protons[%d]", iMomBin, nPart[0][iMomBin], nPart[1][iMomBin], nPart[2][iMomBin], nPart[3][iMomBin], nPart[4][iMomBin]);
 545     Printf("\n");
 546   }
 547 }
 548
 549
 550 //________________________________________________________________________
 551 void AliTRDpidRefMaker::TrainNetworks(Int_t mombin)
 552 {
 553   //
 554   // train the neural networks
 555   //
 556
 557   if(fDebugLevel>=2) Printf("Training momentum bin %d", mombin);
 558
 559   // set variable to monitor the training and to save the development of the networks
 560   Int_t nEpochs = fEpochs/kMoniTrain;
 561   if(fDebugLevel>=2) Printf("Training %d times %d epochs", kMoniTrain, nEpochs);
 562
 563   // make directories to save the networks
 564   gSystem->Exec(Form("rm -r ./Networks_%d/MomBin_%d",fDate, mombin));
 565   gSystem->Exec(Form("mkdir ./Networks_%d/MomBin_%d",fDate, mombin));
 566
 567   // variable to check if network can load weights from previous training
 568   Bool_t bFirstLoop[AliTRDCalPID::kNPlane];
 569   memset(bFirstLoop, kTRUE, AliTRDCalPID::kNPlane*sizeof(Bool_t));
 570
 571   // train networks over several loops and save them after each loop
 572   for(Int_t iLoop = 0; iLoop < kMoniTrain; iLoop++){
 573     // loop over chambers
 574     for(Int_t iChamb = 0; iChamb < AliTRDCalPID::kNPlane; iChamb++){
 575       // set the event lists
 576       fNN -> SetEventList(fTrain[mombin][iChamb]);
 577       fNN -> SetEventList(fTest[mombin][iChamb]);
 578
 579       if(fDebugLevel>=2) Printf("Trainingloop[%d] Chamber[%d]", iLoop, iChamb);
 580
 581       // check if network is already implemented
 582       if(bFirstLoop[iChamb] == kTRUE){
 583         fNet[iChamb] = new TMultiLayerPerceptron("fdEdx[0],fdEdx[1],fdEdx[2],fdEdx[3],fdEdx[4],fdEdx[5],fdEdx[6],fdEdx[7]:15:7:fv0pid[0],fv0pid[1],fv0pid[2],fv0pid[3],fv0pid[4]!",fNN,fTrain[mombin][iChamb],fTest[mombin][iChamb]);
 584         fNet[iChamb] -> SetLearningMethod(TMultiLayerPerceptron::kStochastic);       // set learning method
 585         fNet[iChamb] -> TMultiLayerPerceptron::SetEta(0.001);                        // set learning speed
 586         if(fDebugLevel>=2) fNet[iChamb] -> Train(nEpochs,"text update=10, graph");
 587         else fNet[iChamb] -> Train(nEpochs,"");
 588         bFirstLoop[iChamb] = kFALSE;
 589       }
 590       else{
 591         if(fDebugLevel>=2) fNet[iChamb] -> Train(nEpochs,"text update=10, graph");
 592         else fNet[iChamb] -> Train(nEpochs,"+");
 593       }
 594
 595       // save weights for monitoring of the training
 596       fNet[iChamb] -> DumpWeights(Form("./Networks_%d/MomBin_%d/Net%d_%d",fDate, mombin, iChamb, iLoop));
 597     } // end chamber loop
 598   }   // end training loop
 599 }
 600
 601
 602 //________________________________________________________________________
 603 void AliTRDpidRefMaker::BuildLQRefs(Int_t mombin)
 604 {
 605   //
 606   // build the 2-dim LQ reference histograms
 607   //
 608
 609   if(fDebugLevel>=2) Printf("Building LQRefs for momentum bin %d", mombin);
 610 }
 611
 612
 613 //________________________________________________________________________
 614 void AliTRDpidRefMaker::MonitorTraining(Int_t /*mombin*/)
 615 {
 616   //
 617   // train the neural networks
 618   //
 619
 620 }