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 "TGraphErrors.h"
   7 #include "TTree.h"
   8 #include "TTreeStream.h"
   9 #include "TEventList.h"
  10 #include "TMultiLayerPerceptron.h"
  11
  12 #include "AliPID.h"
  13 #include "AliESDEvent.h"
  14 #include "AliESDInputHandler.h"
  15 #include "AliTrackReference.h"
  16
  17 #include "AliAnalysisTask.h"
  18
  19 #include "AliTRDtrackV1.h"
  20 #include "AliTRDReconstructor.h"
  21 #include "../Cal/AliTRDCalPID.h"
  22 #include "../Cal/AliTRDCalPIDNN.h"
  23
  24 #include "AliTRDpidUtil.h"
  25
  26 #include "AliTRDpidRefMaker.h"
  27 #include "info/AliTRDtrackInfo.h"
  28 #include "info/AliTRDv0Info.h"
  29
  30 // builds the reference tree for the training of neural networks
  31
  32
  33 ClassImp(AliTRDpidRefMaker)
  34
  35 //________________________________________________________________________
  36 AliTRDpidRefMaker::AliTRDpidRefMaker()
  37   :AliTRDrecoTask("PidRefMaker", "PID(NN) Reference Maker")
  38   ,fReconstructor(0x0)
  39   ,fV0s(0x0)
  40   ,fNN(0x0)
  41   ,fLQ(0x0)
  42   ,fLayer(0xff)
  43   ,fTrainMomBin(kAll)
  44   ,fEpochs(1000)
  45   ,fMinTrain(100)
  46   ,fDate(0)
  47   ,fMom(0.)
  48   ,fDoTraining(0)
  49   ,fContinueTraining(0)
  50   ,fTrainPath(0x0)
  51 {
  52   //
  53   // Default constructor
  54   //
  55
  56   fReconstructor = new AliTRDReconstructor();
  57   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
  58   memset(fv0pid, 0, AliPID::kSPECIES*sizeof(Float_t));
  59   memset(fdEdx, 0, 10*sizeof(Float_t));
  60
  61   const Int_t nnSize = AliTRDCalPID::kNMom * AliTRDgeometry::kNlayer;
  62   memset(fTrain, 0, nnSize*sizeof(TEventList*));
  63   memset(fTest, 0, nnSize*sizeof(TEventList*));
  64   memset(fNet, 0, AliTRDgeometry::kNlayer*sizeof(TMultiLayerPerceptron*));
  65
  66   TDatime datime;
  67   fDate = datime.GetDate();
  68
  69   DefineInput(1, TObjArray::Class());
  70   DefineOutput(1, TTree::Class());
  71   DefineOutput(2, TTree::Class());
  72 }
  73
  74
  75 //________________________________________________________________________
  76 AliTRDpidRefMaker::~AliTRDpidRefMaker()
  77 {
  78   if(fReconstructor) delete fReconstructor;
  79   //if(fNN) delete fNN;
  80   //if(fLQ) delete fLQ;
  81 }
  82
  83
  84 //________________________________________________________________________
  85 void AliTRDpidRefMaker::ConnectInputData(Option_t *opt)
  86 {
  87   AliTRDrecoTask::ConnectInputData(opt);
  88   fV0s = dynamic_cast<TObjArray*>(GetInputData(1));
  89 }
  90
  91 //________________________________________________________________________
  92 void AliTRDpidRefMaker::CreateOutputObjects()
  93 {
  94   // Create histograms
  95   // Called once
  96
  97   OpenFile(0, "RECREATE");
  98   fContainer = new TObjArray();
  99   fContainer->AddAt(new TH1F("hPDG","hPDG",AliPID::kSPECIES,-0.5,5.5),0);
 100
 101   TGraphErrors *gEffisTrain = new TGraphErrors(kMoniTrain);
 102   gEffisTrain -> SetLineColor(4);
 103   gEffisTrain -> SetMarkerColor(4);
 104   gEffisTrain -> SetMarkerStyle(29);
 105   gEffisTrain -> SetMarkerSize(1);
 106
 107   TGraphErrors *gEffisTest = new TGraphErrors(kMoniTrain);
 108   gEffisTest -> SetLineColor(2);
 109   gEffisTest -> SetMarkerColor(2);
 110   gEffisTest -> SetMarkerStyle(29);
 111   gEffisTest -> SetMarkerSize(1);
 112
 113   fContainer -> AddAt(gEffisTrain,kGraphTrain);
 114   fContainer -> AddAt(gEffisTest,kGraphTest);
 115
 116   // open reference TTree for NN
 117   OpenFile(1, "RECREATE");
 118   fNN = new TTree("NN", "Reference data for NN");
 119   fNN->Branch("fLayer", &fLayer, "fLayer/I");
 120   fNN->Branch("fMom", &fMom, "fMom/F");
 121   fNN->Branch("fv0pid", fv0pid, Form("fv0pid[%d]/F", AliPID::kSPECIES));
 122   fNN->Branch("fdEdx", fdEdx, Form("fdEdx[%d]/F", AliTRDpidUtil::kNNslices));
 123
 124   // open reference TTree for LQ
 125   OpenFile(2, "RECREATE");
 126   fLQ = new TTree("LQ", "Reference data for LQ");
 127   fLQ->Branch("fLayer", &fLayer, "fLayer/I");
 128   fLQ->Branch("fMom", &fMom, "fMom/F");
 129   fLQ->Branch("fv0pid", fv0pid, Form("fv0pid[%d]/F", AliPID::kSPECIES));
 130   fLQ->Branch("fdEdx", fdEdx, Form("fdEdx[%d]/F", AliTRDpidUtil::kLQslices));
 131 }
 132
 133
 134 //________________________________________________________________________
 135 void AliTRDpidRefMaker::Exec(Option_t *)
 136 {
 137   // Main loop
 138   // Called for each event
 139
 140   Int_t labelsacc[10000];
 141   memset(labelsacc, 0, sizeof(Int_t) * 10000);
 142
 143   Float_t mom;
 144   ULong_t status;
 145   Int_t nTRD = 0;
 146
 147   AliTRDtrackInfo     *track = 0x0;
 148   AliTRDv0Info           *v0 = 0x0;
 149   AliTRDtrackV1    *TRDtrack = 0x0;
 150   AliTrackReference     *ref = 0x0;
 151   AliExternalTrackParam *esd = 0x0;
 152   AliTRDseedV1 *TRDtracklet = 0x0;
 153
 154   for(Int_t iv0=0; iv0<fV0s->GetEntriesFast(); iv0++){
 155     v0 = dynamic_cast<AliTRDv0Info*>(fV0s->At(iv0));
 156     v0->Print();
 157   }
 158
 159   for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
 160
 161     // reset the pid information
 162     for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++)
 163       fv0pid[iPart] = 0.;
 164
 165     track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
 166     if(!track->HasESDtrack()) continue;
 167     status = track->GetStatus();
 168     if(!(status&AliESDtrack::kTPCout)) continue;
 169
 170     if(!(TRDtrack = track->GetTrack())) continue;
 171     //&&(track->GetNumberOfClustersRefit()
 172
 173     // use only tracks that hit 6 chambers
 174     if(!(TRDtrack->GetNumberOfTracklets() == AliTRDgeometry::kNlayer)) continue;
 175
 176     ref = track->GetTrackRef(0);
 177     esd = track->GetESDinfo()->GetOuterParam();
 178     mom = ref ? ref->P(): esd->P();
 179     fMom = mom;
 180
 181
 182     labelsacc[nTRD] = track->GetLabel();
 183     nTRD++;
 184
 185     // if no monte carlo data available -> use V0 information
 186     if(!HasMCdata()){
 187       GetV0info(TRDtrack,fv0pid);
 188     }
 189     // else use the MC info
 190     else{
 191       switch(track -> GetPDG()){
 192       case kElectron:
 193       case kPositron:
 194         fv0pid[AliPID::kElectron] = 1.;
 195         break;
 196       case kMuonPlus:
 197       case kMuonMinus:
 198         fv0pid[AliPID::kMuon] = 1.;
 199         break;
 200       case kPiPlus:
 201       case kPiMinus:
 202         fv0pid[AliPID::kPion] = 1.;
 203         break;
 204       case kKPlus:
 205       case kKMinus:
 206         fv0pid[AliPID::kKaon] = 1.;
 207         break;
 208       case kProton:
 209       case kProtonBar:
 210         fv0pid[AliPID::kProton] = 1.;
 211         break;
 212       }
 213     }
 214
 215     // set reconstructor
 216     Float_t *dedx;
 217     TRDtrack -> SetReconstructor(fReconstructor);
 218
 219     // fill the dE/dx information for NN
 220     fReconstructor -> SetOption("nn");
 221     for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){
 222       if(!(TRDtracklet = TRDtrack -> GetTracklet(ily))) continue;
 223       TRDtracklet->CookdEdx(AliTRDpidUtil::kNNslices);
 224       dedx = TRDtracklet->GetdEdx();
 225       for(Int_t iSlice = 0; iSlice < AliTRDpidUtil::kNNslices; iSlice++)
 226         dedx[iSlice] = dedx[iSlice]/AliTRDCalPIDNN::kMLPscale;
 227       memcpy(fdEdx, dedx, AliTRDpidUtil::kNNslices*sizeof(Float_t));
 228       if(fDebugLevel>=2) Printf("LayerNN : %d", ily);
 229       fLayer = ily;
 230       fNN->Fill();
 231     }
 232
 233
 234     // fill the dE/dx information for LQ
 235     fReconstructor -> SetOption("!nn");
 236     for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){
 237       if(!(TRDtracklet = TRDtrack -> GetTracklet(ily))) continue;
 238       TRDtracklet->CookdEdx(AliTRDpidUtil::kLQslices);
 239       dedx = TRDtracklet->GetdEdx();
 240       memcpy(fdEdx, dedx, AliTRDpidUtil::kLQslices*sizeof(Float_t));
 241       if(fDebugLevel>=2) Printf("LayerLQ : %d", ily);
 242       fLayer = ily;
 243       fLQ->Fill();
 244     }
 245
 246
 247     for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
 248       if(fDebugLevel>=4) Printf("PDG is %d %f", iPart, fv0pid[iPart]);
 249     }
 250   }
 251
 252   PostData(0, fContainer);
 253   PostData(1, fNN);
 254   PostData(2, fLQ);
 255 }
 256
 257
 258 //________________________________________________________________________
 259 Bool_t AliTRDpidRefMaker::PostProcess()
 260 {
 261   // Draw result to the screen
 262   // Called once at the end of the query
 263
 264   // build the training andthe test list for the neural networks
 265   MakeTrainingLists();
 266   if(!fDoTraining) return kTRUE;
 267
 268   // train the neural networks and build the refrence histos for 2-dim LQ
 269   gSystem->Exec(Form("mkdir ./Networks_%d/",fDate));
 270   if(fDebugLevel>=2) Printf("TrainMomBin [%d] [%d]", fTrainMomBin, kAll);
 271
 272   // train single network for a single momentum (recommended)
 273   if(!(fTrainMomBin == kAll)){
 274     if(fTrain[fTrainMomBin][0] -> GetN() < fMinTrain){
 275       if(fDebugLevel>=2) Printf("Warning in AliTRDpidRefMaker::PostProcess : Not enough events for training available! Please check Data sample!");
 276       return kFALSE;
 277     }
 278     TrainNetworks(fTrainMomBin);
 279     BuildLQRefs(fTrainMomBin);
 280     MonitorTraining(fTrainMomBin);
 281   }
 282   // train all momenta
 283   else{
 284     for(Int_t iMomBin = 0; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
 285       if(fTrain[iMomBin][0] -> GetN() < fMinTrain){
 286         if(fDebugLevel>=2) Printf("Warning in AliTRDpidRefMaker::PostProcess : Not enough events for training available for momentum bin [%d]! Please check Data sample!", iMomBin);
 287         continue;
 288       }
 289       TrainNetworks(iMomBin);
 290       BuildLQRefs(fTrainMomBin);
 291       MonitorTraining(iMomBin);
 292     }
 293   }
 294
 295   return kTRUE; // testing protection
 296 }
 297
 298
 299 //________________________________________________________________________
 300 void AliTRDpidRefMaker::Terminate(Option_t *)
 301 {
 302   // Draw result to the screen
 303   // Called once at the end of the query
 304
 305   fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
 306   if (!fContainer) {
 307     Printf("ERROR: list not available");
 308     return;
 309   }
 310 }
 311
 312
 313 //________________________________________________________________________
 314 void AliTRDpidRefMaker::GetV0info(AliTRDtrackV1 *TRDtrack, Float_t *v0pid)
 315 {
 316   // !!!! PREMILMINARY FUNCTION !!!!
 317   //
 318   // this is the place for the V0 procedure
 319   // as long as there is no one implemented,
 320   // just the probabilities
 321   // of the TRDtrack are used!
 322
 323   TRDtrack -> SetReconstructor(fReconstructor);
 324   fReconstructor -> SetOption("nn");
 325   TRDtrack -> CookPID();
 326   for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
 327     v0pid[iPart] = TRDtrack -> GetPID(iPart);
 328     if(fDebugLevel>=4) Printf("PDG is (in V0info) %d %f", iPart, v0pid[iPart]);
 329   }
 330 }
 331
 332
 333 //________________________________________________________________________
 334 void AliTRDpidRefMaker::MakeTrainingLists()
 335 {
 336   //
 337   // build the training lists for the neural networks
 338   //
 339
 340   if (!fNN) {
 341     LoadFiles("TRD.TaskPidRefMakerNN.root", "TRD.TaskPidRefMakerLQ.root");
 342   }
 343
 344   if (!fNN) {
 345     Printf("ERROR tree for training list not available");
 346     return;
 347   }
 348
 349   if(fDebugLevel>=2) Printf("\n Making training lists! \n");
 350
 351   Int_t nPart[AliPID::kSPECIES][AliTRDCalPID::kNMom];
 352   memset(nPart, 0, AliPID::kSPECIES*AliTRDCalPID::kNMom*sizeof(Int_t));
 353
 354   // set needed branches
 355   fNN -> SetBranchAddress("fv0pid", &fv0pid);
 356   fNN -> SetBranchAddress("fMom", &fMom);
 357   fNN -> SetBranchAddress("fLayer", &fLayer);
 358
 359   AliTRDpidUtil *util = new AliTRDpidUtil();
 360
 361   // start first loop to check total number of each particle type
 362   for(Int_t iEv=0; iEv < fNN -> GetEntries(); iEv++){
 363     fNN -> GetEntry(iEv);
 364
 365     // use only events with goes through 6 layers TRD
 366     if(!fLayer == 0)
 367       continue;
 368
 369     // set the 11 momentum bins
 370     Int_t iMomBin = -1;
 371     iMomBin = util -> GetMomentumBin(fMom);
 372
 373     // check PID information and count particle types per momentum interval
 374     if(fv0pid[AliPID::kElectron] == 1)
 375       nPart[AliPID::kElectron][iMomBin]++;
 376     else if(fv0pid[AliPID::kMuon] == 1)
 377       nPart[AliPID::kMuon][iMomBin]++;
 378     else if(fv0pid[AliPID::kPion] == 1)
 379       nPart[AliPID::kPion][iMomBin]++;
 380     else if(fv0pid[AliPID::kKaon] == 1)
 381       nPart[AliPID::kKaon][iMomBin]++;
 382     else if(fv0pid[AliPID::kProton] == 1)
 383       nPart[AliPID::kProton][iMomBin]++;
 384   }
 385
 386   if(fDebugLevel>=2){
 387     Printf("Particle multiplicities:");
 388     for(Int_t iMomBin = 0; iMomBin <AliTRDCalPID::kNMom; iMomBin++)
 389       Printf("Momentum[%d]  Elecs[%d] Muons[%d] Pions[%d] Kaons[%d] Protons[%d]", iMomBin, nPart[AliPID::kElectron][iMomBin], nPart[AliPID::kMuon][iMomBin], nPart[AliPID::kPion][iMomBin], nPart[AliPID::kKaon][iMomBin], nPart[AliPID::kProton][iMomBin]);
 390     Printf("\n");
 391   }
 392
 393   // implement counter of training and test sample size
 394   Int_t iTrain[AliTRDCalPID::kNMom], iTest[AliTRDCalPID::kNMom];
 395   memset(iTrain, 0, AliTRDCalPID::kNMom*sizeof(Int_t));
 396   memset(iTest, 0, AliTRDCalPID::kNMom*sizeof(Int_t));
 397
 398   // set training sample size per momentum interval to 2/3
 399   // of smallest particle counter and test sample to 1/3
 400   for(Int_t iMomBin = 0; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
 401     iTrain[iMomBin] = nPart[0][iMomBin];
 402     for(Int_t iPart = 1; iPart < AliPID::kSPECIES; iPart++){
 403       if(iTrain[iMomBin] > nPart[iPart][iMomBin])
 404         iTrain[iMomBin] = nPart[iPart][iMomBin];
 405     }
 406     iTrain[iMomBin] = Int_t(iTrain[iMomBin] * .66);
 407     iTest[iMomBin] = Int_t( iTrain[iMomBin] * .5);
 408     if(fDebugLevel>=2) Printf("Momentum[%d]  Train[%d] Test[%d]", iMomBin, iTrain[iMomBin], iTest[iMomBin]);
 409   }
 410   if(fDebugLevel>=2) Printf("\n");
 411
 412
 413   // reset couters
 414   memset(nPart, 0, AliPID::kSPECIES*AliTRDCalPID::kNMom*sizeof(Int_t));
 415
 416   // start second loop to set the event lists
 417   for(Int_t iEv = 0; iEv < fNN -> GetEntries(); iEv++){
 418     fNN -> GetEntry(iEv);
 419
 420     // use only events with goes through 6 layers TRD
 421     if(!fLayer == 0)
 422       continue;
 423
 424     // set the 11 momentum bins
 425     Int_t iMomBin = -1;
 426     iMomBin = util -> GetMomentumBin(fMom);
 427
 428     // set electrons
 429     if(fv0pid[AliPID::kElectron] == 1){
 430       if(nPart[AliPID::kElectron][iMomBin] < iTrain[iMomBin]){
 431         for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
 432           fTrain[iMomBin][ily] -> Enter(iEv + ily);
 433         nPart[AliPID::kElectron][iMomBin]++;
 434       }
 435       else if(nPart[AliPID::kElectron][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){
 436         for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
 437           fTest[iMomBin][ily] -> Enter(iEv + ily);
 438         nPart[AliPID::kElectron][iMomBin]++;
 439       }
 440       else
 441         continue;
 442     }
 443     // set muons
 444     else if(fv0pid[AliPID::kMuon] == 1){
 445       if(nPart[AliPID::kMuon][iMomBin] < iTrain[iMomBin]){
 446         for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
 447           fTrain[iMomBin][ily] -> Enter(iEv + ily);
 448         nPart[AliPID::kMuon][iMomBin]++;
 449       }
 450       else if(nPart[AliPID::kMuon][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){
 451         for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
 452           fTest[iMomBin][ily] -> Enter(iEv + ily);
 453         nPart[AliPID::kMuon][iMomBin]++;
 454       }
 455       else
 456         continue;
 457     }
 458     // set pions
 459     else if(fv0pid[AliPID::kPion] == 1){
 460       if(nPart[AliPID::kPion][iMomBin] < iTrain[iMomBin]){
 461         for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
 462           fTrain[iMomBin][ily] -> Enter(iEv + ily);
 463         nPart[AliPID::kPion][iMomBin]++;
 464       }
 465       else if(nPart[AliPID::kPion][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){
 466         for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
 467           fTest[iMomBin][ily] -> Enter(iEv + ily);
 468         nPart[AliPID::kPion][iMomBin]++;
 469       }
 470       else
 471         continue;
 472     }
 473     // set kaons
 474     else if(fv0pid[AliPID::kKaon] == 1){
 475       if(nPart[AliPID::kKaon][iMomBin] < iTrain[iMomBin]){
 476         for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
 477           fTrain[iMomBin][ily] -> Enter(iEv + ily);
 478         nPart[AliPID::kKaon][iMomBin]++;
 479       }
 480       else if(nPart[AliPID::kKaon][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){
 481         for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
 482           fTest[iMomBin][ily] -> Enter(iEv + ily);
 483         nPart[AliPID::kKaon][iMomBin]++;
 484       }
 485       else
 486         continue;
 487     }
 488     // set protons
 489     else if(fv0pid[AliPID::kProton] == 1){
 490       if(nPart[AliPID::kProton][iMomBin] < iTrain[iMomBin]){
 491         for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
 492           fTrain[iMomBin][ily] -> Enter(iEv + ily);
 493         nPart[AliPID::kProton][iMomBin]++;
 494       }
 495       else if(nPart[AliPID::kProton][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){
 496         for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
 497           fTest[iMomBin][ily] -> Enter(iEv + ily);
 498         nPart[AliPID::kProton][iMomBin]++;
 499       }
 500       else
 501         continue;
 502     }
 503   }
 504
 505   if(fDebugLevel>=2){
 506     Printf("Particle multiplicities in both lists:");
 507     for(Int_t iMomBin = 0; iMomBin <AliTRDCalPID::kNMom; iMomBin++)
 508       Printf("Momentum[%d]  Elecs[%d] Muons[%d] Pions[%d] Kaons[%d] Protons[%d]", iMomBin, nPart[AliPID::kElectron][iMomBin], nPart[AliPID::kMuon][iMomBin], nPart[AliPID::kPion][iMomBin], nPart[AliPID::kKaon][iMomBin], nPart[AliPID::kProton][iMomBin]);
 509     Printf("\n");
 510   }
 511
 512   util -> Delete();
 513 }
 514
 515
 516 //________________________________________________________________________
 517 void AliTRDpidRefMaker::TrainNetworks(Int_t mombin)
 518 {
 519   //
 520   // train the neural networks
 521   //
 522
 523
 524   if (!fNN) {
 525     LoadFiles("TRD.TaskPidRefMakerNN.root", "TRD.TaskPidRefMakerLQ.root");
 526   }
 527
 528   if (!fNN) {
 529     Printf("ERROR tree for training list not available");
 530     return;
 531   }
 532
 533   TDatime datime;
 534   fDate = datime.GetDate();
 535
 536   if(fDebugLevel>=2) Printf("Training momentum bin %d", mombin);
 537
 538   // set variable to monitor the training and to save the development of the networks
 539   Int_t nEpochs = fEpochs/kMoniTrain;
 540   if(fDebugLevel>=2) Printf("Training %d times %d epochs", kMoniTrain, nEpochs);
 541
 542   // make directories to save the networks
 543   gSystem->Exec(Form("rm -r ./Networks_%d/MomBin_%d",fDate, mombin));
 544   gSystem->Exec(Form("mkdir ./Networks_%d/MomBin_%d",fDate, mombin));
 545
 546   // variable to check if network can load weights from previous training
 547   Bool_t bFirstLoop[AliTRDgeometry::kNlayer];
 548   memset(bFirstLoop, kTRUE, AliTRDgeometry::kNlayer*sizeof(Bool_t));
 549
 550   // train networks over several loops and save them after each loop
 551   for(Int_t iLoop = 0; iLoop < kMoniTrain; iLoop++){
 552     // loop over chambers
 553     for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 554       // set the event lists
 555       fNN -> SetEventList(fTrain[mombin][iChamb]);
 556       fNN -> SetEventList(fTest[mombin][iChamb]);
 557
 558       if(fDebugLevel>=2) Printf("Trainingloop[%d] Chamber[%d]", iLoop, iChamb);
 559
 560       // check if network is already implemented
 561       if(bFirstLoop[iChamb] == kTRUE){
 562         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]);
 563         fNet[iChamb] -> SetLearningMethod(TMultiLayerPerceptron::kStochastic);       // set learning method
 564         fNet[iChamb] -> TMultiLayerPerceptron::SetEta(0.001);                        // set learning speed
 565         if(!fContinueTraining){
 566           if(fDebugLevel>=2) fNet[iChamb] -> Train(nEpochs,"text update=10, graph");
 567           else fNet[iChamb] -> Train(nEpochs,"");
 568         }
 569         else{
 570           fNet[iChamb] -> LoadWeights(Form("./Networks_%d/MomBin_%d/Net%d_%d",fTrainPath, mombin, iChamb, kMoniTrain - 1));
 571           if(fDebugLevel>=2) fNet[iChamb] -> Train(nEpochs,"text update=10, graph+");
 572           else fNet[iChamb] -> Train(nEpochs,"+");
 573         }
 574         bFirstLoop[iChamb] = kFALSE;
 575       }
 576       else{
 577         if(fDebugLevel>=2) fNet[iChamb] -> Train(nEpochs,"text update=10, graph+");
 578         else fNet[iChamb] -> Train(nEpochs,"+");
 579       }
 580
 581       // save weights for monitoring of the training
 582       fNet[iChamb] -> DumpWeights(Form("./Networks_%d/MomBin_%d/Net%d_%d",fDate, mombin, iChamb, iLoop));
 583     } // end chamber loop
 584   }   // end training loop
 585 }
 586
 587
 588 //________________________________________________________________________
 589 void AliTRDpidRefMaker::BuildLQRefs(Int_t mombin)
 590 {
 591   //
 592   // build the 2-dim LQ reference histograms
 593   //
 594
 595   if(fDebugLevel>=2) Printf("Building LQRefs for momentum bin %d", mombin);
 596 }
 597
 598
 599 //________________________________________________________________________
 600 void AliTRDpidRefMaker::MonitorTraining(Int_t mombin)
 601 {
 602   //
 603   // train the neural networks
 604   //
 605
 606   if(!fContainer){
 607     LoadContainer("TRD.TaskPidRefMaker.root");
 608   }
 609   if(!fContainer){
 610     Printf("ERROR container not available");
 611     return;
 612   }
 613
 614   if (!fNN) {
 615     LoadFiles("TRD.TaskPidRefMakerNN.root", "TRD.TaskPidRefMakerLQ.root");
 616   }
 617   if (!fNN) {
 618     Printf("ERROR tree for training list not available");
 619     return;
 620   }
 621
 622   // init networks and set event list
 623   for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 624         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]);
 625         fNN -> SetEventList(fTrain[mombin][iChamb]);
 626         fNN -> SetEventList(fTest[mombin][iChamb]);
 627   }
 628
 629   // implement variables for likelihoods
 630   Float_t Like[AliPID::kSPECIES][AliTRDgeometry::kNlayer];
 631   memset(Like, 0, AliPID::kSPECIES*AliTRDgeometry::kNlayer*sizeof(Float_t));
 632   Float_t LikeAll[AliPID::kSPECIES], TotProb;
 633
 634   Double_t PionEffiTrain[kMoniTrain], PionEffiErrTrain[kMoniTrain];
 635   Double_t PionEffiTest[kMoniTrain], PionEffiErrTest[kMoniTrain];
 636   memset(PionEffiTrain, 0, kMoniTrain*sizeof(Double_t));
 637   memset(PionEffiErrTrain, 0, kMoniTrain*sizeof(Double_t));
 638   memset(PionEffiTest, 0, kMoniTrain*sizeof(Double_t));
 639   memset(PionEffiErrTest, 0, kMoniTrain*sizeof(Double_t));
 640
 641   // init histos
 642   const Float_t epsilon = 1/(2*(AliTRDpidUtil::kBins-1));     // get nice histos with bin center at 0 and 1
 643   TH1F *hElecs, *hPions;
 644   hElecs = new TH1F("hElecs","Likelihood for electrons", AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon);
 645   hPions = new TH1F("hPions","Likelihood for pions", AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon);
 646
 647   TGraphErrors *gEffisTrain=0x0, *gEffisTest=0x0;
 648   gEffisTrain = (TGraphErrors*)fContainer->At(kGraphTrain);
 649   gEffisTest = (TGraphErrors*)fContainer->At(kGraphTest);
 650
 651   AliTRDpidUtil *util = new AliTRDpidUtil();
 652
 653   // monitor training progress
 654   for(Int_t iLoop = 0; iLoop < kMoniTrain; iLoop++){
 655
 656     // load weights
 657     for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 658       fNet[iChamb] -> LoadWeights(Form("./Networks_%d/MomBin_%d/Net%d_%d",fDate, mombin, iChamb, iLoop));
 659     }
 660
 661     // event loop training list
 662     for(Int_t iEvent = 0; iEvent < fTrain[mombin][0] -> GetN(); iEvent++ ){
 663
 664       // reset particle probabilities
 665       for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
 666         LikeAll[iPart] = 1./AliPID::kSPECIES;
 667       }
 668       TotProb = 0.;
 669
 670       fNN -> GetEntry(fTrain[mombin][0] -> GetEntry(iEvent));
 671       // use event only if it is electron or pion
 672       if(!((fv0pid[AliPID::kElectron] == 1.0) || (fv0pid[AliPID::kPion] == 1.0))) continue;
 673
 674       // get the probabilities for each particle type in each chamber
 675       for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 676         for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
 677           Like[iPart][iChamb] = fNet[iChamb] -> Result(fTrain[mombin][iChamb] -> GetEntry(iEvent), iPart);
 678           LikeAll[iPart] *=  Like[iPart][iChamb];
 679         }
 680       }
 681
 682       // get total probability and normalize it
 683       for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
 684         TotProb += LikeAll[iPart];
 685       }
 686       for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
 687         LikeAll[iPart] /= TotProb;
 688       }
 689
 690       // fill likelihood distributions
 691       if(fv0pid[AliPID::kElectron] == 1)
 692         hElecs -> Fill(LikeAll[AliPID::kElectron]);
 693       if(fv0pid[AliPID::kPion] == 1)
 694         hPions -> Fill(LikeAll[AliPID::kElectron]);
 695     } // end event loop
 696
 697
 698     // calculate the pion efficiency and fill the graph
 699     util -> CalculatePionEffi(hElecs, hPions);
 700     PionEffiTrain[iLoop] = util -> GetPionEfficiency();
 701     PionEffiErrTrain[iLoop] = util -> GetError();
 702
 703     gEffisTrain -> SetPoint(iLoop, iLoop+1, PionEffiTrain[iLoop]);
 704     gEffisTrain -> SetPointError(iLoop, 0, PionEffiErrTrain[iLoop]);
 705     hElecs -> Reset();
 706     hPions -> Reset();
 707     if(fDebugLevel>=2) Printf("TrainingLoop[%d] PionEfficiency[%f +/- %f]", iLoop, PionEffiTrain[iLoop], PionEffiErrTrain[iLoop]);
 708     // end training loop
 709
 710
 711
 712     // event loop test list
 713     for(Int_t iEvent = 0; iEvent < fTest[mombin][0] -> GetN(); iEvent++ ){
 714
 715       // reset particle probabilities
 716       for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
 717         LikeAll[iPart] = 1./AliTRDgeometry::kNlayer;
 718       }
 719       TotProb = 0.;
 720
 721       fNN -> GetEntry(fTest[mombin][0] -> GetEntry(iEvent));
 722       // use event only if it is electron or pion
 723       if(!((fv0pid[AliPID::kElectron] == 1.0) || (fv0pid[AliPID::kPion] == 1.0))) continue;
 724
 725       // get the probabilities for each particle type in each chamber
 726       for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 727         for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
 728           Like[iPart][iChamb] = fNet[iChamb] -> Result(fTest[mombin][iChamb] -> GetEntry(iEvent), iPart);
 729           LikeAll[iPart] *=  Like[iPart][iChamb];
 730         }
 731       }
 732
 733       // get total probability and normalize it
 734       for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
 735         TotProb += LikeAll[iPart];
 736       }
 737       for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
 738         LikeAll[iPart] /= TotProb;
 739       }
 740
 741       // fill likelihood distributions
 742       if(fv0pid[AliPID::kElectron] == 1)
 743         hElecs -> Fill(LikeAll[AliPID::kElectron]);
 744       if(fv0pid[AliPID::kPion] == 1)
 745         hPions -> Fill(LikeAll[AliPID::kElectron]);
 746     } // end event loop
 747
 748     // calculate the pion efficiency and fill the graph
 749     util -> CalculatePionEffi(hElecs, hPions);
 750     PionEffiTest[iLoop] = util -> GetPionEfficiency();
 751     PionEffiErrTest[iLoop] = util -> GetError();
 752
 753     gEffisTest -> SetPoint(iLoop, iLoop+1, PionEffiTest[iLoop]);
 754     gEffisTest -> SetPointError(iLoop, 0, PionEffiErrTest[iLoop]);
 755     hElecs -> Reset();
 756     hPions -> Reset();
 757     if(fDebugLevel>=2) Printf("TestLoop[%d] PionEfficiency[%f +/- %f] \n", iLoop, PionEffiTest[iLoop], PionEffiErrTest[iLoop]);
 758
 759   } //   end training loop
 760
 761   util -> Delete();
 762
 763   gEffisTest -> Draw("PAL");
 764   gEffisTrain -> Draw("PL");
 765
 766 }
 767
 768
 769 //________________________________________________________________________
 770 void AliTRDpidRefMaker::LoadFiles(const Char_t *InFileNN, const Char_t *InFileLQ)
 771 {
 772   //
 773   // Loads the files and sets the event list
 774   // for neural network training and
 775   // building of the 2-dim reference histograms.
 776   // Useable for training outside of the makeResults.C macro
 777   //
 778
 779   TFile *fInFileNN;
 780   fInFileNN = new TFile(InFileNN, "READ");
 781   fNN = (TTree*)fInFileNN -> Get("NN");
 782
 783   TFile *fInFileLQ;
 784   fInFileLQ = new TFile(InFileLQ, "READ");
 785   fLQ = (TTree*)fInFileLQ -> Get("LQ");
 786
 787   for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
 788     for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){
 789       fTrain[iMom][ily] = new TEventList(Form("fTrainMom%d_%d", iMom, ily), Form("Training list for momentum intervall %d and plane %d", iMom, ily));
 790       fTest[iMom][ily] = new TEventList(Form("fTestMom%d_%d", iMom, ily), Form("Test list for momentum intervall %d and plane %d", iMom, ily));
 791     }
 792   }
 793 }
 794
 795
 796 //________________________________________________________________________
 797 void AliTRDpidRefMaker::LoadContainer(const Char_t *InFileCont)
 798 {
 799
 800   //
 801   // Loads the container if no container is there.
 802   // Useable for training outside of the makeResults.C macro
 803   //
 804
 805   TFile *fInFileCont;
 806   fInFileCont = new TFile(InFileCont, "READ");
 807   fContainer = (TObjArray*)fInFileCont -> Get("PidRefMaker");
 808
 809 }
 810
 811
 812 // //________________________________________________________________________
 813 // void AliTRDpidRefMaker::CreateGraphs()
 814 // {
 815 //   // Create histograms
 816 //   // Called once
 817
 818 //   OpenFile(0, "RECREATE");
 819 //   fContainer = new TObjArray();
 820 //   fContainer->AddAt(new TH1F("hPDG","hPDG",AliPID::kSPECIES,-0.5,5.5),0);
 821
 822 //   TGraphErrors *gEffisTrain = new TGraphErrors(kMoniTrain);
 823 //   gEffisTrain -> SetLineColor(4);
 824 //   gEffisTrain -> SetMarkerColor(4);
 825 //   gEffisTrain -> SetMarkerStyle(29);
 826 //   gEffisTrain -> SetMarkerSize(2);
 827
 828 //   TGraphErrors *gEffisTest = new TGraphErrors(kMoniTrain);
 829 //   gEffisTest -> SetLineColor(2);
 830 //   gEffisTest -> SetMarkerColor(2);
 831 //   gEffisTest -> SetMarkerSize(2);
 832
 833 //   fContainer -> AddAt(gEffisTrain,kGraphTrain);
 834 //   fContainer -> AddAt(gEffisTest,kGraphTest);
 835 // }
 836