[u/mrichter/AliRoot.git] / TRD / qaRec / AliTRDpidRefMakerNN.cxx

#include "TSystem.h"
#include "TDatime.h"
#include "TEventList.h"
#include "TGraphErrors.h"
#include "TTreeStream.h"
#include "TH1F.h"
#include "TMultiLayerPerceptron.h"

#include "AliLog.h"
#include "AliPID.h"
#include "AliESDEvent.h"

#include "../Cal/AliTRDCalPIDNN.h"

#include "AliTRDpidUtil.h"

#include "AliTRDpidRefMakerNN.h"
#include "info/AliTRDtrackInfo.h"

// builds the reference tree for the training of neural networks


ClassImp(AliTRDpidRefMakerNN)

//________________________________________________________________________
AliTRDpidRefMakerNN::AliTRDpidRefMakerNN() 
  :AliTRDpidRefMaker(UChar_t(AliTRDpidUtil::kNNslices), "PidRefMakerNN", "PID(NN) Reference Maker")
  ,fTrainMomBin(kAll)
  ,fEpochs(1000)
  ,fMinTrain(100)
  ,fDate(0)
  ,fDoTraining(0)
  ,fContinueTraining(0)
  ,fTrainPath(0x0)
{
  //
  // Default constructor
  //
  SetAbundance(.67, .33);
  SetScaledEdx(Float_t(AliTRDCalPIDNN::kMLPscale));
  TDatime datime;
  fDate = datime.GetDate();
}


//________________________________________________________________________
AliTRDpidRefMakerNN::~AliTRDpidRefMakerNN() 
{
}


//________________________________________________________________________
void AliTRDpidRefMakerNN::CreateOutputObjects()
{
  // Create histograms
  // Called once
  AliTRDpidRefMaker::CreateOutputObjects();
  TGraphErrors *gEffisTrain = new TGraphErrors(kMoniTrain);
  gEffisTrain -> SetLineColor(4);
  gEffisTrain -> SetMarkerColor(4);
  gEffisTrain -> SetMarkerStyle(29);
  gEffisTrain -> SetMarkerSize(1);

  TGraphErrors *gEffisTest = new TGraphErrors(kMoniTrain);
  gEffisTest -> SetLineColor(2);
  gEffisTest -> SetMarkerColor(2);
  gEffisTest -> SetMarkerStyle(29);
  gEffisTest -> SetMarkerSize(1);

  fContainer -> AddAt(gEffisTrain,kGraphTrain);
  fContainer -> AddAt(gEffisTest,kGraphTest);
}


//________________________________________________________________________
Bool_t AliTRDpidRefMakerNN::PostProcess()
{
  // Draw result to the screen
  // Called once at the end of the query

  // build the training andthe test list for the neural networks
  MakeTrainingList();        
  if(!fDoTraining) return kTRUE;

  // train the neural networks and build the refrence histos for 2-dim LQ
  gSystem->Exec(Form("mkdir ./Networks_%d/",fDate));
  AliDebug(3, Form("TrainMomBin [%d] [%d]", fTrainMomBin, kAll));

  // train single network for a single momentum (recommended)
  if(!(fTrainMomBin == kAll)){
    if(fTrain[fTrainMomBin][0] -> GetN() < fMinTrain){
      AliWarning("Not enough events for training available! Please check Data sample!");
      return kFALSE;
    }
    MakeRefs(fTrainMomBin);
    MonitorTraining(fTrainMomBin);
  }
  // train all momenta
  else{
    for(Int_t iMomBin = 0; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
      if(fTrain[iMomBin][0] -> GetN() < fMinTrain){
        AliWarning(Form("Not enough events for training available for momentum bin [%d]! Please check Data sample!", iMomBin));
        continue;
      }
      MakeRefs(fTrainMomBin);
      MonitorTraining(iMomBin);
    }
  }

  return kTRUE; // testing protection
}


