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

Last change on this file since 7273 was 7178, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 9.6 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 <TH1.h>
37#include <TH2.h>
38#include <TStyle.h>
39#include <TCanvas.h>
40#include <TMath.h>
41#include <TVector.h>
42
43#include "MHMatrix.h"
44
45#include "MLog.h"
46#include "MLogManip.h"
47
48#include "MParList.h"
49#include "MTaskList.h"
50#include "MEvtLoop.h"
51
52#include "MRanTree.h"
53#include "MRanForest.h"
54#include "MRanForestGrow.h"
55
56#include "MData.h"
57#include "MParameters.h"
58
59ClassImp(MRFEnergyEst);
60
61using namespace std;
62
63static const TString gsDefName = "MRFEnergyEst";
64static const TString gsDefTitle = "RF for energy estimation";
65
66// --------------------------------------------------------------------------
67//
68// Default constructor. Set the name and title of this object
69//
70MRFEnergyEst::MRFEnergyEst(const char *name, const char *title)
71 : fDebug(kFALSE), fData(0), fEnergyEst(0),
72 fNumTrees(-1), fNumTry(-1), fNdSize(-1),
73 fTestMatrix(0), fEstimationMode(kMean)
74{
75 fName = name ? name : gsDefName.Data();
76 fTitle = title ? title : gsDefTitle.Data();
77}
78
79MRFEnergyEst::~MRFEnergyEst()
80{
81 fEForests.Delete();
82}
83
84// --------------------------------------------------------------------------
85//
86// Train the RF with the goven MHMatrix. The last column must contain the
87// True energy.
88//
89Int_t MRFEnergyEst::Train(const MHMatrix &matrixtrain, const TArrayD &grid)
90{
91 gLog.Separator("MRFEnergyEst - Train");
92
93 if (!matrixtrain.GetColumns())
94 {
95 *fLog << err << "ERROR - MHMatrix does not contain rules... abort." << endl;
96 return kFALSE;
97 }
98
99 const Int_t ncols = matrixtrain.GetM().GetNcols();
100 const Int_t nrows = matrixtrain.GetM().GetNrows();
101 if (ncols<=0 || nrows <=0)
102 {
103 *fLog << err << "ERROR - No. of columns or no. of rows of matrixtrain equal 0 ... abort." << endl;
104 return kFALSE;
105 }
106
107 const Int_t nbins = grid.GetSize()-1;
108 if (nbins<=0)
109 {
110 *fLog << err << "ERROR - Energy grid not vaild... abort." << endl;
111 return kFALSE;
112 }
113
114 // rules (= combination of image par) to be used for energy estimation
115 TFile fileRF(fFileName, "recreate");
116 if (!fileRF.IsOpen())
117 {
118 *fLog << err << "ERROR - File to store RFs could not be opened... abort." << endl;
119 return kFALSE;
120 }
121
122 const Int_t nobs = 3; // Number of obsolete columns
123
124 MDataArray usedrules;
125 for(Int_t i=0; i<ncols-nobs; i++) // -3 is important!!!
126 usedrules.AddEntry((*matrixtrain.GetColumns())[i].GetRule());
127
128 // training of RF for each energy bin
129 for (Int_t ie=0; ie<nbins; ie++)
130 {
131 TMatrix mat1(nrows, ncols);
132 TMatrix mat0(nrows, ncols);
133
134 // prepare matrix for current energy bin
135 Int_t irow1=0;
136 Int_t irow0=0;
137
138 const TMatrix &m = matrixtrain.GetM();
139 for (Int_t j=0; j<nrows; j++)
140 {
141 const Double_t energy = m(j,ncols-1);
142
143 if (energy>grid[ie] && energy<=grid[ie+1])
144 TMatrixFRow(mat1, irow1++) = TMatrixFRow_const(m,j);
145 else
146 TMatrixFRow(mat0, irow0++) = TMatrixFRow_const(m,j);
147 }
148
149 const Bool_t invalid = irow1==0 || irow0==0;
150
151 if (invalid)
152 *fLog << warn << "WARNING - Skipping";
153 else
154 *fLog << inf << "Training RF for";
155
156 *fLog << " energy bin " << ie << " (" << grid[ie] << ", " << grid[ie+1] << ") " << irow0 << " " << irow1 << endl;
157
158 if (invalid)
159 continue;
160
161 if (fDebug)
162 gLog.SetNullOutput(kTRUE);
163
164 // last column must be removed (true energy col.)
165 mat1.ResizeTo(irow1, ncols-nobs);
166 mat0.ResizeTo(irow0, ncols-nobs);
167
168 // create MHMatrix as input for RF
169 MHMatrix matrix1(mat1, "MatrixHadrons");
170 MHMatrix matrix0(mat0, "MatrixGammas");
171
172 // training of RF
173 MTaskList tlist;
174 MParList plist;
175 plist.AddToList(&tlist);
176 plist.AddToList(&matrix0);
177 plist.AddToList(&matrix1);
178
179 MRanForestGrow rfgrow;
180 rfgrow.SetNumTrees(fNumTrees); // number of trees
181 rfgrow.SetNumTry(fNumTry); // number of trials in random split selection
182 rfgrow.SetNdSize(fNdSize); // limit for nodesize
183
184 tlist.AddToList(&rfgrow);
185
186 MEvtLoop evtloop;
187 evtloop.SetDisplay(fDisplay);
188 evtloop.SetParList(&plist);
189
190 if (!evtloop.Eventloop())
191 return kFALSE;
192
193 if (fDebug)
194 gLog.SetNullOutput(kFALSE);
195
196 // Calculate bin center
197 const Double_t E = (TMath::Log10(grid[ie])+TMath::Log10(grid[ie+1]))/2;
198
199 // save whole forest
200 MRanForest *forest=(MRanForest*)plist.FindObject("MRanForest");
201 forest->SetUserVal(E);
202 forest->Write(Form("%.10f", E));
203 }
204
205 // save rules
206 usedrules.Write("rules");
207
208 return kTRUE;
209}
210
211Int_t MRFEnergyEst::ReadForests(MParList &plist)
212{
213 TFile fileRF(fFileName, "read");
214 if (!fileRF.IsOpen())
215 {
216 *fLog << err << dbginf << "File containing RFs could not be opened... aborting." << endl;
217 return kFALSE;
218 }
219
220 fEForests.Delete();
221
222 TIter Next(fileRF.GetListOfKeys());
223 TObject *o=0;
224 while ((o=Next()))
225 {
226 MRanForest *forest;
227 fileRF.GetObject(o->GetName(), forest);
228 if (!forest)
229 continue;
230
231 forest->SetUserVal(atof(o->GetName()));
232
233 fEForests.Add(forest);
234 }
235
236 if (fData->Read("rules")<=0)
237 {
238 *fLog << err << "ERROR - Reading 'rules' from file " << fFileName << endl;
239 return kFALSE;
240 }
241
242 return kTRUE;
243}
244
245Int_t MRFEnergyEst::PreProcess(MParList *plist)
246{
247 fEnergyEst = (MParameterD*)plist->FindCreateObj("MParameterD", "MEnergyEst");
248 if (!fEnergyEst)
249 return kFALSE;
250
251 fData = (MDataArray*)plist->FindCreateObj("MDataArray");
252 if (!fData)
253 return kFALSE;
254
255 if (!ReadForests(*plist))
256 {
257 *fLog << err << "Reading RFs failed... aborting." << endl;
258 return kFALSE;
259 }
260
261 *fLog << inf << "RF read from " << fFileName << endl;
262
263 if (fTestMatrix)
264 return kTRUE;
265
266 fData->Print();
267
268 if (!fData->PreProcess(plist))
269 {
270 *fLog << err << "PreProcessing of the MDataArray failed... aborting." << endl;
271 return kFALSE;
272 }
273
274 return kTRUE;
275}
276
277// --------------------------------------------------------------------------
278//
279//
280#include <TGraph.h>
281#include <TF1.h>
282Int_t MRFEnergyEst::Process()
283{
284 static TF1 f1("f1", "gaus");
285
286 TVector event;
287 if (fTestMatrix)
288 *fTestMatrix >> event;
289 else
290 *fData >> event;
291
292 Double_t sume = 0;
293 Double_t sumh = 0;
294 Double_t maxh = 0;
295 Double_t maxe = 0;
296
297 Double_t max = -1e10;
298 Double_t min = 1e10;
299
300 //TH1C h("", "", fEForests.GetSize(), 0, 1);
301
302 TIter Next(&fEForests);
303 MRanForest *rf = 0;
304
305 TGraph g;
306 while ((rf=(MRanForest*)Next()))
307 {
308 const Double_t h = rf->CalcHadroness(event);
309 const Double_t e = rf->GetUserVal();
310 g.SetPoint(g.GetN(), e, h);
311 sume += e*h;
312 sumh += h;
313 if (h>maxh)
314 {
315 maxh = h;
316 maxe = e;
317 }
318 if (e>max)
319 max = e;
320 if (e<min)
321 min = e;
322 }
323
324 switch (fEstimationMode)
325 {
326 case kMean:
327 fEnergyEst->SetVal(pow(10, sume/sumh));
328 break;
329 case kMaximum:
330 fEnergyEst->SetVal(pow(10, maxe));
331 break;
332 case kFit:
333 f1.SetParameter(0, maxh);
334 f1.SetParameter(1, maxe);
335 f1.SetParameter(2, 0.125);
336 g.Fit(&f1, "Q0N");
337 fEnergyEst->SetVal(pow(10, f1.GetParameter(1)));
338 break;
339
340 }
341 fEnergyEst->SetReadyToSave();
342
343 return kTRUE;
344}
345
346// --------------------------------------------------------------------------
347//
348//
349Int_t MRFEnergyEst::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
350{
351 Bool_t rc = kFALSE;
352 if (IsEnvDefined(env, prefix, "FileName", print))
353 {
354 rc = kTRUE;
355 SetFileName(GetEnvValue(env, prefix, "FileName", fFileName));
356 }
357 if (IsEnvDefined(env, prefix, "Debug", print))
358 {
359 rc = kTRUE;
360 SetDebug(GetEnvValue(env, prefix, "Debug", fDebug));
361 }
362 if (IsEnvDefined(env, prefix, "EstimationMode", print))
363 {
364 TString txt = GetEnvValue(env, prefix, "EstimationMode", "");
365 txt = txt.Strip(TString::kBoth);
366 txt.ToLower();
367 if (txt==(TString)"mean")
368 fEstimationMode = kMean;
369 if (txt==(TString)"maximum")
370 fEstimationMode = kMaximum;
371 if (txt==(TString)"fit")
372 fEstimationMode = kFit;
373 rc = kTRUE;
374 }
375 return rc;
376}
Note: See TracBrowser for help on using the repository browser.