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