//________________________________________________________________________
void AliTRDpidRefMakerNN::MakeRefs(Int_t mombin) 
{
  //
  // train the neural networks
  //
  
  
  if (!fData) LoadFile(Form("TRD.%s.root", GetName()));

  if (!fData) {
    AliError("Tree for training list not available");
    return;
  }

  TDatime datime;
  fDate = datime.GetDate();

  AliDebug(2, Form("Training momentum bin %d", mombin));

  // set variable to monitor the training and to save the development of the networks
  Int_t nEpochs = fEpochs/kMoniTrain;       
  AliDebug(2, Form("Training %d times %d epochs", kMoniTrain, nEpochs));

  // make directories to save the networks 
  gSystem->Exec(Form("rm -r ./Networks_%d/MomBin_%d",fDate, mombin));
  gSystem->Exec(Form("mkdir ./Networks_%d/MomBin_%d",fDate, mombin));

  // variable to check if network can load weights from previous training
  Bool_t bFirstLoop[AliTRDgeometry::kNlayer];
  memset(bFirstLoop, kTRUE, AliTRDgeometry::kNlayer*sizeof(Bool_t));
 
  // train networks over several loops and save them after each loop
  for(Int_t iLoop = 0; iLoop < kMoniTrain; iLoop++){
    // loop over chambers
    for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
      // set the event lists
      fData -> SetEventList(fTrain[mombin][iChamb]);
      fData -> SetEventList(fTest[mombin][iChamb]);
      
      AliDebug(2, Form("Trainingloop[%d] Chamber[%d]", iLoop, iChamb));
      
      // check if network is already implemented
      if(bFirstLoop[iChamb] == kTRUE){
	fNet[iChamb] = new TMultiLayerPerceptron("fdEdx[0],fdEdx[1],fdEdx[2],fdEdx[3],fdEdx[4],fdEdx[5],fdEdx[6],fdEdx[7]:15:7:fPID[0],fPID[1],fPID[2],fPID[3],fPID[4]!",fData,fTrain[mombin][iChamb],fTest[mombin][iChamb]);
	fNet[iChamb] -> SetLearningMethod(TMultiLayerPerceptron::kStochastic);       // set learning method
	fNet[iChamb] -> TMultiLayerPerceptron::SetEta(0.001);                        // set learning speed
	if(!fContinueTraining){
	  if(fDebugLevel>=2) fNet[iChamb] -> Train(nEpochs,"text update=10, graph");
	  else fNet[iChamb] -> Train(nEpochs,"");
	}
	else{
	  fNet[iChamb] -> LoadWeights(Form("./Networks_%d/MomBin_%d/Net%d_%d",fTrainPath, mombin, iChamb, kMoniTrain - 1));
	  if(fDebugLevel>=2) fNet[iChamb] -> Train(nEpochs,"text update=10, graph+");      
	  else fNet[iChamb] -> Train(nEpochs,"+");                   
	}
	bFirstLoop[iChamb] = kFALSE;
      }
      else{    
	if(fDebugLevel>=2) fNet[iChamb] -> Train(nEpochs,"text update=10, graph+");      
	else fNet[iChamb] -> Train(nEpochs,"+");                   
      }
      
      // save weights for monitoring of the training
      fNet[iChamb] -> DumpWeights(Form("./Networks_%d/MomBin_%d/Net%d_%d",fDate, mombin, iChamb, iLoop));
    } // end chamber loop
  }   // end training loop
}


