source: trunk/MagicSoft/Mars/mranforest/MRFEnergyEst.cc@ 6724

Last change on this file since 6724 was 6540, checked in by hengsteb, 20 years ago
*** empty log message ***
File size: 11.7 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Hengstebeck 2/2005 <mailto:hengsteb@physik.hu-berlin.de>
19!
20! Copyright: MAGIC Software Development, 2000-2005
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MRFEnergyEst
28//
29//
30////////////////////////////////////////////////////////////////////////////
31#include "MRFEnergyEst.h"
32
33#include <TFile.h>
34#include <TList.h>
35
36#include <TH1F.h>
37#include <TH2F.h>
38#include <TStyle.h>
39#include <TCanvas.h>
40
41#include "MHMatrix.h"
42
43#include "MLog.h"
44#include "MLogManip.h"
45
46#include "MParList.h"
47#include "MTaskList.h"
48#include "MEvtLoop.h"
49
50#include "MRanTree.h"
51#include "MRanForest.h"
52#include "MRanForestGrow.h"
53
54#include "MData.h"
55
56#include "TFile.h"
57#include "TList.h"
58
59#include "TH1F.h"
60#include "TH2F.h"
61#include "TStyle.h"
62#include "TCanvas.h"
63#include "TVector.h"
64
65#include "MEnergyEst.h"
66
67
68ClassImp(MRFEnergyEst);
69
70using namespace std;
71
72static const TString gsDefName = "MRFEnergyEst";
73static const TString gsDefTitle = "RF for energy estimation";
74
75// --------------------------------------------------------------------------
76//
77//
78MRFEnergyEst::MRFEnergyEst(const char *name, const char *title):fNumTrees(-1)
79{
80 //
81 // set the name and title of this object
82 //
83 fName = name ? name : gsDefName.Data();
84 fTitle = title ? title : gsDefTitle.Data();
85}
86
87// --------------------------------------------------------------------------
88//
89// Delete the data chains
90//
91MRFEnergyEst::~MRFEnergyEst()
92{
93 // delete fData;
94}
95
96// --------------------------------------------------------------------------
97Int_t MRFEnergyEst::Train()
98{
99 if(!fMatrixTrain)
100 {
101 *fLog << err << dbginf << "fMatrixTrain not found... aborting." << endl;
102 return kFALSE;
103 }
104
105 if(!fMatrixTrain->GetColumns())
106 {
107 *fLog << err << dbginf << "fMatrixTrain does not contain rules... aborting." << endl;
108 return kFALSE;
109 }
110
111 const Int_t ncols = (fMatrixTrain->GetM()).GetNcols();
112 const Int_t nrows = (fMatrixTrain->GetM()).GetNrows();
113
114 cout<<"ncols="<<ncols<<" nrows="<<nrows<<endl<<flush;
115
116 if(ncols<=0 || nrows <=0)
117 {
118 *fLog << err << dbginf << "No. of columns or no. of rows of fMatrixTrain equal 0 ... aborting." << endl;
119 return kFALSE;
120 }
121
122 // rules (= combination of image par) to be used for energy estimation
123 MDataArray used_rules;
124 TString energy_rule;
125 for(Int_t i=0;i<ncols;i++)
126 {
127 MDataArray *rules=fMatrixTrain->GetColumns();
128 MData &data=(*rules)[i];
129
130 if(i<ncols-1)
131 used_rules.AddEntry(data.GetRule());
132 else
133 energy_rule=data.GetRule();
134 }
135
136 if(energy_rule.CompareTo("MMcEvt.fEnergy"))
137 {
138 *fLog << err << dbginf << "Can not accept "<<energy_rule<<" as true energy ... aborting." << endl;
139 return kFALSE;
140 }
141
142 TFile fileRF(fRFfileName,"recreate");
143 if(!fileRF.IsOpen())
144 {
145 *fLog << err << dbginf << "File to store RFs could not be opened... aborting." << endl;
146 return kFALSE;
147 }
148
149 TMatrix *mptr=(TMatrix*)&(fMatrixTrain->GetM());
150 const Int_t nbins = fEnergyGrid.GetSize()-1;
151 if(nbins<=0)
152 {
153 *fLog << err << dbginf << "Energy grid not vaild... aborting." << endl;
154 return kFALSE;
155 }
156
157 // training of RF for each energy bin
158 Int_t numbins=0;
159 for(Int_t ie=0;ie<nbins;ie++)
160 {
161 TMatrix mat1(nrows,ncols);
162 TMatrix mat0(nrows,ncols);
163
164 // prepare matrix for current energy bin
165 Int_t irow1=0;
166 Int_t irow0=0;
167
168 for(Int_t j=0;j<nrows;j++)
169 {
170 Double_t energy=(*mptr)(j,ncols-1);
171 if(energy>pow(10.,fEnergyGrid[ie]) && energy<=pow(10.,fEnergyGrid[ie+1]))
172 {
173 TMatrixRow(mat1, irow1) = TMatrixRow(*mptr,j);
174 irow1++;
175 }else{
176 TMatrixRow(mat0, irow0) = TMatrixRow(*mptr,j);
177 irow0++;
178 }
179 }
180 if(irow1*irow0==0)continue;
181
182 *fLog << inf << dbginf << "Training RF for energy bin "<<ie<< endl;
183
184 // last column must be removed (true energy col.)
185 mat1.ResizeTo(irow1,ncols-1);
186 mat0.ResizeTo(irow0,ncols-1);
187
188 // create MHMatrix as input for RF
189 MHMatrix matrix1(mat1,"MatrixHadrons");
190 MHMatrix matrix0(mat0,"MatrixGammas");
191
192 MDataArray *rules1=matrix1.GetColumns();
193 MDataArray *rules0=matrix0.GetColumns();
194 // rules of new matrices be re-set
195 if(rules1)delete rules1; rules1=new MDataArray();
196 if(rules0)delete rules0; rules0=new MDataArray();
197
198 for(Int_t i=0;i<ncols-1;i++)
199 {
200 //MDataArray *rules=fMatrixTrain->GetColumns();
201 //MData &data=(*rules)[i];
202 MData &data=used_rules[i];
203 rules1->AddEntry(data.GetRule());
204 rules0->AddEntry(data.GetRule());
205 }
206
207 // training of RF
208 MParList plist;
209 MTaskList tlist;
210 plist.AddToList(&tlist);
211 plist.AddToList(&matrix0);
212 plist.AddToList(&matrix1);
213
214 MRanForestGrow rfgrow;
215 rfgrow.SetNumTrees(fNumTrees); // number of trees
216 rfgrow.SetNumTry(fNumTry); // number of trials in random split selection
217 rfgrow.SetNdSize(fNdSize); // limit for nodesize
218
219 tlist.AddToList(&rfgrow);
220
221 MEvtLoop evtloop;
222 evtloop.SetParList(&plist);
223
224 //gLog.SetNullOutput(kTRUE);
225
226 if (!evtloop.Eventloop()) return kFALSE;
227
228 //gLog.SetNullOutput(kFALSE);
229
230 // save whole forest
231 MRanForest *forest=(MRanForest*)plist.FindObject("MRanForest");
232 forest->SetTitle(Form("%f",0.5*(fEnergyGrid[ie]+fEnergyGrid[ie+1])));
233 forest->Write(Form("%d",numbins));
234 numbins++;
235 }
236
237 // save rules
238 used_rules.Write("rules");
239
240 fileRF.Close();
241
242 return kTRUE;
243}
244
245Int_t MRFEnergyEst::Test()
246{
247 if(!fMatrixTest)
248 {
249 *fLog << err << dbginf << "fMatrixTest not found... aborting." << endl;
250 return kFALSE;
251 }
252
253 const Int_t ncols = (fMatrixTest->GetM()).GetNcols();
254 const Int_t nrows = (fMatrixTest->GetM()).GetNrows();
255
256 if(ncols<=0 || nrows <=0)
257 {
258 *fLog << err << dbginf << "No. of columns or no. of rows of fMatrixTrain equal 0 ... aborting." << endl;
259 return kFALSE;
260 }
261
262 TMatrix *mptr=(TMatrix*)&(fMatrixTest->GetM());
263
264 if(!ReadForests())
265 {
266 *fLog << err << dbginf << "Reading RFs failed... aborting." << endl;
267 return kFALSE;
268 }
269
270 const Int_t nbins=fEForests.GetSize();
271
272 Double_t e_low = 100;
273 Double_t e_up = 0;
274
275 for(Int_t j=0;j<nbins;j++)
276 {
277 e_low = min(atof((fEForests[j])->GetTitle()),e_low);
278 e_up = max(atof((fEForests[j])->GetTitle()),e_up);
279 }
280
281 TH1F hres("hres","",1000,-10,10);
282 TH2F henergy("henergy","",100,e_low,e_up,100,e_low,e_up);
283
284 for(Int_t i=0;i<nrows;i++)
285 {
286 Double_t e_true = (*mptr)(i,ncols-1);
287 Double_t e_est = 0;
288 Double_t hmax = 0;
289 Double_t hsum = 0;
290
291 for(Int_t j=0;j<nbins;j++)
292 {
293 const TVector v=TMatrixRow(*mptr,i);
294 Double_t h = ((MRanForest*) (fEForests[j]))->CalcHadroness(v);
295 Double_t e = atof((fEForests[j])->GetTitle());
296 /*if(h>=hmax)
297 {
298 hmax=h;
299 e_est=pow(10.,e);
300 }*/
301 hsum+=h;
302 e_est+=h*e;
303 }
304 e_est/=hsum;
305 e_est=pow(10.,e_est);
306
307 if(e_true>80.) hres.Fill((e_est-e_true)/e_true);
308 henergy.Fill(log10(e_true),log10(e_est));
309 }
310
311 if(TestBit(kEnableGraphicalOutput))
312 {
313 if(gStyle) gStyle->SetOptFit(1);
314 // show results
315 TCanvas *c=MH::MakeDefCanvas("c","",300,900);
316
317 c->Divide(1,3);
318 c->cd(1);
319 henergy.SetTitle("Estimated vs true energy");
320 henergy.GetXaxis()->SetTitle("log10(E_{true}[GeV])");
321 henergy.GetYaxis()->SetTitle("log10(E_{est}[GeV])");
322 henergy.DrawCopy();
323
324 c->cd(2);
325
326 TH1F *hptr=(TH1F*)henergy.ProfileX("_px",1,100,"S");
327 hptr->SetTitle("Estimated vs true energy - profile");
328 hptr->GetXaxis()->SetTitle("log10(E_{true}[GeV])");
329 hptr->GetYaxis()->SetTitle("log10(E_{est}[GeV])");
330 hptr->DrawCopy();
331
332 c->cd(3);
333 hres.Fit("gaus");
334 hres.SetTitle("Energy resolution for E_{true}>80Gev");
335 hres.GetXaxis()->SetTitle("(E_{estimated}-E_{true})/E_{true})");
336 hres.GetYaxis()->SetTitle("counts");
337 hres.DrawCopy();
338
339
340 c->GetPad(1)->SetGrid();
341 c->GetPad(2)->SetGrid();
342 c->GetPad(3)->SetGrid();
343
344 }
345
346 return kTRUE;
347}
348
349Int_t MRFEnergyEst::ReadForests(MParList *plist)
350{
351 TFile fileRF(fRFfileName,"read");
352 if(!fileRF.IsOpen())
353 {
354 *fLog << err << dbginf << "File containing RFs could not be opened... aborting." << endl;
355 return kFALSE;
356 }
357
358 TList *list=(TList*)fileRF.GetListOfKeys();
359 const Int_t n=list->GetSize()-1;// subtract 1 since 1 key belongs to MDataArray
360
361 fEForests.Expand(n);
362
363 MRanForest forest;
364 for(Int_t i=0;i<n;i++)
365 {
366 forest.Read(Form("%d",i));
367 MRanForest *curforest=(MRanForest*)forest.Clone();
368 const char *energy=list->FindObject(Form("%d",i))->GetTitle();
369 const char *bin =list->FindObject(Form("%d",i))->GetName();
370
371 curforest->SetTitle(energy);
372 curforest->SetName(bin);
373
374 fEForests[i]=curforest;
375 }
376
377 if(plist)
378 {
379 fData=(MDataArray*)plist->FindCreateObj("MDataArray");
380 fData->Read("rules");
381 }
382 fileRF.Close();
383
384 return kTRUE;
385}
386
387Int_t MRFEnergyEst::PreProcess(MParList *plist)
388{
389 fEnergyEst = (MEnergyEst*)plist->FindCreateObj("MEnergyEst");
390 if (!fEnergyEst)
391 {
392 *fLog << err << dbginf << "MEnergyEst not found... aborting." << endl;
393 return kFALSE;
394 }
395
396 if(!ReadForests(plist))
397 {
398 *fLog << err << dbginf << "Reading RFs failed... aborting." << endl;
399 return kFALSE;
400 }
401
402 if (!fData)
403 {
404 *fLog << err << dbginf << "MDataArray not found... aborting." << endl;
405 return kFALSE;
406 }
407
408 if (!fData->PreProcess(plist))
409 {
410 *fLog << err << dbginf << "PreProcessing of the MDataArray failed for the columns failed... aborting." << endl;
411 return kFALSE;
412 }
413
414 return kTRUE;
415}
416
417// --------------------------------------------------------------------------
418//
419//
420Int_t MRFEnergyEst::Process()
421{
422 TVector event;
423 *fData >> event;
424
425 Double_t e_est = 0;
426 Double_t hmax = 0;
427 Double_t hsum = 0;
428
429 for(Int_t j=0;j<fEForests.GetSize();j++)
430 {
431 Double_t h = ((MRanForest*) (fEForests[j]))->CalcHadroness(event);
432 Double_t e = atof((fEForests[j])->GetTitle());
433 /*if(h>=hmax)
434 {
435 hmax=h;
436 e_est=pow(10.,e);
437 }*/
438 hsum+=h;
439 e_est+=h*e;
440 }
441 e_est/=hsum;
442 e_est=pow(10.,e_est);
443
444 fEnergyEst->SetEnergy(e_est);
445 fEnergyEst->SetReadyToSave();
446
447 return kTRUE;
448}
449
450Int_t MRFEnergyEst::PostProcess()
451{
452
453 return kTRUE;
454}
Note: See TracBrowser for help on using the repository browser.