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

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