//________________________________________________________________________
void AliTRDpidRefMakerNN::MonitorTraining(Int_t mombin) 
{
  //
  // train the neural networks
  //
  
  if(!fContainer) LoadContainer(Form("TRD.Task%s.root", GetName()));

  if(!fContainer){
    AliError("ERROR container not available");
    return;
  }

  if (!fData) LoadFile(Form("TRD.Task%s.root", GetName()));
  if (!fData) {
    AliError("Tree for training list not available");
    return;
  }

  // init networks and set event list
  for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
	fNet[iChamb] = new TMultiLayerPerceptron("fdEdx[0],fdEdx[1],fdEdx[2],fdEdx[3],fdEdx[4],fdEdx[5],fdEdx[6],fdEdx[7]:15:7:fPID[0],fPID[1],fPID[2],fPID[3],fPID[4]!",fData,fTrain[mombin][iChamb],fTest[mombin][iChamb]);   
	fData -> SetEventList(fTrain[mombin][iChamb]);
	fData -> SetEventList(fTest[mombin][iChamb]);
  }

  // implement variables for likelihoods
  Float_t Like[AliPID::kSPECIES][AliTRDgeometry::kNlayer];
  memset(Like, 0, AliPID::kSPECIES*AliTRDgeometry::kNlayer*sizeof(Float_t));
  Float_t LikeAll[AliPID::kSPECIES], TotProb;

  Double_t PionEffiTrain[kMoniTrain], PionEffiErrTrain[kMoniTrain];
  Double_t PionEffiTest[kMoniTrain], PionEffiErrTest[kMoniTrain];
  memset(PionEffiTrain, 0, kMoniTrain*sizeof(Double_t));
  memset(PionEffiErrTrain, 0, kMoniTrain*sizeof(Double_t));
  memset(PionEffiTest, 0, kMoniTrain*sizeof(Double_t));
  memset(PionEffiErrTest, 0, kMoniTrain*sizeof(Double_t));

  // init histos
  const Float_t epsilon = 1/(2*(AliTRDpidUtil::kBins-1));     // get nice histos with bin center at 0 and 1
  TH1F *hElecs, *hPions;
  hElecs = new TH1F("hElecs","Likelihood for electrons", AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon);
  hPions = new TH1F("hPions","Likelihood for pions", AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon);

  TGraphErrors *gEffisTrain=0x0, *gEffisTest=0x0;
  gEffisTrain = (TGraphErrors*)fContainer->At(kGraphTrain);
  gEffisTest = (TGraphErrors*)fContainer->At(kGraphTest);

  AliTRDpidUtil *util = new AliTRDpidUtil();
  
  // monitor training progress
  for(Int_t iLoop = 0; iLoop < kMoniTrain; iLoop++){

    // load weights
    for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
      fNet[iChamb] -> LoadWeights(Form("./Networks_%d/MomBin_%d/Net%d_%d",fDate, mombin, iChamb, iLoop));
    }

    // event loop training list
    for(Int_t iEvent = 0; iEvent < fTrain[mombin][0] -> GetN(); iEvent++ ){

      // reset particle probabilities
      for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
	LikeAll[iPart] = 1./AliPID::kSPECIES;
      }
      TotProb = 0.;

      fData -> GetEntry(fTrain[mombin][0] -> GetEntry(iEvent));
      // use event only if it is electron or pion
      if(!((fPID[AliPID::kElectron] == 1.0) || (fPID[AliPID::kPion] == 1.0))) continue;

      // get the probabilities for each particle type in each chamber
      for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
	for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
	  Like[iPart][iChamb] = fNet[iChamb] -> Result(fTrain[mombin][iChamb] -> GetEntry(iEvent), iPart);
	  LikeAll[iPart] *=  Like[iPart][iChamb];
	}
      }

      // get total probability and normalize it
      for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
	TotProb += LikeAll[iPart];
      }
      for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
	LikeAll[iPart] /= TotProb;
      }

      // fill likelihood distributions
      if(fPID[AliPID::kElectron] == 1)      
	hElecs -> Fill(LikeAll[AliPID::kElectron]);
      if(fPID[AliPID::kPion] == 1)      
	hPions -> Fill(LikeAll[AliPID::kElectron]);
    } // end event loop


    // calculate the pion efficiency and fill the graph
    util -> CalculatePionEffi(hElecs, hPions);
    PionEffiTrain[iLoop] = util -> GetPionEfficiency();
    PionEffiErrTrain[iLoop] = util -> GetError();

    gEffisTrain -> SetPoint(iLoop, iLoop+1, PionEffiTrain[iLoop]);
    gEffisTrain -> SetPointError(iLoop, 0, PionEffiErrTrain[iLoop]);
    hElecs -> Reset();
    hPions -> Reset();
    if(fDebugLevel>=2) Printf("TrainingLoop[%d] PionEfficiency[%f +/- %f]", iLoop, PionEffiTrain[iLoop], PionEffiErrTrain[iLoop]);
    // end training loop
    

    // event loop test list
    for(Int_t iEvent = 0; iEvent < fTest[mombin][0] -> GetN(); iEvent++ ){

      // reset particle probabilities
      for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
	LikeAll[iPart] = 1./AliTRDgeometry::kNlayer;
      }
      TotProb = 0.;

      fData -> GetEntry(fTest[mombin][0] -> GetEntry(iEvent));
      // use event only if it is electron or pion
      if(!((fPID[AliPID::kElectron] == 1.0) || (fPID[AliPID::kPion] == 1.0))) continue;

      // get the probabilities for each particle type in each chamber
      for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
	for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
	  Like[iPart][iChamb] = fNet[iChamb] -> Result(fTest[mombin][iChamb] -> GetEntry(iEvent), iPart);
	  LikeAll[iPart] *=  Like[iPart][iChamb];
	}
      }

      // get total probability and normalize it
      for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
	TotProb += LikeAll[iPart];
      }
      for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
	LikeAll[iPart] /= TotProb;
      }

      // fill likelihood distributions
      if(fPID[AliPID::kElectron] == 1)      
	hElecs -> Fill(LikeAll[AliPID::kElectron]);
      if(fPID[AliPID::kPion] == 1)      
	hPions -> Fill(LikeAll[AliPID::kElectron]);
    } // end event loop

    // calculate the pion efficiency and fill the graph
    util -> CalculatePionEffi(hElecs, hPions);
    PionEffiTest[iLoop] = util -> GetPionEfficiency();
    PionEffiErrTest[iLoop] = util -> GetError();

    gEffisTest -> SetPoint(iLoop, iLoop+1, PionEffiTest[iLoop]);
    gEffisTest -> SetPointError(iLoop, 0, PionEffiErrTest[iLoop]);
    hElecs -> Reset();
    hPions -> Reset();
    if(fDebugLevel>=2) Printf("TestLoop[%d] PionEfficiency[%f +/- %f] \n", iLoop, PionEffiTest[iLoop], PionEffiErrTest[iLoop]);
    
  } //   end training loop
 
  util -> Delete();

  gEffisTest -> Draw("PAL");
  gEffisTrain -> Draw("PL");

}
Commit	Line	Data
72f5f463	1	#include "TSystem.h"
	2	#include "TDatime.h"
	3	#include "TEventList.h"
	4	#include "TGraphErrors.h"
	5	#include "TTreeStream.h"
	6	#include "TH1F.h"
	7	#include "TMultiLayerPerceptron.h"
	8
	9	#include "AliLog.h"
	10	#include "AliPID.h"
	11	#include "AliESDEvent.h"
	12
	13	#include "../Cal/AliTRDCalPIDNN.h"
	14
	15	#include "AliTRDpidUtil.h"
	16
	17	#include "AliTRDpidRefMakerNN.h"
	18	#include "info/AliTRDtrackInfo.h"
	19
	20	// builds the reference tree for the training of neural networks
	21
	22
	23	ClassImp(AliTRDpidRefMakerNN)
	24
	25	//________________________________________________________________________
	26	AliTRDpidRefMakerNN::AliTRDpidRefMakerNN()
	27	:AliTRDpidRefMaker(UChar_t(AliTRDpidUtil::kNNslices), "PidRefMakerNN", "PID(NN) Reference Maker")
	28	,fTrainMomBin(kAll)
	29	,fEpochs(1000)
	30	,fMinTrain(100)
	31	,fDate(0)
	32	,fDoTraining(0)
	33	,fContinueTraining(0)
	34	,fTrainPath(0x0)
	35	{
	36	//
	37	// Default constructor
	38	//
	39	SetAbundance(.67, .33);
	40	SetScaledEdx(Float_t(AliTRDCalPIDNN::kMLPscale));
	41	TDatime datime;
	42	fDate = datime.GetDate();
	43	}
	44
	45
	46	//________________________________________________________________________
	47	AliTRDpidRefMakerNN::~AliTRDpidRefMakerNN()
	48	{
	49	}
	50
	51
	52	//________________________________________________________________________
	53	void AliTRDpidRefMakerNN::CreateOutputObjects()
	54	{
	55	// Create histograms
	56	// Called once
	57	AliTRDpidRefMaker::CreateOutputObjects();
	58	TGraphErrors *gEffisTrain = new TGraphErrors(kMoniTrain);
	59	gEffisTrain -> SetLineColor(4);
	60	gEffisTrain -> SetMarkerColor(4);
	61	gEffisTrain -> SetMarkerStyle(29);
	62	gEffisTrain -> SetMarkerSize(1);
	63
	64	TGraphErrors *gEffisTest = new TGraphErrors(kMoniTrain);
65	gEffisTest -> SetLineColor(2);
66	gEffisTest -> SetMarkerColor(2);
67	gEffisTest -> SetMarkerStyle(29);
68	gEffisTest -> SetMarkerSize(1);
69
70	fContainer -> AddAt(gEffisTrain,kGraphTrain);
71	fContainer -> AddAt(gEffisTest,kGraphTest);
72	}
73
74
75	//________________________________________________________________________
76	Bool_t AliTRDpidRefMakerNN::PostProcess()
77	{
78	// Draw result to the screen
79	// Called once at the end of the query
80
81	// build the training andthe test list for the neural networks
82	MakeTrainingList();
83	if(!fDoTraining) return kTRUE;
84
85	// train the neural networks and build the refrence histos for 2-dim LQ
86	gSystem->Exec(Form("mkdir ./Networks_%d/",fDate));
87	AliDebug(3, Form("TrainMomBin [%d] [%d]", fTrainMomBin, kAll));
88
89	// train single network for a single momentum (recommended)
90	if(!(fTrainMomBin == kAll)){
91	if(fTrain[fTrainMomBin][0] -> GetN() < fMinTrain){
92	AliWarning("Not enough events for training available! Please check Data sample!");
93	return kFALSE;
94	}
95	MakeRefs(fTrainMomBin);
96	MonitorTraining(fTrainMomBin);
97	}
98	// train all momenta
99	else{
100	for(Int_t iMomBin = 0; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
101	if(fTrain[iMomBin][0] -> GetN() < fMinTrain){
102	AliWarning(Form("Not enough events for training available for momentum bin [%d]! Please check Data sample!", iMomBin));
103	continue;
104	}
105	MakeRefs(fTrainMomBin);
106	MonitorTraining(iMomBin);
107	}
108	}
109
110	return kTRUE; // testing protection
111	}
112
113
114
115	//________________________________________________________________________
116	void AliTRDpidRefMakerNN::MakeRefs(Int_t mombin)
117	{
118	//
119	// train the neural networks
120	//
121
122
123	if (!fData) LoadFile(Form("TRD.%s.root", GetName()));
124
125	if (!fData) {
126	AliError("Tree for training list not available");
127	return;
128	}
129
130	TDatime datime;
131	fDate = datime.GetDate();
132
133	AliDebug(2, Form("Training momentum bin %d", mombin));
134
135	// set variable to monitor the training and to save the development of the networks
136	Int_t nEpochs = fEpochs/kMoniTrain;
137	AliDebug(2, Form("Training %d times %d epochs", kMoniTrain, nEpochs));
138
139	// make directories to save the networks
140	gSystem->Exec(Form("rm -r ./Networks_%d/MomBin_%d",fDate, mombin));
141	gSystem->Exec(Form("mkdir ./Networks_%d/MomBin_%d",fDate, mombin));
142
143	// variable to check if network can load weights from previous training
144	Bool_t bFirstLoop[AliTRDgeometry::kNlayer];
145	memset(bFirstLoop, kTRUE, AliTRDgeometry::kNlayer*sizeof(Bool_t));
146
147	// train networks over several loops and save them after each loop
148	for(Int_t iLoop = 0; iLoop < kMoniTrain; iLoop++){
149	// loop over chambers
150	for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
151	// set the event lists
152	fData -> SetEventList(fTrain[mombin][iChamb]);
153	fData -> SetEventList(fTest[mombin][iChamb]);
154
155	AliDebug(2, Form("Trainingloop[%d] Chamber[%d]", iLoop, iChamb));
156
157	// check if network is already implemented
158	if(bFirstLoop[iChamb] == kTRUE){
159	fNet[iChamb] = new TMultiLayerPerceptron("fdEdx[0],fdEdx[1],fdEdx[2],fdEdx[3],fdEdx[4],fdEdx[5],fdEdx[6],fdEdx[7]:15:7:fPID[0],fPID[1],fPID[2],fPID[3],fPID[4]!",fData,fTrain[mombin][iChamb],fTest[mombin][iChamb]);
160	fNet[iChamb] -> SetLearningMethod(TMultiLayerPerceptron::kStochastic); // set learning method
161	fNet[iChamb] -> TMultiLayerPerceptron::SetEta(0.001); // set learning speed
162	if(!fContinueTraining){
163	if(fDebugLevel>=2) fNet[iChamb] -> Train(nEpochs,"text update=10, graph");
164	else fNet[iChamb] -> Train(nEpochs,"");
165	}
166	else{
167	fNet[iChamb] -> LoadWeights(Form("./Networks_%d/MomBin_%d/Net%d_%d",fTrainPath, mombin, iChamb, kMoniTrain - 1));
168	if(fDebugLevel>=2) fNet[iChamb] -> Train(nEpochs,"text update=10, graph+");
169	else fNet[iChamb] -> Train(nEpochs,"+");
170	}
171	bFirstLoop[iChamb] = kFALSE;
172	}
173	else{
174	if(fDebugLevel>=2) fNet[iChamb] -> Train(nEpochs,"text update=10, graph+");
175	else fNet[iChamb] -> Train(nEpochs,"+");
176	}
177
178	// save weights for monitoring of the training
179	fNet[iChamb] -> DumpWeights(Form("./Networks_%d/MomBin_%d/Net%d_%d",fDate, mombin, iChamb, iLoop));
180	} // end chamber loop
181	} // end training loop
182	}
183
184
185	//________________________________________________________________________
186	void AliTRDpidRefMakerNN::MonitorTraining(Int_t mombin)
187	{
188	//
189	// train the neural networks
190	//
191
192	if(!fContainer) LoadContainer(Form("TRD.Task%s.root", GetName()));
193
194	if(!fContainer){
195	AliError("ERROR container not available");
196	return;
197	}
198
199	if (!fData) LoadFile(Form("TRD.Task%s.root", GetName()));
200	if (!fData) {
201	AliError("Tree for training list not available");
202	return;
203	}
204
205	// init networks and set event list
206	for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
207	fNet[iChamb] = new TMultiLayerPerceptron("fdEdx[0],fdEdx[1],fdEdx[2],fdEdx[3],fdEdx[4],fdEdx[5],fdEdx[6],fdEdx[7]:15:7:fPID[0],fPID[1],fPID[2],fPID[3],fPID[4]!",fData,fTrain[mombin][iChamb],fTest[mombin][iChamb]);
208	fData -> SetEventList(fTrain[mombin][iChamb]);
209	fData -> SetEventList(fTest[mombin][iChamb]);
210	}
211
212	// implement variables for likelihoods
213	Float_t Like[AliPID::kSPECIES][AliTRDgeometry::kNlayer];
214	memset(Like, 0, AliPID::kSPECIESAliTRDgeometry::kNlayersizeof(Float_t));
215	Float_t LikeAll[AliPID::kSPECIES], TotProb;
216
217	Double_t PionEffiTrain[kMoniTrain], PionEffiErrTrain[kMoniTrain];
218	Double_t PionEffiTest[kMoniTrain], PionEffiErrTest[kMoniTrain];
219	memset(PionEffiTrain, 0, kMoniTrain*sizeof(Double_t));
220	memset(PionEffiErrTrain, 0, kMoniTrain*sizeof(Double_t));
221	memset(PionEffiTest, 0, kMoniTrain*sizeof(Double_t));
222	memset(PionEffiErrTest, 0, kMoniTrain*sizeof(Double_t));
223
224	// init histos
225	const Float_t epsilon = 1/(2*(AliTRDpidUtil::kBins-1)); // get nice histos with bin center at 0 and 1
226	TH1F hElecs, hPions;
227	hElecs = new TH1F("hElecs","Likelihood for electrons", AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon);
228	hPions = new TH1F("hPions","Likelihood for pions", AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon);
229
230	TGraphErrors gEffisTrain=0x0, gEffisTest=0x0;
231	gEffisTrain = (TGraphErrors*)fContainer->At(kGraphTrain);
232	gEffisTest = (TGraphErrors*)fContainer->At(kGraphTest);
233
234	AliTRDpidUtil *util = new AliTRDpidUtil();
235
236	// monitor training progress
237	for(Int_t iLoop = 0; iLoop < kMoniTrain; iLoop++){
238
239	// load weights
240	for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
241	fNet[iChamb] -> LoadWeights(Form("./Networks_%d/MomBin_%d/Net%d_%d",fDate, mombin, iChamb, iLoop));
242	}
243
244	// event loop training list
245	for(Int_t iEvent = 0; iEvent < fTrain[mombin][0] -> GetN(); iEvent++ ){
246
247	// reset particle probabilities
248	for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
249	LikeAll[iPart] = 1./AliPID::kSPECIES;
250	}
251	TotProb = 0.;
252
253	fData -> GetEntry(fTrain[mombin][0] -> GetEntry(iEvent));
254	// use event only if it is electron or pion
255	if(!((fPID[AliPID::kElectron] == 1.0) \|\| (fPID[AliPID::kPion] == 1.0))) continue;
256
257	// get the probabilities for each particle type in each chamber
258	for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
259	for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
260	Like[iPart][iChamb] = fNet[iChamb] -> Result(fTrain[mombin][iChamb] -> GetEntry(iEvent), iPart);
261	LikeAll[iPart] *= Like[iPart][iChamb];
262	}
263	}
264
265	// get total probability and normalize it
266	for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
267	TotProb += LikeAll[iPart];
268	}
269	for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
270	LikeAll[iPart] /= TotProb;
271	}
272
273	// fill likelihood distributions
274	if(fPID[AliPID::kElectron] == 1)
275	hElecs -> Fill(LikeAll[AliPID::kElectron]);
276	if(fPID[AliPID::kPion] == 1)
277	hPions -> Fill(LikeAll[AliPID::kElectron]);
278	} // end event loop
279
280
281	// calculate the pion efficiency and fill the graph
282	util -> CalculatePionEffi(hElecs, hPions);
283	PionEffiTrain[iLoop] = util -> GetPionEfficiency();
284	PionEffiErrTrain[iLoop] = util -> GetError();
285
286	gEffisTrain -> SetPoint(iLoop, iLoop+1, PionEffiTrain[iLoop]);
287	gEffisTrain -> SetPointError(iLoop, 0, PionEffiErrTrain[iLoop]);
288	hElecs -> Reset();
289	hPions -> Reset();
290	if(fDebugLevel>=2) Printf("TrainingLoop[%d] PionEfficiency[%f +/- %f]", iLoop, PionEffiTrain[iLoop], PionEffiErrTrain[iLoop]);
291	// end training loop
292
293
294
295	// event loop test list
296	for(Int_t iEvent = 0; iEvent < fTest[mombin][0] -> GetN(); iEvent++ ){
297
298	// reset particle probabilities
299	for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
300	LikeAll[iPart] = 1./AliTRDgeometry::kNlayer;
301	}
302	TotProb = 0.;
303
304	fData -> GetEntry(fTest[mombin][0] -> GetEntry(iEvent));
305	// use event only if it is electron or pion
306	if(!((fPID[AliPID::kElectron] == 1.0) \|\| (fPID[AliPID::kPion] == 1.0))) continue;
307
308	// get the probabilities for each particle type in each chamber
309	for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
310	for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
311	Like[iPart][iChamb] = fNet[iChamb] -> Result(fTest[mombin][iChamb] -> GetEntry(iEvent), iPart);
312	LikeAll[iPart] *= Like[iPart][iChamb];
313	}
314	}
315
316	// get total probability and normalize it
317	for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
318	TotProb += LikeAll[iPart];
319	}
320	for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
321	LikeAll[iPart] /= TotProb;
322	}
323
324	// fill likelihood distributions
325	if(fPID[AliPID::kElectron] == 1)
326	hElecs -> Fill(LikeAll[AliPID::kElectron]);
327	if(fPID[AliPID::kPion] == 1)
328	hPions -> Fill(LikeAll[AliPID::kElectron]);
329	} // end event loop
330
331	// calculate the pion efficiency and fill the graph
332	util -> CalculatePionEffi(hElecs, hPions);
333	PionEffiTest[iLoop] = util -> GetPionEfficiency();
334	PionEffiErrTest[iLoop] = util -> GetError();
335
336	gEffisTest -> SetPoint(iLoop, iLoop+1, PionEffiTest[iLoop]);
337	gEffisTest -> SetPointError(iLoop, 0, PionEffiErrTest[iLoop]);
338	hElecs -> Reset();
339	hPions -> Reset();
340	if(fDebugLevel>=2) Printf("TestLoop[%d] PionEfficiency[%f +/- %f] \n", iLoop, PionEffiTest[iLoop], PionEffiErrTest[iLoop]);
341
342	} // end training loop
343
344	util -> Delete();
345
346	gEffisTest -> Draw("PAL");
347	gEffisTrain -> Draw("PL");
348
349	}
350
351