]>
Commit | Line | Data |
---|---|---|
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 | ||
32 | ClassImp(AliTRDpidRefMaker) | |
33 | ||
34 | //________________________________________________________________________ | |
35 | AliTRDpidRefMaker::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 | //________________________________________________________________________ | |
71 | AliTRDpidRefMaker::~AliTRDpidRefMaker() | |
72 | { | |
73 | if(fReconstructor) delete fReconstructor; | |
b718144c | 74 | //if(fNN) delete fNN; |
75 | //if(fLQ) delete fLQ; | |
28efdace | 76 | } |
77 | ||
78 | ||
79 | //________________________________________________________________________ | |
80 | void 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 | //________________________________________________________________________ | |
123 | void 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); | |
162 | esd = track->GetOuterParam(); | |
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 | ||
243 | //________________________________________________________ | |
a391a274 | 244 | void AliTRDpidRefMaker::GetRefFigure(Int_t /*ifig*/) |
28efdace | 245 | { |
246 | ||
247 | } | |
248 | ||
249 | ||
250 | //________________________________________________________________________ | |
251 | Bool_t AliTRDpidRefMaker::PostProcess() | |
252 | { | |
253 | // Draw result to the screen | |
254 | // Called once at the end of the query | |
255 | ||
aea8dd24 | 256 | // build the training andthe test list for the neural networks |
257 | MakeTrainingLists(); | |
82719f10 | 258 | if(!fDoTraining) return kTRUE; |
aea8dd24 | 259 | |
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); | |
263 | ||
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!"); | |
268 | return kFALSE; | |
269 | } | |
270 | TrainNetworks(fTrainMomBin); | |
271 | BuildLQRefs(fTrainMomBin); | |
82719f10 | 272 | MonitorTraining(fTrainMomBin); |
aea8dd24 | 273 | } |
274 | // train all momenta | |
275 | else{ | |
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); | |
279 | continue; | |
280 | } | |
281 | TrainNetworks(iMomBin); | |
282 | BuildLQRefs(fTrainMomBin); | |
82719f10 | 283 | MonitorTraining(iMomBin); |
aea8dd24 | 284 | } |
285 | } | |
286 | ||
28efdace | 287 | return kTRUE; // testing protection |
288 | } | |
289 | ||
290 | ||
291 | //________________________________________________________________________ | |
292 | void AliTRDpidRefMaker::Terminate(Option_t *) | |
293 | { | |
294 | // Draw result to the screen | |
295 | // Called once at the end of the query | |
296 | ||
297 | fContainer = dynamic_cast<TObjArray*>(GetOutputData(0)); | |
298 | if (!fContainer) { | |
299 | Printf("ERROR: list not available"); | |
300 | return; | |
301 | } | |
302 | } | |
303 | ||
304 | ||
305 | //________________________________________________________________________ | |
aea8dd24 | 306 | void AliTRDpidRefMaker::GetV0info(AliTRDtrackV1 *TRDtrack, Float_t *v0pid) |
28efdace | 307 | { |
28efdace | 308 | // !!!! PREMILMINARY FUNCTION !!!! |
309 | // | |
310 | // this is the place for the V0 procedure | |
aea8dd24 | 311 | // as long as there is no one implemented, |
312 | // just the probabilities | |
313 | // of the TRDtrack are used! | |
28efdace | 314 | |
315 | TRDtrack -> SetReconstructor(fReconstructor); | |
316 | fReconstructor -> SetOption("nn"); | |
317 | TRDtrack -> CookPID(); | |
318 | for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){ | |
aea8dd24 | 319 | v0pid[iPart] = TRDtrack -> GetPID(iPart); |
320 | if(fDebugLevel>=4) Printf("PDG is (in V0info) %d %f", iPart, v0pid[iPart]); | |
28efdace | 321 | } |
322 | } | |
aea8dd24 | 323 | |
324 | ||
325 | //________________________________________________________________________ | |
326 | void AliTRDpidRefMaker::MakeTrainingLists() | |
327 | { | |
328 | // | |
329 | // build the training lists for the neural networks | |
330 | // | |
331 | ||
82719f10 | 332 | if (!fNN) { |
333 | LoadFiles("TRD.TaskPidRefMakerNN.root", "TRD.TaskPidRefMakerLQ.root"); | |
334 | } | |
335 | ||
336 | if (!fNN) { | |
337 | Printf("ERROR tree for training list not available"); | |
338 | return; | |
339 | } | |
340 | ||
aea8dd24 | 341 | if(fDebugLevel>=2) Printf("\n Making training lists! \n"); |
342 | ||
343 | Int_t nPart[AliPID::kSPECIES][AliTRDCalPID::kNMom]; | |
344 | memset(nPart, 0, AliPID::kSPECIES*AliTRDCalPID::kNMom*sizeof(Int_t)); | |
345 | ||
346 | // set needed branches | |
347 | fNN -> SetBranchAddress("fv0pid", &fv0pid); | |
348 | fNN -> SetBranchAddress("fMom", &fMom); | |
349 | fNN -> SetBranchAddress("fLayer", &fLayer); | |
350 | ||
82719f10 | 351 | AliTRDpidUtil *util = new AliTRDpidUtil(); |
352 | ||
aea8dd24 | 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); | |
356 | ||
357 | // use only events with goes through 6 layers TRD | |
358 | if(!fLayer == 0) | |
359 | continue; | |
360 | ||
361 | // set the 11 momentum bins | |
362 | Int_t iMomBin = -1; | |
82719f10 | 363 | iMomBin = util -> GetMomentumBin(fMom); |
aea8dd24 | 364 | |
365 | // check PID information and count particle types per momentum interval | |
82719f10 | 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]++; | |
aea8dd24 | 376 | } |
377 | ||
378 | if(fDebugLevel>=2){ | |
379 | Printf("Particle multiplicities:"); | |
380 | for(Int_t iMomBin = 0; iMomBin <AliTRDCalPID::kNMom; iMomBin++) | |
82719f10 | 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]); |
aea8dd24 | 382 | Printf("\n"); |
383 | } | |
384 | ||
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)); | |
389 | ||
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]; | |
397 | } | |
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]); | |
401 | } | |
402 | if(fDebugLevel>=2) Printf("\n"); | |
403 | ||
404 | ||
405 | // reset couters | |
406 | memset(nPart, 0, AliPID::kSPECIES*AliTRDCalPID::kNMom*sizeof(Int_t)); | |
407 | ||
408 | // start second loop to set the event lists | |
409 | for(Int_t iEv = 0; iEv < fNN -> GetEntries(); iEv++){ | |
410 | fNN -> GetEntry(iEv); | |
411 | ||
412 | // use only events with goes through 6 layers TRD | |
413 | if(!fLayer == 0) | |
414 | continue; | |
415 | ||
416 | // set the 11 momentum bins | |
417 | Int_t iMomBin = -1; | |
82719f10 | 418 | iMomBin = util -> GetMomentumBin(fMom); |
aea8dd24 | 419 | |
420 | // set electrons | |
82719f10 | 421 | if(fv0pid[AliPID::kElectron] == 1){ |
422 | if(nPart[AliPID::kElectron][iMomBin] < iTrain[iMomBin]){ | |
5d6dc395 | 423 | for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++) |
aea8dd24 | 424 | fTrain[iMomBin][ily] -> Enter(iEv + ily); |
82719f10 | 425 | nPart[AliPID::kElectron][iMomBin]++; |
aea8dd24 | 426 | } |
82719f10 | 427 | else if(nPart[AliPID::kElectron][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){ |
5d6dc395 | 428 | for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++) |
aea8dd24 | 429 | fTest[iMomBin][ily] -> Enter(iEv + ily); |
82719f10 | 430 | nPart[AliPID::kElectron][iMomBin]++; |
aea8dd24 | 431 | } |
432 | else | |
433 | continue; | |
434 | } | |
435 | // set muons | |
82719f10 | 436 | else if(fv0pid[AliPID::kMuon] == 1){ |
437 | if(nPart[AliPID::kMuon][iMomBin] < iTrain[iMomBin]){ | |
5d6dc395 | 438 | for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++) |
aea8dd24 | 439 | fTrain[iMomBin][ily] -> Enter(iEv + ily); |
82719f10 | 440 | nPart[AliPID::kMuon][iMomBin]++; |
aea8dd24 | 441 | } |
82719f10 | 442 | else if(nPart[AliPID::kMuon][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){ |
5d6dc395 | 443 | for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++) |
aea8dd24 | 444 | fTest[iMomBin][ily] -> Enter(iEv + ily); |
82719f10 | 445 | nPart[AliPID::kMuon][iMomBin]++; |
aea8dd24 | 446 | } |
447 | else | |
448 | continue; | |
449 | } | |
450 | // set pions | |
82719f10 | 451 | else if(fv0pid[AliPID::kPion] == 1){ |
452 | if(nPart[AliPID::kPion][iMomBin] < iTrain[iMomBin]){ | |
5d6dc395 | 453 | for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++) |
aea8dd24 | 454 | fTrain[iMomBin][ily] -> Enter(iEv + ily); |
82719f10 | 455 | nPart[AliPID::kPion][iMomBin]++; |
aea8dd24 | 456 | } |
82719f10 | 457 | else if(nPart[AliPID::kPion][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){ |
5d6dc395 | 458 | for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++) |
aea8dd24 | 459 | fTest[iMomBin][ily] -> Enter(iEv + ily); |
82719f10 | 460 | nPart[AliPID::kPion][iMomBin]++; |
aea8dd24 | 461 | } |
462 | else | |
463 | continue; | |
464 | } | |
465 | // set kaons | |
82719f10 | 466 | else if(fv0pid[AliPID::kKaon] == 1){ |
467 | if(nPart[AliPID::kKaon][iMomBin] < iTrain[iMomBin]){ | |
5d6dc395 | 468 | for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++) |
aea8dd24 | 469 | fTrain[iMomBin][ily] -> Enter(iEv + ily); |
82719f10 | 470 | nPart[AliPID::kKaon][iMomBin]++; |
aea8dd24 | 471 | } |
82719f10 | 472 | else if(nPart[AliPID::kKaon][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){ |
5d6dc395 | 473 | for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++) |
aea8dd24 | 474 | fTest[iMomBin][ily] -> Enter(iEv + ily); |
82719f10 | 475 | nPart[AliPID::kKaon][iMomBin]++; |
aea8dd24 | 476 | } |
477 | else | |
478 | continue; | |
479 | } | |
480 | // set protons | |
82719f10 | 481 | else if(fv0pid[AliPID::kProton] == 1){ |
482 | if(nPart[AliPID::kProton][iMomBin] < iTrain[iMomBin]){ | |
5d6dc395 | 483 | for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++) |
aea8dd24 | 484 | fTrain[iMomBin][ily] -> Enter(iEv + ily); |
82719f10 | 485 | nPart[AliPID::kProton][iMomBin]++; |
aea8dd24 | 486 | } |
82719f10 | 487 | else if(nPart[AliPID::kProton][iMomBin] < iTest[iMomBin]+iTrain[iMomBin]){ |
5d6dc395 | 488 | for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++) |
aea8dd24 | 489 | fTest[iMomBin][ily] -> Enter(iEv + ily); |
82719f10 | 490 | nPart[AliPID::kProton][iMomBin]++; |
aea8dd24 | 491 | } |
492 | else | |
493 | continue; | |
494 | } | |
495 | } | |
496 | ||
497 | if(fDebugLevel>=2){ | |
498 | Printf("Particle multiplicities in both lists:"); | |
499 | for(Int_t iMomBin = 0; iMomBin <AliTRDCalPID::kNMom; iMomBin++) | |
82719f10 | 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]); |
aea8dd24 | 501 | Printf("\n"); |
502 | } | |
82719f10 | 503 | |
504 | util -> Delete(); | |
aea8dd24 | 505 | } |
506 | ||
507 | ||
508 | //________________________________________________________________________ | |
509 | void AliTRDpidRefMaker::TrainNetworks(Int_t mombin) | |
510 | { | |
511 | // | |
512 | // train the neural networks | |
513 | // | |
514 | ||
82719f10 | 515 | |
516 | if (!fNN) { | |
517 | LoadFiles("TRD.TaskPidRefMakerNN.root", "TRD.TaskPidRefMakerLQ.root"); | |
518 | } | |
519 | ||
520 | if (!fNN) { | |
521 | Printf("ERROR tree for training list not available"); | |
522 | return; | |
523 | } | |
524 | ||
525 | TDatime datime; | |
526 | fDate = datime.GetDate(); | |
527 | ||
aea8dd24 | 528 | if(fDebugLevel>=2) Printf("Training momentum bin %d", mombin); |
529 | ||
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); | |
533 | ||
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)); | |
537 | ||
538 | // variable to check if network can load weights from previous training | |
5d6dc395 | 539 | Bool_t bFirstLoop[AliTRDgeometry::kNlayer]; |
540 | memset(bFirstLoop, kTRUE, AliTRDgeometry::kNlayer*sizeof(Bool_t)); | |
aea8dd24 | 541 | |
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 | |
5d6dc395 | 545 | for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){ |
aea8dd24 | 546 | // set the event lists |
547 | fNN -> SetEventList(fTrain[mombin][iChamb]); | |
548 | fNN -> SetEventList(fTest[mombin][iChamb]); | |
549 | ||
550 | if(fDebugLevel>=2) Printf("Trainingloop[%d] Chamber[%d]", iLoop, iChamb); | |
551 | ||
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; | |
560 | } | |
561 | else{ | |
82719f10 | 562 | if(fDebugLevel>=2) fNet[iChamb] -> Train(nEpochs,"text update=10, graph+"); |
aea8dd24 | 563 | else fNet[iChamb] -> Train(nEpochs,"+"); |
564 | } | |
565 | ||
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 | |
570 | } | |
571 | ||
572 | ||
573 | //________________________________________________________________________ | |
574 | void AliTRDpidRefMaker::BuildLQRefs(Int_t mombin) | |
575 | { | |
576 | // | |
577 | // build the 2-dim LQ reference histograms | |
578 | // | |
579 | ||
580 | if(fDebugLevel>=2) Printf("Building LQRefs for momentum bin %d", mombin); | |
581 | } | |
582 | ||
583 | ||
584 | //________________________________________________________________________ | |
82719f10 | 585 | void AliTRDpidRefMaker::MonitorTraining(Int_t mombin) |
aea8dd24 | 586 | { |
587 | // | |
588 | // train the neural networks | |
589 | // | |
590 | ||
82719f10 | 591 | if(!fContainer){ |
592 | LoadContainer("TRD.TaskPidRefMaker.root"); | |
593 | } | |
594 | if(!fContainer){ | |
595 | Printf("ERROR container not available"); | |
596 | return; | |
597 | } | |
598 | ||
599 | if (!fNN) { | |
600 | LoadFiles("TRD.TaskPidRefMakerNN.root", "TRD.TaskPidRefMakerLQ.root"); | |
601 | } | |
602 | if (!fNN) { | |
603 | Printf("ERROR tree for training list not available"); | |
604 | return; | |
605 | } | |
606 | ||
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]); | |
612 | } | |
613 | ||
614 | // implement variables for likelihoods | |
615 | Float_t Like[AliPID::kSPECIES][AliTRDgeometry::kNlayer]; | |
a391a274 | 616 | memset(Like, 0, AliPID::kSPECIES*AliTRDgeometry::kNlayer*sizeof(Float_t)); |
82719f10 | 617 | Float_t LikeAll[AliPID::kSPECIES], TotProb; |
618 | ||
619 | Double_t PionEffiTrain[kMoniTrain], PionEffiErrTrain[kMoniTrain]; | |
620 | Double_t PionEffiTest[kMoniTrain], PionEffiErrTest[kMoniTrain]; | |
a391a274 | 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)); | |
82719f10 | 625 | |
626 | // init histos | |
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); | |
631 | ||
632 | TGraphErrors *gEffisTrain=0x0, *gEffisTest=0x0; | |
633 | gEffisTrain = (TGraphErrors*)fContainer->At(kGraphTrain); | |
634 | gEffisTest = (TGraphErrors*)fContainer->At(kGraphTest); | |
635 | ||
636 | AliTRDpidUtil *util = new AliTRDpidUtil(); | |
637 | ||
638 | // monitor training progress | |
639 | for(Int_t iLoop = 0; iLoop < kMoniTrain; iLoop++){ | |
640 | ||
641 | // load weights | |
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)); | |
644 | } | |
645 | ||
646 | // // debug loop | |
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)); | |
654 | ||
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]++; | |
665 | // } | |
666 | // Printf("Numbers [%d] [%d] [%d] [%d] [%d]", iType[AliPID::kElectron], iType[AliPID::kMuon], iType[AliPID::kPion], iType[AliPID::kKaon], iType[AliPID::kProton]); | |
667 | // } | |
668 | ||
669 | ||
670 | // event loop training list | |
671 | for(Int_t iEvent = 0; iEvent < fTrain[mombin][0] -> GetN(); iEvent++ ){ | |
672 | ||
673 | // reset particle probabilities | |
674 | for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){ | |
675 | LikeAll[iPart] = 1./AliPID::kSPECIES; | |
676 | } | |
677 | TotProb = 0.; | |
678 | ||
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; | |
682 | ||
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]; | |
688 | } | |
689 | } | |
690 | ||
691 | // get total probability and normalize it | |
692 | for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){ | |
693 | TotProb += LikeAll[iPart]; | |
694 | } | |
695 | for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){ | |
696 | LikeAll[iPart] /= TotProb; | |
697 | } | |
698 | ||
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]); | |
704 | } // end event loop | |
705 | ||
706 | ||
707 | // calculate the pion efficiency and fill the graph | |
708 | util -> CalculatePionEffi(hElecs, hPions); | |
709 | PionEffiTrain[iLoop] = util -> GetPionEfficiency(); | |
710 | PionEffiErrTrain[iLoop] = util -> GetError(); | |
711 | ||
712 | gEffisTrain -> SetPoint(iLoop, iLoop+1, PionEffiTrain[iLoop]); | |
713 | gEffisTrain -> SetPointError(iLoop, 0, PionEffiErrTrain[iLoop]); | |
714 | hElecs -> Reset(); | |
715 | hPions -> Reset(); | |
716 | if(fDebugLevel>=2) Printf("TrainingLoop[%d] PionEfficiency[%f +/- %f]", iLoop, PionEffiTrain[iLoop], PionEffiErrTrain[iLoop]); | |
717 | // end training loop | |
718 | ||
719 | ||
720 | ||
721 | // event loop test list | |
722 | for(Int_t iEvent = 0; iEvent < fTest[mombin][0] -> GetN(); iEvent++ ){ | |
723 | ||
724 | // reset particle probabilities | |
725 | for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){ | |
726 | LikeAll[iPart] = 1./AliTRDgeometry::kNlayer; | |
727 | } | |
728 | TotProb = 0.; | |
729 | ||
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; | |
733 | ||
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]; | |
739 | } | |
740 | } | |
741 | ||
742 | // get total probability and normalize it | |
743 | for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){ | |
744 | TotProb += LikeAll[iPart]; | |
745 | } | |
746 | for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){ | |
747 | LikeAll[iPart] /= TotProb; | |
748 | } | |
749 | ||
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]); | |
755 | } // end event loop | |
756 | ||
757 | // calculate the pion efficiency and fill the graph | |
758 | util -> CalculatePionEffi(hElecs, hPions); | |
759 | PionEffiTest[iLoop] = util -> GetPionEfficiency(); | |
760 | PionEffiErrTest[iLoop] = util -> GetError(); | |
761 | ||
762 | gEffisTest -> SetPoint(iLoop, iLoop+1, PionEffiTest[iLoop]); | |
763 | gEffisTest -> SetPointError(iLoop, 0, PionEffiErrTest[iLoop]); | |
764 | hElecs -> Reset(); | |
765 | hPions -> Reset(); | |
766 | if(fDebugLevel>=2) Printf("TestLoop[%d] PionEfficiency[%f +/- %f] \n", iLoop, PionEffiTest[iLoop], PionEffiErrTest[iLoop]); | |
767 | ||
768 | } // end training loop | |
769 | ||
770 | util -> Delete(); | |
771 | ||
772 | gEffisTest -> Draw("PAL"); | |
773 | gEffisTrain -> Draw("PL"); | |
774 | ||
aea8dd24 | 775 | } |
82719f10 | 776 | |
777 | ||
778 | //________________________________________________________________________ | |
779 | void AliTRDpidRefMaker::LoadFiles(const Char_t *InFileNN, const Char_t *InFileLQ) | |
780 | { | |
781 | // | |
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 | |
786 | // | |
787 | ||
788 | TFile *fInFileNN; | |
789 | fInFileNN = new TFile(InFileNN, "READ"); | |
790 | fNN = (TTree*)fInFileNN -> Get("NN"); | |
791 | ||
792 | TFile *fInFileLQ; | |
793 | fInFileLQ = new TFile(InFileLQ, "READ"); | |
794 | fLQ = (TTree*)fInFileLQ -> Get("LQ"); | |
795 | ||
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)); | |
800 | } | |
801 | } | |
802 | } | |
803 | ||
804 | ||
805 | //________________________________________________________________________ | |
806 | void AliTRDpidRefMaker::LoadContainer(const Char_t *InFileCont) | |
807 | { | |
808 | ||
809 | // | |
810 | // Loads the container if no container is there. | |
811 | // Useable for training outside of the makeResults.C macro | |
812 | // | |
813 | ||
814 | TFile *fInFileCont; | |
815 | fInFileCont = new TFile(InFileCont, "READ"); | |
816 | fContainer = (TObjArray*)fInFileCont -> Get("PidRefMaker"); | |
817 | ||
818 | } | |
819 | ||
820 | ||
821 | // //________________________________________________________________________ | |
822 | // void AliTRDpidRefMaker::CreateGraphs() | |
823 | // { | |
824 | // // Create histograms | |
825 | // // Called once | |
826 | ||
827 | // OpenFile(0, "RECREATE"); | |
828 | // fContainer = new TObjArray(); | |
829 | // fContainer->AddAt(new TH1F("hPDG","hPDG",AliPID::kSPECIES,-0.5,5.5),0); | |
830 | ||
831 | // TGraphErrors *gEffisTrain = new TGraphErrors(kMoniTrain); | |
832 | // gEffisTrain -> SetLineColor(4); | |
833 | // gEffisTrain -> SetMarkerColor(4); | |
834 | // gEffisTrain -> SetMarkerStyle(29); | |
835 | // gEffisTrain -> SetMarkerSize(2); | |
836 | ||
837 | // TGraphErrors *gEffisTest = new TGraphErrors(kMoniTrain); | |
838 | // gEffisTest -> SetLineColor(2); | |
839 | // gEffisTest -> SetMarkerColor(2); | |
840 | // gEffisTest -> SetMarkerSize(2); | |
841 | ||
842 | // fContainer -> AddAt(gEffisTrain,kGraphTrain); | |
843 | // fContainer -> AddAt(gEffisTest,kGraphTest); | |
844 | // } | |
845 |