source: trunk/MagicSoft/Mars/macros/CT1EgyEst.C@ 4608

Last change on this file since 4608 was 2357, checked in by wittek, 21 years ago
*** empty log message ***
File size: 9.5 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 Bretz, 09/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
19! Abelardo Moralejo, 03/2003 <mailto:moralejo@pd.infn.it>
20! Wolfgang Wittek, 04/2003 <mailto:wittek@mppmu.mpg.de>
21!
22! Copyright: MAGIC Software Development, 2000-2003
23!
24!
25\* ======================================================================== */
26
27#include "MParList.h"
28#include "MTaskList.h"
29#include "MGeomCamCT1.h"
30#include "MFEventSelector.h"
31#include "MReadTree.h"
32#include "MFCT1SelFinal.h"
33#include "MHMatrix.h"
34#include "MEnergyEstParam.h"
35//#include "MEnergyEstParamDanielMkn421.h"
36#include "MMatrixLoop.h"
37#include "MChisqEval.h"
38
39/*
40#include "MLog.h"
41#include "MLogManip.h"
42*/
43
44void CT1EgyEst()
45{
46 // TString inPath = "~magican/ct1test/wittek/marsoutput/optionA/";
47 TString inPath = "./";
48 TString fileNameIn = "MC1.root";
49
50 // TString outPath = "~magican/ct1test/wittek/marsoutput/optionA/";
51 TString outPath = "./";
52 TString energyParName = "energyest.root";
53
54 TString hilName = "MHillas";
55 TString hilSrcName = "MHillasSrc";
56
57 // TString hadronnessName = "MHadronness";
58 TString hadronnessName = "HadNN";
59
60 // Int_t howMany = 10000;
61 Int_t howMany = 2000;
62
63 Float_t maxhadronness = 0.3;
64 Float_t maxalpha = 20.0;
65 Float_t maxdist = 1.;
66
67 MBinning BinningE("BinningE");
68 MBinning BinningTheta("BinningTheta");
69 BinningE.SetEdgesLog(50, 500, 50000.);
70
71 const Double_t yedge[7] = {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50.};
72 const TArrayD yed(7,yedge);
73 BinningTheta.SetEdges(yed);
74
75 CT1EEst(inPath, fileNameIn,
76 outPath, energyParName,
77 hilName, hilSrcName, hadronnessName,
78 howMany, maxhadronness, maxalpha, maxdist,
79 &BinningE, &BinningTheta);
80}
81
82//------------------------------------------------------------------------
83//
84// Arguments of CT1EEst :
85// ------------------------
86//
87// inPath, fileNameIn path and name of input file (MC gamma events)
88// outPath, energyParName path and name of file containing the parameters
89// of the energy estimator
90// hilName, hilSrcName names of "MHillas" and "MHillasSrc" containers
91// hadronnessName name of container holding the hadronness
92// howMany no.of events to be read from input file
93// maxhadronness maximum value of hadronness to be accepted
94// maxalpha maximum value of |ALPHA| to be accepted
95// maxdist maximum value of DIST (degrees) to be accepted
96//
97
98void CT1EEst(TString inPath, TString fileNameIn,
99 TString outPath, TString energyParName,
100 TString hilName, TString hilSrcName, TString hadronnessName,
101 Int_t howMany,
102 Float_t maxhadronness, Float_t maxalpha, Float_t maxdist,
103 MBinning* binsE, MBinning* binsTheta)
104{
105 cout << "================================================================"
106 << endl;
107 cout << "Start of energy estimation part" << endl;
108 cout << "" << endl;
109 cout << "CT1EEst input values : " << endl;
110 cout << "---------------------- " << endl;
111 cout << "inPath, fileNameIn = "
112 << inPath << ", " << fileNameIn << endl;
113 cout << "outPath, energyParName = "
114 << outPath << ", " << energyParName << endl;
115 cout << "hilName, hilSrcName, hadronnessName = " << hilName << ", "
116 << hilSrcName << ", " << hadronnessName << endl;
117 cout << "howMany, maxhadronness, maxalpha, maxdist = " << howMany << ", "
118 << maxhadronness << ", " << maxalpha << ", " << maxdist << endl;
119 cout << "" << endl;
120
121
122 //==========================================================================
123 // Start of Part 1 (determination of the parameters of the energy estimator)
124 //
125
126 //-----------------------------------------------------------------------
127 // Fill events into a MHMatrix,
128 // to be used later in the minimization by MINUIT
129 //
130
131 MMcEnergyEst Optimize;
132
133 TString inputfile(outPath);
134 inputfile += fileNameIn;
135 Optimize.SetInFile(inputfile);
136
137 TString paramout(outPath);
138 paramout += energyParName;
139 Optimize.SetOutFile(paramout);
140
141
142 MFCT1SelFinal filterhadrons;
143 filterhadrons.SetHadronnessName(hadronnessName);
144 filterhadrons.SetCuts(maxhadronness, maxalpha, maxdist);
145 filterhadrons.SetInverted();
146 Optimize.SetEventFilter(&filterhadrons);
147
148
149 Optimize.SetNevents(howMany);
150
151 //
152 // Find the optimal parameters for the energy estimator:
153 //
154 Optimize.FindParams();
155
156 cout << "--------------------------------------" << endl;
157 cout << "Write parameters of energy estimator onto file '"
158 << paramout << endl;
159 cout << "--------------------------------------" << endl;
160
161 Optimize.Print();
162
163 TFile out(paramout, "recreate");
164 Optimize.SetReadyToSave();
165 Optimize.Write();
166 out.Close();
167
168 //
169 // END of Part 1
170 //==========================================================================
171
172
173
174 //==========================================================================
175 // Start of Part 2 (test quality of energy estimation)
176 //
177 //
178
179 //--------------------------------------------
180 // Read the parameters of the energy estimator
181 //
182
183 TString paramin(outPath);
184 paramin += energyParName;
185
186 TFile enparam(paramin);
187 MMcEnergyEst mcest;
188 mcest.Read("MMcEnergyEst");
189 enparam.Close();
190
191 cout << "--------------------------------------" << endl;
192 cout << "Read parameters of energy estimator from file '"
193 << paramin << endl;
194
195 TArrayD parA(5);
196 TArrayD parB(7);
197
198 for (Int_t i=0; i<parA.GetSize(); i++)
199 parA[i] = mcest.GetCoeff(i);
200 for (Int_t i=0; i<parB.GetSize(); i++)
201 parB[i] = mcest.GetCoeff(i+parA.GetSize());
202
203 //-----------------------------------------------
204 // Create energy estimator task, add necessary containers and
205 // initialize parameters read from file:
206 //
207
208 MEnergyEstParam eest2(hilName);
209 eest2.Add(hilSrcName);
210
211 eest2.SetCoeffA(parA);
212 eest2.SetCoeffB(parB);
213
214 // estimate energy using Daniel's parameters
215 //MEnergyEstParamDanielMkn421 eest2(hilName);
216 //eest2.Add(hilSrcName);
217
218
219 //=======================================================================
220 // Setup the event loop
221 //
222 cout << "--------------------------------------" << endl;
223 cout << "Setup event loop for checking quality of energy estimation" << endl;
224
225
226 MTaskList tlist2;
227 MParList parlist2;
228
229 //-----------------------------------------------
230 // Read events
231 //
232
233 TString inputfile(inPath);
234 inputfile += fileNameIn;
235
236 cout << "Read events from file '" << inputfile << "'" << endl;
237 MReadTree read2("Events");
238 read2.AddFile(inputfile);
239 read2.DisableAutoScheme();
240
241
242 //-----------------------------------------------
243 // Select events
244 //
245
246 cout << "Select events with hadronness < " << maxhadronness
247 << " and |alpha| < " << maxalpha << endl;
248 MFCT1SelFinal hcut2;
249 hcut2.SetHadronnessName(hadronnessName); hcut2;
250 hcut2.SetCuts(maxhadronness, maxalpha, maxdist);
251
252 MContinue cont(&hcut2);
253
254 //-----------------------------------------------
255 // Create Object MHMcEnergyMigration containing useful histograms,
256 // and task MHMcEnergyMigration to fill them:
257 //
258
259 //MHMcEnergyMigration mighist;
260
261 MFillH migfill("MHMcEnergyMigration", "MMcEvt");
262
263 MBinning BinningDE("BinningDE");
264 MBinning BinningImpact("BinningImpact");
265
266 BinningDE.SetEdges(60, -1.2, 1.2);
267 BinningImpact.SetEdges(50, 0., 400.);
268
269 //********************************
270 // Entries in MParList
271
272 parlist2.AddToList(&tlist2);
273 //parlist2.AddToList(&mighist);
274 parlist2.AddToList(binsE);
275 parlist2.AddToList(binsTheta);
276 parlist2.AddToList(&BinningDE);
277 parlist2.AddToList(&BinningImpact);
278
279 //********************************
280 // Entries in MTaskList
281
282 tlist2.AddToList(&read2);
283 tlist2.AddToList(&cont);
284 tlist2.AddToList(&eest2);
285 tlist2.AddToList(&migfill);
286
287 //-----------------------------------------------
288 cout << "Event loop was setup" << endl;
289 MProgressBar bar;
290
291 MEvtLoop evtloop2;
292 evtloop2.SetProgressBar(&bar);
293 evtloop2.SetParList(&parlist2);
294
295 if (!evtloop2.Eventloop())
296 return;
297
298 tlist2.PrintStatistics(0, kTRUE);
299
300 MHMcEnergyMigration &mighist =
301 *(MHMcEnergyMigration*)parlist2->FindObject("MHMcEnergyMigration");
302 if (!&mighist)
303 {
304 cout << "CT1EgyEst.C : object 'MHMcEnergyMigration' not found ... aborting"
305 << endl;
306 return;
307 }
308
309 TString paramout(outPath);
310 paramout += energyParName;
311
312 TFile out2(paramout, "UPDATE");
313 mighist.Write();
314 out2.Close();
315
316 TCanvas *c = new TCanvas;
317 c->cd();
318 mighist.Draw();
319
320 cout << "Quality histograms and migration matrix were added onto the file '" << paramout << endl;
321 cout << endl;
322 cout << "End of energy estimation part" << endl;
323 cout << "================================================================" << endl;
324 //
325 // End of Part 2 (test quality of energy estimation)
326 //==========================================================================
327 //
328
329 return;
330
331}
Note: See TracBrowser for help on using the repository browser.