]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/qaRec/AliTRDpidRefMaker.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
updates of the PID reference data calibration (Alex W)
[u/mrichter/AliRoot.git] / TRD / qaRec / AliTRDpidRefMaker.cxx
CommitLineData
aea8dd24 1#include "TSystem.h"
2#include "TDatime.h"
28efdace 3#include "TPDGCode.h"
4#include "TH1F.h"
aea8dd24 5#include "TFile.h"
82719f10 6#include "TGraphErrors.h"
aea8dd24 7#include "TTree.h"
28efdace 8#include "TTreeStream.h"
aea8dd24 9#include "TEventList.h"
10#include "TMultiLayerPerceptron.h"
28efdace 11
82719f10 12#include "AliPID.h"
28efdace 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
82719f10 24#include "AliTRDpidUtil.h"
25
28efdace 26#include "AliTRDpidRefMaker.h"
27#include "AliTRDtrackInfo/AliTRDtrackInfo.h"
28
29// builds the reference tree for the training of neural networks
30
31
32ClassImp(AliTRDpidRefMaker)
33
34//________________________________________________________________________
35AliTRDpidRefMaker::AliTRDpidRefMaker()
aea8dd24 36 :AliTRDrecoTask("PidRefMaker", "PID(NN) Reference Maker")
28efdace 37 ,fReconstructor(0x0)
aea8dd24 38 ,fNN(0x0)
39 ,fLQ(0x0)
40 ,fLayer(0xff)
82719f10 41 ,fTrainMomBin(kAll)
42 ,fEpochs(1000)
43 ,fMinTrain(100)
aea8dd24 44 ,fDate(0)
45 ,fMom(0.)
82719f10 46 ,fDoTraining(0)
b46633bb 47 ,fContinueTraining(0)
48 ,fTrainPath(0x0)
28efdace 49{
50 //
51 // Default constructor
52 //
53
54 fReconstructor = new AliTRDReconstructor();
55 fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
aea8dd24 56 memset(fv0pid, 0, AliPID::kSPECIES*sizeof(Float_t));
57 memset(fdEdx, 0, 10*sizeof(Float_t));
58
5d6dc395 59 const Int_t nnSize = AliTRDCalPID::kNMom * AliTRDgeometry::kNlayer;
aea8dd24 60 memset(fTrain, 0, nnSize*sizeof(TEventList*));
61 memset(fTest, 0, nnSize*sizeof(TEventList*));
5d6dc395 62 memset(fNet, 0, AliTRDgeometry::kNlayer*sizeof(TMultiLayerPerceptron*));
aea8dd24 63
aea8dd24 64 TDatime datime;
65 fDate = datime.GetDate();
66
82719f10 67 DefineOutput(1, TTree::Class());
68 DefineOutput(2, TTree::Class());
28efdace 69}
70
71
72//________________________________________________________________________
73AliTRDpidRefMaker::~AliTRDpidRefMaker()
74{
75 if(fReconstructor) delete fReconstructor;
b718144c 76 //if(fNN) delete fNN;
77 //if(fLQ) delete fLQ;
28efdace 78}
79
80
81//________________________________________________________________________
82void AliTRDpidRefMaker::CreateOutputObjects()
83{
84 // Create histograms
85 // Called once
86
87 OpenFile(0, "RECREATE");
88 fContainer = new TObjArray();
28efdace 89 fContainer->AddAt(new TH1F("hPDG","hPDG",AliPID::kSPECIES,-0.5,5.5),0);
aea8dd24 90
82719f10 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
aea8dd24 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));
0d83b3a5 112 fNN->Branch("fdEdx", fdEdx, Form("fdEdx[%d]/F", AliTRDpidUtil::kNNslices));
aea8dd24 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));
0d83b3a5 120 fLQ->Branch("fdEdx", fdEdx, Form("fdEdx[%d]/F", AliTRDpidUtil::kLQslices));
28efdace 121}
122
123
124//________________________________________________________________________
125void 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
28efdace 133 Float_t mom;
134 ULong_t status;
135 Int_t nTRD = 0;
28efdace 136
137 AliTRDtrackInfo *track = 0x0;
138 AliTRDtrackV1 *TRDtrack = 0x0;
139 AliTrackReference *ref = 0x0;
140 AliExternalTrackParam *esd = 0x0;
141
aea8dd24 142 AliTRDseedV1 *TRDtracklet = 0x0;
28efdace 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++)
aea8dd24 150 fv0pid[iPart] = 0.;
28efdace 151
152 track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
153 if(!track->HasESDtrack()) continue;
154 status = track->GetStatus();
155 if(!(status&AliESDtrack::kTPCout)) continue;
156
b718144c 157 if(!(TRDtrack = track->GetTrack())) continue;
28efdace 158 //&&(track->GetNumberOfClustersRefit()
159
160 // use only tracks that hit 6 chambers
5d6dc395 161 if(!(TRDtrack->GetNumberOfTracklets() == AliTRDgeometry::kNlayer)) continue;
28efdace 162
163 ref = track->GetTrackRef(0);
b38d0daf 164 esd = track->GetESDinfo()->GetOuterParam();
28efdace 165 mom = ref ? ref->P(): esd->P();
aea8dd24 166 fMom = mom;
167
28efdace 168
169 labelsacc[nTRD] = track->GetLabel();
170 nTRD++;
171
172 // if no monte carlo data available -> use V0 information
173 if(!HasMCdata()){
aea8dd24 174 GetV0info(TRDtrack,fv0pid);
28efdace 175 }
176 // else use the MC info
177 else{
178 switch(track -> GetPDG()){
179 case kElectron:
180 case kPositron:
aea8dd24 181 fv0pid[AliPID::kElectron] = 1.;
28efdace 182 break;
183 case kMuonPlus:
184 case kMuonMinus:
aea8dd24 185 fv0pid[AliPID::kMuon] = 1.;
28efdace 186 break;
187 case kPiPlus:
188 case kPiMinus:
aea8dd24 189 fv0pid[AliPID::kPion] = 1.;
28efdace 190 break;
191 case kKPlus:
192 case kKMinus:
aea8dd24 193 fv0pid[AliPID::kKaon] = 1.;
28efdace 194 break;
195 case kProton:
196 case kProtonBar:
aea8dd24 197 fv0pid[AliPID::kProton] = 1.;
28efdace 198 break;
199 }
200 }
201
28efdace 202 // set reconstructor
aea8dd24 203 Float_t *dedx;
28efdace 204 TRDtrack -> SetReconstructor(fReconstructor);
28efdace 205
aea8dd24 206 // fill the dE/dx information for NN
207 fReconstructor -> SetOption("nn");
5d6dc395 208 for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){
aea8dd24 209 if(!(TRDtracklet = TRDtrack -> GetTracklet(ily))) continue;
0d83b3a5 210 TRDtracklet->CookdEdx(AliTRDpidUtil::kNNslices);
aea8dd24 211 dedx = TRDtracklet->GetdEdx();
0d83b3a5 212 for(Int_t iSlice = 0; iSlice < AliTRDpidUtil::kNNslices; iSlice++)
aea8dd24 213 dedx[iSlice] = dedx[iSlice]/AliTRDCalPIDNN::kMLPscale;
0d83b3a5 214 memcpy(fdEdx, dedx, AliTRDpidUtil::kNNslices*sizeof(Float_t));
aea8dd24 215 if(fDebugLevel>=2) Printf("LayerNN : %d", ily);
216 fLayer = ily;
217 fNN->Fill();
28efdace 218 }
219
220
aea8dd24 221 // fill the dE/dx information for LQ
222 fReconstructor -> SetOption("!nn");
5d6dc395 223 for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){
aea8dd24 224 if(!(TRDtracklet = TRDtrack -> GetTracklet(ily))) continue;
0d83b3a5 225 TRDtracklet->CookdEdx(AliTRDpidUtil::kLQslices);
aea8dd24 226 dedx = TRDtracklet->GetdEdx();
0d83b3a5 227 memcpy(fdEdx, dedx, AliTRDpidUtil::kLQslices*sizeof(Float_t));
aea8dd24 228 if(fDebugLevel>=2) Printf("LayerLQ : %d", ily);
229 fLayer = ily;
230 fLQ->Fill();
28efdace 231 }
aea8dd24 232
28efdace 233
234 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
aea8dd24 235 if(fDebugLevel>=4) Printf("PDG is %d %f", iPart, fv0pid[iPart]);
28efdace 236 }
237 }
238
239 PostData(0, fContainer);
aea8dd24 240 PostData(1, fNN);
241 PostData(2, fLQ);
28efdace 242}
243
244
28efdace 245//________________________________________________________________________
246Bool_t AliTRDpidRefMaker::PostProcess()
247{
248 // Draw result to the screen
249 // Called once at the end of the query
250
aea8dd24 251 // build the training andthe test list for the neural networks
252 MakeTrainingLists();
82719f10 253 if(!fDoTraining) return kTRUE;
aea8dd24 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);
82719f10 267 MonitorTraining(fTrainMomBin);
aea8dd24 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);
82719f10 278 MonitorTraining(iMomBin);
aea8dd24 279 }
280 }
281
28efdace 282 return kTRUE; // testing protection
283}
284
285
286//________________________________________________________________________
287void 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//________________________________________________________________________
aea8dd24 301void AliTRDpidRefMaker::GetV0info(AliTRDtrackV1 *TRDtrack, Float_t *v0pid)
28efdace 302{
28efdace 303 // !!!! PREMILMINARY FUNCTION !!!!
304 //
305 // this is the place for the V0 procedure
aea8dd24 306 // as long as there is no one implemented,
307 // just the probabilities
308 // of the TRDtrack are used!
28efdace 309
310 TRDtrack -> SetReconstructor(fReconstructor);
311 fReconstructor -> SetOption("nn");
312 TRDtrack -> CookPID();
313 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
aea8dd24 314 v0pid[iPart] = TRDtrack -> GetPID(iPart);
315 if(fDebugLevel>=4) Printf("PDG is (in V0info) %d %f", iPart, v0pid[iPart]);
28efdace 316 }
317}
aea8dd24 318
319
320//________________________________________________________________________
321void AliTRDpidRefMaker::MakeTrainingLists()
322{
323 //
324 // build the training lists for the neural networks
325 //
326
82719f10 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
aea8dd24 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
82719f10 346 AliTRDpidUtil *util = new AliTRDpidUtil();
347
aea8dd24 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;
82719f10 358 iMomBin = util -> GetMomentumBin(fMom);
aea8dd24 359
360 // check PID information and count particle types per momentum interval
82719f10 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]++;
aea8dd24 371 }
372
373 if(fDebugLevel>=2){
374 Printf("Particle multiplicities:");
375 for(Int_t iMomBin = 0; iMomBin <AliTRDCalPID::kNMom; iMomBin++)
82719f10 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]);
aea8dd24 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;
82719f10 413 iMomBin = util -> GetMomentumBin(fMom);
aea8dd24 414
415 // set electrons
82719f10 416 if(fv0pid[AliPID::kElectron] == 1){
417 if(nPart[AliPID::kElectron][iMomBin] < iTrain[iMomBin]){
5d6dc395 418 for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
aea8dd24 419 fTrain[iMomBin][ily] -> Enter(iEv + ily);
82719f10 420 nPart[AliPID::kElectron][iMomBin]++;
aea8dd24 421 }
82719f10 422 else if(nPart[AliPID::kElectron][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){
5d6dc395 423 for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
aea8dd24 424 fTest[iMomBin][ily] -> Enter(iEv + ily);
82719f10 425 nPart[AliPID::kElectron][iMomBin]++;
aea8dd24 426 }
427 else
428 continue;
429 }
430 // set muons
82719f10 431 else if(fv0pid[AliPID::kMuon] == 1){
432 if(nPart[AliPID::kMuon][iMomBin] < iTrain[iMomBin]){
5d6dc395 433 for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
aea8dd24 434 fTrain[iMomBin][ily] -> Enter(iEv + ily);
82719f10 435 nPart[AliPID::kMuon][iMomBin]++;
aea8dd24 436 }
82719f10 437 else if(nPart[AliPID::kMuon][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){
5d6dc395 438 for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
aea8dd24 439 fTest[iMomBin][ily] -> Enter(iEv + ily);
82719f10 440 nPart[AliPID::kMuon][iMomBin]++;
aea8dd24 441 }
442 else
443 continue;
444 }
445 // set pions
82719f10 446 else if(fv0pid[AliPID::kPion] == 1){
447 if(nPart[AliPID::kPion][iMomBin] < iTrain[iMomBin]){
5d6dc395 448 for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
aea8dd24 449 fTrain[iMomBin][ily] -> Enter(iEv + ily);
82719f10 450 nPart[AliPID::kPion][iMomBin]++;
aea8dd24 451 }
82719f10 452 else if(nPart[AliPID::kPion][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){
5d6dc395 453 for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
aea8dd24 454 fTest[iMomBin][ily] -> Enter(iEv + ily);
82719f10 455 nPart[AliPID::kPion][iMomBin]++;
aea8dd24 456 }
457 else
458 continue;
459 }
460 // set kaons
82719f10 461 else if(fv0pid[AliPID::kKaon] == 1){
462 if(nPart[AliPID::kKaon][iMomBin] < iTrain[iMomBin]){
5d6dc395 463 for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
aea8dd24 464 fTrain[iMomBin][ily] -> Enter(iEv + ily);
82719f10 465 nPart[AliPID::kKaon][iMomBin]++;
aea8dd24 466 }
82719f10 467 else if(nPart[AliPID::kKaon][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){
5d6dc395 468 for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
aea8dd24 469 fTest[iMomBin][ily] -> Enter(iEv + ily);
82719f10 470 nPart[AliPID::kKaon][iMomBin]++;
aea8dd24 471 }
472 else
473 continue;
474 }
475 // set protons
82719f10 476 else if(fv0pid[AliPID::kProton] == 1){
477 if(nPart[AliPID::kProton][iMomBin] < iTrain[iMomBin]){
5d6dc395 478 for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
aea8dd24 479 fTrain[iMomBin][ily] -> Enter(iEv + ily);
82719f10 480 nPart[AliPID::kProton][iMomBin]++;
aea8dd24 481 }
82719f10 482 else if(nPart[AliPID::kProton][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){
5d6dc395 483 for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++)
aea8dd24 484 fTest[iMomBin][ily] -> Enter(iEv + ily);
82719f10 485 nPart[AliPID::kProton][iMomBin]++;
aea8dd24 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++)
82719f10 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]);
aea8dd24 496 Printf("\n");
497 }
82719f10 498
499 util -> Delete();
aea8dd24 500}
501
502
503//________________________________________________________________________
504void AliTRDpidRefMaker::TrainNetworks(Int_t mombin)
505{
506 //
507 // train the neural networks
508 //
509
82719f10 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
aea8dd24 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
5d6dc395 534 Bool_t bFirstLoop[AliTRDgeometry::kNlayer];
535 memset(bFirstLoop, kTRUE, AliTRDgeometry::kNlayer*sizeof(Bool_t));
aea8dd24 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
5d6dc395 540 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
aea8dd24 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
b46633bb 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 }
aea8dd24 561 bFirstLoop[iChamb] = kFALSE;
562 }
563 else{
82719f10 564 if(fDebugLevel>=2) fNet[iChamb] -> Train(nEpochs,"text update=10, graph+");
aea8dd24 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//________________________________________________________________________
576void 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//________________________________________________________________________
82719f10 587void AliTRDpidRefMaker::MonitorTraining(Int_t mombin)
aea8dd24 588{
589 //
590 // train the neural networks
591 //
592
82719f10 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];
a391a274 618 memset(Like, 0, AliPID::kSPECIES*AliTRDgeometry::kNlayer*sizeof(Float_t));
82719f10 619 Float_t LikeAll[AliPID::kSPECIES], TotProb;
620
621 Double_t PionEffiTrain[kMoniTrain], PionEffiErrTrain[kMoniTrain];
622 Double_t PionEffiTest[kMoniTrain], PionEffiErrTest[kMoniTrain];
a391a274 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));
82719f10 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
82719f10 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
aea8dd24 753}
82719f10 754
755
756//________________________________________________________________________
757void 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//________________________________________________________________________
784void 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