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