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

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