source: trunk/MagicSoft/Mars/mhflux/MHThetaSqN.cc@ 9311

Last change on this file since 9311 was 9153, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 12.1 KB
Line 
1/* ======================================================================== *\
2! $Name: not supported by cvs2svn $:$Id: MHThetaSqN.cc,v 1.10 2008-11-11 11:42:14 tbretz Exp $
3! --------------------------------------------------------------------------
4!
5! *
6! * This file is part of MARS, the MAGIC Analysis and Reconstruction
7! * Software. It is distributed to you in the hope that it can be a useful
8! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
9! * It is distributed WITHOUT ANY WARRANTY.
10! *
11! * Permission to use, copy, modify and distribute this software and its
12! * documentation for any purpose is hereby granted without fee,
13! * provided that the above copyright notice appear in all copies and
14! * that both that copyright notice and this permission notice appear
15! * in supporting documentation. It is provided "as is" without express
16! * or implied warranty.
17! *
18!
19!
20! Author(s): Thomas Bretz, 5/2006 <mailto:tbretz@astro.uni-wuerzburg.de>
21!
22! Copyright: MAGIC Software Development, 2000-2006
23!
24!
25\* ======================================================================== */
26
27//////////////////////////////////////////////////////////////////////////////
28//
29// MHThetaSqN
30//
31// For more detailes see MHAlpha.
32//
33// How to use this class:
34// - in your ganymed.rc add:
35// MJCut.NameHist: MHThetaSqN
36// - setup the number of off-source regions in your ganymed.rc by:
37// MHThetaSqN.NumOffSourcePos: 2
38// - switch on/off whether an off-theta cut should be done:
39// MHThetaSqN.DoOffCut: Yes,No
40// - and if necessary switch off the Theta cut in your Magic cuts:
41// Cut1.ThetaCut: None
42//
43//////////////////////////////////////////////////////////////////////////////
44#include "MHThetaSqN.h"
45
46#include <TVector2.h>
47
48#include "MParameters.h"
49#include "MSrcPosCam.h"
50#include "MGeomCam.h"
51#include "MHMatrix.h"
52#include "MHillas.h"
53
54#include "MBinning.h"
55#include "MParList.h"
56#include "MTaskList.h"
57#include "MParameters.h"
58
59#include "MLog.h"
60#include "MLogManip.h"
61
62ClassImp(MHThetaSqN);
63
64using namespace std;
65
66// --------------------------------------------------------------------------
67//
68// Default Constructor
69//
70MHThetaSqN::MHThetaSqN(const char *name, const char *title)
71 : MHAlpha(name, title), fDisp(0), fSrcPosCam(0), fSignificanceCutLevel(1.7),
72 fNumBinsSignal(3), fNumBinsTotal(75), fNumOffSourcePos(3), fDoOffCut(kTRUE),
73 fCounter(4)
74{
75 //
76 // set the name and title of this object
77 //
78 fName = name ? name : "MHThetaSqN";
79 fTitle = title ? title : "Theta Squared plot";
80
81 fNameParameter = "ThetaSquared";
82
83 fHist.SetName("Theta");
84 fHist.SetTitle("Theta");
85 fHist.SetZTitle("\\vartheta^{2} [deg^{2}]");
86 fHist.SetDirectory(NULL);
87
88 // Main histogram
89 fHistTime.SetName("Theta");
90 fHistTime.SetXTitle("\\vartheta^{2} [deg^{2}]");
91 fHistTime.SetDirectory(NULL);
92
93 MBinning binsa, binse, binst;
94 binsa.SetEdges(75, 0, 1.5);
95 //binsa.SetEdges(arr);
96 binse.SetEdgesLog(15, 10, 100000);
97 binst.SetEdgesASin(67, -0.005, 0.665);
98 binsa.Apply(fHistTime);
99
100 MH::SetBinning(&fHist, &binst, &binse, &binsa);
101
102 fParameter = new MParameterD;
103}
104
105MHThetaSqN::~MHThetaSqN()
106{
107 delete fParameter;
108}
109
110// --------------------------------------------------------------------------
111//
112// Overwrites the binning in Alpha (ThetaSq) with a binning for which
113// the upper edge of the 5th bin (bin=5) fit the signal integration window.
114// In total 75 bins are setup.
115//
116// In case of fOffData!=NULL the binnings are taken later from fOffData anyhow.
117//
118Bool_t MHThetaSqN::SetupFill(const MParList *pl)
119{
120 // Default is from default fitter
121 // if a user defined fitter is in the parlist: use this range
122 MAlphaFitter *fit = (MAlphaFitter*)pl->FindObject("MAlphaFitter");
123 if (!fit)
124 fit = &fFit;
125
126 MParameterD *cut = (MParameterD*)pl->FindObject("ThetaSquaredCut", "MParameterD");
127 if (cut)
128 fit->SetSignalIntegralMax(cut->GetVal());
129
130 // Get Histogram binnings
131 MBinning binst, binse;
132 binst.SetEdges(fHist, 'x');
133 binse.SetEdges(fHist, 'y');
134
135 // Calculate bining which fits alpha-cut
136 const Double_t intmax = fit->GetSignalIntegralMax();
137 const UInt_t nbins = fNumBinsTotal;
138 const UInt_t nsig = fNumBinsSignal;
139
140 MBinning binsa(nbins, 0, nbins*intmax/nsig);
141
142 // Apply binning
143 binsa.Apply(fHistTime);
144 MH::SetBinning(&fHist, &binst, &binse, &binsa);
145
146 // Remark: Binnings might be overwritten in MHAlpha::SetupFill
147 if (!MHAlpha::SetupFill(pl))
148 return kFALSE;
149
150 fDisp = (MParameterD*)pl->FindObject("Disp", "MParameterD");
151 if (!fDisp)
152 {
153 *fLog << err << "Disp [MParameterD] not found... abort." << endl;
154 return kFALSE;
155 }
156 fHillas = (MHillas*)pl->FindObject("MHillas");
157 if (!fHillas)
158 {
159 *fLog << err << "MHillas not found... abort." << endl;
160 return kFALSE;
161 }
162 fSrcPosCam = (MSrcPosCam*)pl->FindObject("MSrcPosCam");
163 if (!fSrcPosCam)
164 {
165 *fLog << err << "MSrcPosCam not found... abort." << endl;
166 return kFALSE;
167 }
168
169 MGeomCam *geom = (MGeomCam*)pl->FindObject("MGeomCam");
170 if (!geom)
171 {
172 *fLog << err << "MGeomCam not found... abort." << endl;
173 return kFALSE;
174 }
175 fMm2Deg = geom->GetConvMm2Deg();
176
177 if (fFit.GetScaleMode()==MAlphaFitter::kNone)
178 fFit.SetScaleUser(1./fNumOffSourcePos);
179
180 fThetaSqCut = fSignificanceCutLevel*fFit.GetSignalIntegralMax()/1.7;
181
182 fCounter.Reset();
183
184 return kTRUE;
185}
186
187// --------------------------------------------------------------------------
188//
189// Return the value from fParemeter which should be filled into the plots
190//
191Double_t MHThetaSqN::GetVal() const
192{
193 return static_cast<const MParameterD*>(fParameter)->GetVal();
194}
195
196// --------------------------------------------------------------------------
197//
198// Abbreviation to set the value used by MHAlpha to fill the histogram
199//
200void MHThetaSqN::SetVal(Double_t val)
201{
202 static_cast<MParameterD*>(fParameter)->SetVal(val);
203}
204
205// --------------------------------------------------------------------------
206//
207// Abbreviation to set the value used by MHAlpha to fill the histogram
208//
209TVector2 MHThetaSqN::GetVec(const TVector2 &v, Int_t n1) const
210{
211 if (!fMatrix)
212 return v;
213
214 return TVector2( (*fMatrix)[fMap[n1]], (*fMatrix)[fMap[n1+1]] );
215}
216
217Int_t MHThetaSqN::Fill(const MParContainer *par, const Stat_t weight)
218{
219 const TVector2 norm(GetVec(fHillas->GetNormAxis(), 6));
220 const TVector2 mean(GetVec(fHillas->GetMean(), 8));
221
222 const TVector2 org = mean*fMm2Deg + norm*fDisp->GetVal();
223
224 // In the case we are filling off-data src0 contains the anti-source position
225 TVector2 src0(GetVec(fSrcPosCam->GetXY(), 10)*fMm2Deg);
226 if (!fOffData)
227 src0 *= -1;
228
229 const UInt_t n = fNumOffSourcePos+1;
230 const Float_t rad = TMath::TwoPi()/n;
231
232 // Calculate distance (theta-sq) to all (off-)source regions
233 TArrayD dist(n);
234 for (UInt_t i=0; i<n; i++)
235 {
236 const TVector2 src = const_cast<TVector2&>(src0).Rotate(i*rad);
237 dist[i] = (src-org).Mod2();
238 }
239
240 // Processing off-data
241 // Check if event's origin is in the on-regions
242 if (!fOffData && fDoOffCut && dist[0]<fThetaSqCut)
243 return kTRUE;
244
245 for (UInt_t i=0; i<n; i++)
246 {
247 // off: is in src region on: is in off regions
248 /// if (!fOffData && i==0) || (fOffData && i!=0)
249 if ((bool)fOffData ^ (i==0) )
250 continue;
251
252 Stat_t w = weight;
253
254 // Processing on-data
255 if (fOffData && fDoOffCut)
256 {
257 /*
258 static int cnt=0;
259 if (dist[1+(cnt++%fNumOffSourcePos)]<fFit.GetSignalIntegralMax())
260 continue;
261 */
262
263 // Check if event's origin is in one of the off-regions
264 for (UInt_t j=1; j<n; j++)
265 if (dist[j]<fThetaSqCut)
266 {
267 w *= (float)(fNumOffSourcePos-1)/fNumOffSourcePos;
268 break;
269 }
270 }
271
272 SetVal(dist[i]);
273
274 const Int_t rc = MHAlpha::Fill(NULL, w);
275 if (rc!=kTRUE)
276 return rc;
277 }
278
279 if (!fOffData)
280 {
281 // off 0.4deg 0.3deg
282 // 5: 0.200 ~86% 0.150
283 // 4: 0.235 ~92% 0.176
284 // 3: 0.283 ~96% 0.212
285 // 2: 0.346 ~98% 0.260
286
287 const Double_t dis = src0.Mod()*TMath::Sin(rad/2);
288 const Double_t cut = TMath::Sqrt(fFit.GetSignalIntegralMax());
289 if (dis<cut)
290 {
291 *fLog << warn << "WARNING - Source regions overlap: distance ";
292 *fLog << dis << " less than theta-sq cut " << cut << "!" << endl;
293 if (dis*1.7<cut*1.0)
294 fCounter[3]++;
295 else
296 if (dis*1.7<cut*1.5)
297 fCounter[2]++;
298 else
299 fCounter[1]++;
300
301 }
302 else
303 fCounter[0]++;
304 }
305
306 return kTRUE;
307}
308
309// --------------------------------------------------------------------------
310//
311// Print execution statistics
312//
313Bool_t MHThetaSqN::Finalize()
314{
315 if (GetNumExecutions()>0)
316 {
317 *fLog << inf << endl;
318 *fLog << GetDescriptor() << " execution statistics:" << endl;
319 PrintSkipped(fCounter[3], "Region distance below 1.0sigma (68.3%)");
320 PrintSkipped(fCounter[2], "Region distance below 1.5sigma (86.6%)");
321 PrintSkipped(fCounter[1], "Region distance below 1.7sigma (91.1%)");
322 *fLog << " " << (int)fCounter[0] << " (" << Form("%5.1f", 100.*fCounter[0]/GetNumExecutions()) << "%) Evts ok!" << endl;
323 //*fLog << " " << (int)fCounter[0] << " (" << Form("%5.1f", 100.*fCounter[0]/GetNumExecutions()) << "%) Evts survived calculation!" << endl;
324 *fLog << endl;
325 }
326
327 return MHAlpha::Finalize();
328}
329
330// --------------------------------------------------------------------------
331//
332// You can use this function if you want to use a MHMatrix instead of
333// MMcEvt. This function adds all necessary columns to the
334// given matrix. Afterward you should fill the matrix with the corresponding
335// data (eg from a file by using MHMatrix::Fill). If you now loop
336// through the matrix (eg using MMatrixLoop) MHHadronness::Fill
337// will take the values from the matrix instead of the containers.
338//
339// It takes fSkipHist* into account!
340//
341void MHThetaSqN::InitMapping(MHMatrix *mat, Int_t type)
342{
343 if (fMatrix)
344 return;
345
346 fMatrix = mat;
347
348 fMap[0] = -1;
349 fMap[1] = -1;
350 fMap[2] = -1;
351 fMap[3] = -1;
352 fMap[4] = -1;
353
354 if (!fSkipHistEnergy)
355 if (type==0)
356 {
357 fMap[1] = fMatrix->AddColumn("MEnergyEst.fVal");
358 fMap[2] = -1;
359 }
360 else
361 {
362 fMap[1] = -1;
363 fMap[2] = fMatrix->AddColumn("MHillas.fSize");
364 }
365
366 if (!fSkipHistTheta)
367 fMap[3] = fMatrix->AddColumn("MPointingPos.fZd");
368
369 fMap[5] = -1;
370 fMap[6] = fMatrix->AddColumn("MHillas.fCosDelta");
371 fMap[7] = fMatrix->AddColumn("MHillas.fSinDelta");
372 fMap[8] = fMatrix->AddColumn("MHillas.fMeanX");
373 fMap[9] = fMatrix->AddColumn("MHillas.fMeanY");
374 fMap[10] = fMatrix->AddColumn("MSrcPosCam.fX");
375 fMap[11] = fMatrix->AddColumn("MSrcPosCam.fY");
376}
377
378Int_t MHThetaSqN::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
379{
380 Int_t rc = MHAlpha::ReadEnv(env, prefix, print);
381 if (rc==kERROR)
382 return kERROR;
383
384 if (IsEnvDefined(env, prefix, "NumBinsSignal", print))
385 {
386 SetNumBinsSignal(GetEnvValue(env, prefix, "NumBinsSignal", (Int_t)fNumBinsSignal));
387 rc = kTRUE;
388 }
389 if (IsEnvDefined(env, prefix, "NumBinsTotal", print))
390 {
391 SetNumBinsTotal(GetEnvValue(env, prefix, "NumBinsTotal", (Int_t)fNumBinsTotal));
392 rc = kTRUE;
393 }
394 if (IsEnvDefined(env, prefix, "NumOffSourcePos", print))
395 {
396 SetNumOffSourcePos(GetEnvValue(env, prefix, "NumOffSourcePos", (Int_t)fNumOffSourcePos));
397 rc = kTRUE;
398 }
399 if (IsEnvDefined(env, prefix, "DoOffCut", print))
400 {
401 SetDoOffCut(GetEnvValue(env, prefix, "DoOffCut", fDoOffCut));
402 rc = kTRUE;
403 }
404 if (IsEnvDefined(env, prefix, "SignificanceCutLevel", print))
405 {
406 SetSignificanceCutLevel(GetEnvValue(env, prefix, "SignificanceCutLevel", fSignificanceCutLevel));
407 rc = kTRUE;
408 }
409 return rc;
410}
Note: See TracBrowser for help on using the repository browser.