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

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