6 #include "TGraphErrors.h"
8 #include "TTreeStream.h"
9 #include "TEventList.h"
10 #include "TMultiLayerPerceptron.h"
13 #include "AliESDEvent.h"
14 #include "AliESDInputHandler.h"
15 #include "AliTrackReference.h"
17 #include "AliAnalysisTask.h"
19 #include "AliTRDtrackV1.h"
20 #include "AliTRDReconstructor.h"
21 #include "../Cal/AliTRDCalPID.h"
22 #include "../Cal/AliTRDCalPIDNN.h"
24 #include "AliTRDpidUtil.h"
26 #include "AliTRDpidRefMaker.h"
27 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
29 // builds the reference tree for the training of neural networks
32 ClassImp(AliTRDpidRefMaker)
34 //________________________________________________________________________
35 AliTRDpidRefMaker::AliTRDpidRefMaker()
36 :AliTRDrecoTask("PidRefMaker", "PID(NN) Reference Maker")
49 // Default constructor
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));
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*));
63 fDate = datime.GetDate();
65 DefineOutput(1, TTree::Class());
66 DefineOutput(2, TTree::Class());
70 //________________________________________________________________________
71 AliTRDpidRefMaker::~AliTRDpidRefMaker()
73 if(fReconstructor) delete fReconstructor;
79 //________________________________________________________________________
80 void AliTRDpidRefMaker::CreateOutputObjects()
85 OpenFile(0, "RECREATE");
86 fContainer = new TObjArray();
87 fContainer->AddAt(new TH1F("hPDG","hPDG",AliPID::kSPECIES,-0.5,5.5),0);
89 TGraphErrors *gEffisTrain = new TGraphErrors(kMoniTrain);
90 gEffisTrain -> SetLineColor(4);
91 gEffisTrain -> SetMarkerColor(4);
92 gEffisTrain -> SetMarkerStyle(29);
93 gEffisTrain -> SetMarkerSize(1);
95 TGraphErrors *gEffisTest = new TGraphErrors(kMoniTrain);
96 gEffisTest -> SetLineColor(2);
97 gEffisTest -> SetMarkerColor(2);
98 gEffisTest -> SetMarkerStyle(29);
99 gEffisTest -> SetMarkerSize(1);
101 fContainer -> AddAt(gEffisTrain,kGraphTrain);
102 fContainer -> AddAt(gEffisTest,kGraphTest);
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));
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));
122 //________________________________________________________________________
123 void AliTRDpidRefMaker::Exec(Option_t *)
126 // Called for each event
128 Int_t labelsacc[10000];
129 memset(labelsacc, 0, sizeof(Int_t) * 10000);
135 AliTRDtrackInfo *track = 0x0;
136 AliTRDtrackV1 *TRDtrack = 0x0;
137 AliTrackReference *ref = 0x0;
138 AliExternalTrackParam *esd = 0x0;
140 AliTRDseedV1 *TRDtracklet = 0x0;
142 //AliTRDcluster *TRDcluster = 0x0;
144 for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
146 // reset the pid information
147 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++)
150 track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
151 if(!track->HasESDtrack()) continue;
152 status = track->GetStatus();
153 if(!(status&AliESDtrack::kTPCout)) continue;
155 if(!(TRDtrack = track->GetTrack())) continue;
156 //&&(track->GetNumberOfClustersRefit()
158 // use only tracks that hit 6 chambers
159 if(!(TRDtrack->GetNumberOfTracklets() == AliTRDgeometry::kNlayer)) continue;
161 ref = track->GetTrackRef(0);
162 esd = track->GetOuterParam();
163 mom = ref ? ref->P(): esd->P();
167 labelsacc[nTRD] = track->GetLabel();
170 // if no monte carlo data available -> use V0 information
172 GetV0info(TRDtrack,fv0pid);
174 // else use the MC info
176 switch(track -> GetPDG()){
179 fv0pid[AliPID::kElectron] = 1.;
183 fv0pid[AliPID::kMuon] = 1.;
187 fv0pid[AliPID::kPion] = 1.;
191 fv0pid[AliPID::kKaon] = 1.;
195 fv0pid[AliPID::kProton] = 1.;
202 TRDtrack -> SetReconstructor(fReconstructor);
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);
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);
232 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
233 if(fDebugLevel>=4) Printf("PDG is %d %f", iPart, fv0pid[iPart]);
237 PostData(0, fContainer);
243 //________________________________________________________
244 void AliTRDpidRefMaker::GetRefFigure(Int_t /*ifig*/)
250 //________________________________________________________________________
251 Bool_t AliTRDpidRefMaker::PostProcess()
253 // Draw result to the screen
254 // Called once at the end of the query
256 // build the training andthe test list for the neural networks
258 if(!fDoTraining) return kTRUE;
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);
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!");
270 TrainNetworks(fTrainMomBin);
271 BuildLQRefs(fTrainMomBin);
272 MonitorTraining(fTrainMomBin);
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);
281 TrainNetworks(iMomBin);
282 BuildLQRefs(fTrainMomBin);
283 MonitorTraining(iMomBin);
287 return kTRUE; // testing protection
291 //________________________________________________________________________
292 void AliTRDpidRefMaker::Terminate(Option_t *)
294 // Draw result to the screen
295 // Called once at the end of the query
297 fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
299 Printf("ERROR: list not available");
305 //________________________________________________________________________
306 void AliTRDpidRefMaker::GetV0info(AliTRDtrackV1 *TRDtrack, Float_t *v0pid)
308 // !!!! PREMILMINARY FUNCTION !!!!
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!
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]);
325 //________________________________________________________________________
326 void AliTRDpidRefMaker::MakeTrainingLists()
329 // build the training lists for the neural networks
333 LoadFiles("TRD.TaskPidRefMakerNN.root", "TRD.TaskPidRefMakerLQ.root");
337 Printf("ERROR tree for training list not available");
341 if(fDebugLevel>=2) Printf("\n Making training lists! \n");
343 Int_t nPart[AliPID::kSPECIES][AliTRDCalPID::kNMom];
344 memset(nPart, 0, AliPID::kSPECIES*AliTRDCalPID::kNMom*sizeof(Int_t));
346 // set needed branches
347 fNN -> SetBranchAddress("fv0pid", &fv0pid);
348 fNN -> SetBranchAddress("fMom", &fMom);
349 fNN -> SetBranchAddress("fLayer", &fLayer);
351 AliTRDpidUtil *util = new AliTRDpidUtil();
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);
357 // use only events with goes through 6 layers TRD
361 // set the 11 momentum bins
363 iMomBin = util -> GetMomentumBin(fMom);
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]++;
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]);
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));
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];
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]);
402 if(fDebugLevel>=2) Printf("\n");
406 memset(nPart, 0, AliPID::kSPECIES*AliTRDCalPID::kNMom*sizeof(Int_t));
408 // start second loop to set the event lists
409 for(Int_t iEv = 0; iEv < fNN -> GetEntries(); iEv++){
410 fNN -> GetEntry(iEv);
412 // use only events with goes through 6 layers TRD
416 // set the 11 momentum bins
418 iMomBin = util -> GetMomentumBin(fMom);
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]++;
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]++;
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]++;
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]++;
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]++;
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]++;
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]++;
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]++;
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]++;
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]++;
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]);
508 //________________________________________________________________________
509 void AliTRDpidRefMaker::TrainNetworks(Int_t mombin)
512 // train the neural networks
517 LoadFiles("TRD.TaskPidRefMakerNN.root", "TRD.TaskPidRefMakerLQ.root");
521 Printf("ERROR tree for training list not available");
526 fDate = datime.GetDate();
528 if(fDebugLevel>=2) Printf("Training momentum bin %d", mombin);
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);
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));
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));
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]);
550 if(fDebugLevel>=2) Printf("Trainingloop[%d] Chamber[%d]", iLoop, iChamb);
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;
562 if(fDebugLevel>=2) fNet[iChamb] -> Train(nEpochs,"text update=10, graph+");
563 else fNet[iChamb] -> Train(nEpochs,"+");
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
573 //________________________________________________________________________
574 void AliTRDpidRefMaker::BuildLQRefs(Int_t mombin)
577 // build the 2-dim LQ reference histograms
580 if(fDebugLevel>=2) Printf("Building LQRefs for momentum bin %d", mombin);
584 //________________________________________________________________________
585 void AliTRDpidRefMaker::MonitorTraining(Int_t mombin)
588 // train the neural networks
592 LoadContainer("TRD.TaskPidRefMaker.root");
595 Printf("ERROR container not available");
600 LoadFiles("TRD.TaskPidRefMakerNN.root", "TRD.TaskPidRefMakerLQ.root");
603 Printf("ERROR tree for training list not available");
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]);
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;
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));
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);
632 TGraphErrors *gEffisTrain=0x0, *gEffisTest=0x0;
633 gEffisTrain = (TGraphErrors*)fContainer->At(kGraphTrain);
634 gEffisTest = (TGraphErrors*)fContainer->At(kGraphTest);
636 AliTRDpidUtil *util = new AliTRDpidUtil();
638 // monitor training progress
639 for(Int_t iLoop = 0; iLoop < kMoniTrain; iLoop++){
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));
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));
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]++;
666 // Printf("Numbers [%d] [%d] [%d] [%d] [%d]", iType[AliPID::kElectron], iType[AliPID::kMuon], iType[AliPID::kPion], iType[AliPID::kKaon], iType[AliPID::kProton]);
670 // event loop training list
671 for(Int_t iEvent = 0; iEvent < fTrain[mombin][0] -> GetN(); iEvent++ ){
673 // reset particle probabilities
674 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
675 LikeAll[iPart] = 1./AliPID::kSPECIES;
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;
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];
691 // get total probability and normalize it
692 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
693 TotProb += LikeAll[iPart];
695 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
696 LikeAll[iPart] /= TotProb;
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]);
707 // calculate the pion efficiency and fill the graph
708 util -> CalculatePionEffi(hElecs, hPions);
709 PionEffiTrain[iLoop] = util -> GetPionEfficiency();
710 PionEffiErrTrain[iLoop] = util -> GetError();
712 gEffisTrain -> SetPoint(iLoop, iLoop+1, PionEffiTrain[iLoop]);
713 gEffisTrain -> SetPointError(iLoop, 0, PionEffiErrTrain[iLoop]);
716 if(fDebugLevel>=2) Printf("TrainingLoop[%d] PionEfficiency[%f +/- %f]", iLoop, PionEffiTrain[iLoop], PionEffiErrTrain[iLoop]);
721 // event loop test list
722 for(Int_t iEvent = 0; iEvent < fTest[mombin][0] -> GetN(); iEvent++ ){
724 // reset particle probabilities
725 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
726 LikeAll[iPart] = 1./AliTRDgeometry::kNlayer;
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;
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];
742 // get total probability and normalize it
743 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
744 TotProb += LikeAll[iPart];
746 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
747 LikeAll[iPart] /= TotProb;
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]);
757 // calculate the pion efficiency and fill the graph
758 util -> CalculatePionEffi(hElecs, hPions);
759 PionEffiTest[iLoop] = util -> GetPionEfficiency();
760 PionEffiErrTest[iLoop] = util -> GetError();
762 gEffisTest -> SetPoint(iLoop, iLoop+1, PionEffiTest[iLoop]);
763 gEffisTest -> SetPointError(iLoop, 0, PionEffiErrTest[iLoop]);
766 if(fDebugLevel>=2) Printf("TestLoop[%d] PionEfficiency[%f +/- %f] \n", iLoop, PionEffiTest[iLoop], PionEffiErrTest[iLoop]);
768 } // end training loop
772 gEffisTest -> Draw("PAL");
773 gEffisTrain -> Draw("PL");
778 //________________________________________________________________________
779 void AliTRDpidRefMaker::LoadFiles(const Char_t *InFileNN, const Char_t *InFileLQ)
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
789 fInFileNN = new TFile(InFileNN, "READ");
790 fNN = (TTree*)fInFileNN -> Get("NN");
793 fInFileLQ = new TFile(InFileLQ, "READ");
794 fLQ = (TTree*)fInFileLQ -> Get("LQ");
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));
805 //________________________________________________________________________
806 void AliTRDpidRefMaker::LoadContainer(const Char_t *InFileCont)
810 // Loads the container if no container is there.
811 // Useable for training outside of the makeResults.C macro
815 fInFileCont = new TFile(InFileCont, "READ");
816 fContainer = (TObjArray*)fInFileCont -> Get("PidRefMaker");
821 // //________________________________________________________________________
822 // void AliTRDpidRefMaker::CreateGraphs()
824 // // Create histograms
827 // OpenFile(0, "RECREATE");
828 // fContainer = new TObjArray();
829 // fContainer->AddAt(new TH1F("hPDG","hPDG",AliPID::kSPECIES,-0.5,5.5),0);
831 // TGraphErrors *gEffisTrain = new TGraphErrors(kMoniTrain);
832 // gEffisTrain -> SetLineColor(4);
833 // gEffisTrain -> SetMarkerColor(4);
834 // gEffisTrain -> SetMarkerStyle(29);
835 // gEffisTrain -> SetMarkerSize(2);
837 // TGraphErrors *gEffisTest = new TGraphErrors(kMoniTrain);
838 // gEffisTest -> SetLineColor(2);
839 // gEffisTest -> SetMarkerColor(2);
840 // gEffisTest -> SetMarkerSize(2);
842 // fContainer -> AddAt(gEffisTrain,kGraphTrain);
843 // fContainer -> AddAt(gEffisTest,kGraphTest);