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