source: tags/Mars-V0.9.5/mjoptim/MJOptimizeDisp.cc

Last change on this file was 7413, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 5.8 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, 6/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2005
21!
22!
23\* ======================================================================== */
24/////////////////////////////////////////////////////////////////////////////
25//
26// MJOptimizeDisp
27//
28// Class for otimizing the disp parameters. For more details see
29// MJOptimize.
30//
31// MJOptimizeDisp minimzes the width of the Theta^2 distribution by
32// minimizing <Theta^2> (the RMS^2 of Theta^2), while Theta^2 is
33// calculated as:
34// Theta^2 = d^2 + p^2 - 2*d*p*cos(a)
35// with:
36// d: MHillasSrc.fDist [deg]
37// p: Disp as calculated by the given rule, eg: [0] (a constant)
38// a: MHillasSrc.fAlpha [rad]
39//
40// Example:
41// --------
42// MJOptimizeDisp opt;
43// opt.SetDebug(2);
44// opt.SetOptimizer(MJOptimize::kMigrad);
45// opt.SetNumEvents(20000);
46// opt.EnableTestTrain();
47// opt.AddParameter("1-(MHillas.fWidth/MHillas.fLength)");
48// opt.SetParameter(0, 1, 0, 2);
49// char *r = "[0]*M[0]"; //Rule to calculate disp
50// MStatusDisplay *d = new MStatusDisplay;
51// opt.SetDisplay(d);
52// opt.RunDisp("ganymed-summary.root", r);
53//
54/////////////////////////////////////////////////////////////////////////////
55#include "MJOptimizeDisp.h"
56
57#include "MHMatrix.h"
58
59// environment
60#include "MLog.h"
61#include "MLogManip.h"
62
63// eventloop
64#include "MParList.h"
65#include "MTaskList.h"
66#include "MEvtLoop.h"
67
68// parameters
69#include "MParameters.h"
70#include "MGeomCamMagic.h"
71
72// histograms
73#include "MH3.h"
74#include "MBinning.h"
75#include "../mhflux/MAlphaFitter.h"
76#include "../mhflux/MHThetaSq.h"
77#include "../mtools/MChisqEval.h"
78
79// tasks
80#include "MReadTree.h"
81#include "MMatrixLoop.h"
82#include "MParameterCalc.h"
83#include "MFillH.h"
84
85// filters
86#include "MFDataMember.h"
87
88using namespace std;
89
90ClassImp(MJOptimizeDisp);
91
92//------------------------------------------------------------------------
93//
94// Read all events from file which do match rules and optimize
95// disp estimator.
96//
97Bool_t MJOptimizeDisp::RunDisp(const char *fname, const char *rule, MTask *weights)
98{
99 fLog->Separator("Preparing Disp optimization");
100
101 MParList parlist;
102
103 MParameterI par1("DataType");
104 par1.SetVal(1);
105 parlist.AddToList(&par1);
106
107 MParameterD par2("ThetaSquaredCut");
108 par2.SetVal(0.2);
109 parlist.AddToList(&par2);
110
111 MAlphaFitter fit;
112 fit.SetPolynomOrder(0);
113 fit.SetSignalFitMax(0.8);
114 fit.EnableBackgroundFit(kFALSE);
115 fit.SetSignalFunction(MAlphaFitter::kThetaSq);
116 fit.SetMinimizationStrategy(MAlphaFitter::kGaussSigma);
117 parlist.AddToList(&fit);
118
119 MGeomCamMagic geom; // For GetConvMm2Deg
120 parlist.AddToList(&geom);
121
122 MHMatrix m("M");
123 AddRulesToMatrix(m);
124 parlist.AddToList(&m);
125
126 const Int_t num1 = m.AddColumn("MHillasSrc.fDist*MGeomCam.fConvMm2Deg");
127 const Int_t num2 = m.AddColumn("MHillasSrc.fAlpha*kDegToRad");
128 const Int_t num3 = m.AddColumn("MHillas.fSize");
129
130 MHThetaSq hist;
131 hist.SkipHistTime();
132 hist.SkipHistTheta();
133 //hist.SkipHistEnergy();
134 //hist.ForceUsingSize();
135 hist.InitMapping(&m, 1);
136
137 MFDataMember filter("DataType.fVal", '>', 0.5);
138 fPreCuts.Add(&filter);
139
140 MParameterCalc calc1(rule, "MParameters");
141 calc1.SetNameParameter("Disp");
142 parlist.AddToList(&calc1);
143
144 const char *disp = "Disp.fVal";
145 const char *n1 = Form("M[%d]", num1);
146
147 const char *rule2 = Form("(%s*%s) + (%s*%s) - (2*%s*%s*cos(M[%d]))",
148 n1, n1, disp, disp, n1, disp, num2);
149
150 MParameterCalc calc2(rule2, "MThetaSqCalc");
151 calc2.SetNameParameter("ThetaSquared");
152
153 MReadTree read("Events");
154 // NECESSARY BECAUSE OF MDataFormula GetRules missing
155 read.DisableAutoScheme();
156 if (fname)
157 read.AddFile(fname);
158 else
159 AddSequences(read, fNamesOn);
160 if (!FillMatrix(read, parlist, kTRUE))
161 return kFALSE;
162
163 fPreCuts.Remove(&filter);
164
165 MTaskList tasklist;
166 parlist.Replace(&tasklist);
167
168 MFillH fill(&hist);
169 if (weights)
170 fill.SetWeight();
171
172 MChisqEval eval;
173 eval.SetY1("sqrt(ThetaSquared.fVal)");
174
175 MMatrixLoop loop(&m);
176
177 const char *n3 = Form("M[%d]", num3);
178 MH3 hdisp(n3, "sqrt(ThetaSquared.fVal)");
179 hdisp.SetTitle("\\vartheta^{2} distribution vs. Size:Size [phe]:\\vartheta^{2} [\\circ^{2}]");
180
181 MBinning binsx(100, 10, 100000, "BinningMH3X", "log");
182 MBinning binsy(100, 0, 2, "BinningMH3Y", "lin");
183
184 parlist.AddToList(&binsx);
185 parlist.AddToList(&binsy);
186
187 MFillH fillh2(&hdisp);
188 fillh2.SetDrawOption("blue profx");
189
190 tasklist.AddToList(&loop);
191 tasklist.AddToList(&calc1);
192 tasklist.AddToList(&calc2);
193 if (weights)
194 tasklist.AddToList(weights);
195 tasklist.AddToList(&fill);
196 tasklist.AddToList(&fillh2);
197 tasklist.AddToList(&eval);
198
199 // Optimize with the tasklist in this parameterlist
200 if (!Optimize(parlist))
201 return kFALSE;
202
203 // Print the result
204 *fLog << inf << "Rule: " << rule << endl;
205 hist.GetAlphaFitter().Print("result");
206
207 // Store result if requested
208 TObjArray cont;
209 cont.Add(&calc1);
210 return WriteContainer(cont, fNameOut);
211}
Note: See TracBrowser for help on using the repository browser.