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