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

Last change on this file since 8154 was 8106, checked in by tbretz, 18 years ago
*** empty log message ***
File size: 12.2 KB
Line 
1/* ======================================================================== *\
2! $Name: not supported by cvs2svn $:$Id: MHThetaSqN.cc,v 1.8 2006-10-17 17:16:00 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)[n1], (*fMatrix)[n1+1] );
215}
216
217Bool_t MHThetaSqN::Fill(const MParContainer *par, const Stat_t weight)
218{
219 const TVector2 mean(GetVec(fHillas->GetMean(), 6));
220 const TVector2 norm(GetVec(fHillas->GetNormAxis(), 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 if (!MHAlpha::Fill(NULL, w))
275 return kFALSE;
276 }
277
278 if (!fOffData)
279 {
280 // off 0.4deg 0.3deg
281 // 5: 0.200 ~86% 0.150
282 // 4: 0.235 ~92% 0.176
283 // 3: 0.283 ~96% 0.212
284 // 2: 0.346 ~98% 0.260
285
286 const Double_t dis = src0.Mod()*TMath::Sin(rad/2);
287 const Double_t cut = TMath::Sqrt(fFit.GetSignalIntegralMax());
288 if (dis<cut)
289 {
290 *fLog << warn << "WARNING - Source regions overlap: distance ";
291 *fLog << dis << " less than theta-sq cut " << cut << "!" << endl;
292 if (dis*1.7<cut*1.0)
293 fCounter[3]++;
294 else
295 if (dis*1.7<cut*1.5)
296 fCounter[2]++;
297 else
298 fCounter[1]++;
299
300 }
301 else
302 fCounter[0]++;
303 }
304
305 return kTRUE;
306}
307
308// --------------------------------------------------------------------------
309//
310// Print execution statistics
311//
312Bool_t MHThetaSqN::Finalize()
313{
314 if (GetNumExecutions()>0)
315 {
316 *fLog << inf << endl;
317 *fLog << GetDescriptor() << " execution statistics:" << endl;
318 PrintSkipped(fCounter[3], "Region distance below 1.0sigma (68.3%)");
319 PrintSkipped(fCounter[2], "Region distance below 1.5sigma (86.6%)");
320 PrintSkipped(fCounter[1], "Region distance below 1.7sigma (91.1%)");
321 *fLog << " " << (int)fCounter[0] << " (" << Form("%5.1f", 100.*fCounter[0]/GetNumExecutions()) << "%) Evts ok!" << endl;
322 //*fLog << " " << (int)fCounter[0] << " (" << Form("%5.1f", 100.*fCounter[0]/GetNumExecutions()) << "%) Evts survived calculation!" << endl;
323 *fLog << endl;
324 }
325
326 return MHAlpha::Finalize();
327}
328
329// --------------------------------------------------------------------------
330//
331// You can use this function if you want to use a MHMatrix instead of
332// MMcEvt. This function adds all necessary columns to the
333// given matrix. Afterward you should fill the matrix with the corresponding
334// data (eg from a file by using MHMatrix::Fill). If you now loop
335// through the matrix (eg using MMatrixLoop) MHHadronness::Fill
336// will take the values from the matrix instead of the containers.
337//
338// It takes fSkipHist* into account!
339//
340void MHThetaSqN::InitMapping(MHMatrix *mat, Int_t type)
341{
342 if (fMatrix)
343 return;
344
345 fMatrix = mat;
346
347 fMap[0] = -1;
348 fMap[1] = -1;
349 fMap[2] = -1;
350 fMap[3] = -1;
351 fMap[4] = -1;
352
353 if (!fSkipHistEnergy)
354 if (type==0)
355 {
356 fMap[1] = fMatrix->AddColumn("MEnergyEst.fVal");
357 fMap[2] = -1;
358 }
359 else
360 {
361 fMap[1] = -1;
362 fMap[2] = fMatrix->AddColumn("MHillas.fSize");
363 }
364
365 if (!fSkipHistTheta)
366 fMap[3] = fMatrix->AddColumn("MPointingPos.fZd");
367
368 fMap[5] = fMatrix->AddColumn("Disp.fVal");
369 fMap[6] = fMatrix->AddColumn("MHillas.fCosDelta");
370 fMap[7] = fMatrix->AddColumn("MHillas.fSinDelta");
371 fMap[8] = fMatrix->AddColumn("MHillas.fMeanX");
372 fMap[9] = fMatrix->AddColumn("MHillas.fMeanY");
373 fMap[10] = fMatrix->AddColumn("MSrcPosCam.fX");
374 fMap[11] = fMatrix->AddColumn("MSrcPosCan.fY");
375
376 // if (!fSkipHistTime)
377 // fMap[4] = fMatrix->AddColumn("MTime.GetAxisTime");
378}
379
380Int_t MHThetaSqN::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
381{
382 Int_t rc = MHAlpha::ReadEnv(env, prefix, print);
383 if (rc==kERROR)
384 return kERROR;
385
386 if (IsEnvDefined(env, prefix, "NumBinsSignal", print))
387 {
388 SetNumBinsSignal(GetEnvValue(env, prefix, "NumBinsSignal", (Int_t)fNumBinsSignal));
389 rc = kTRUE;
390 }
391 if (IsEnvDefined(env, prefix, "NumBinsTotal", print))
392 {
393 SetNumBinsTotal(GetEnvValue(env, prefix, "NumBinsTotal", (Int_t)fNumBinsTotal));
394 rc = kTRUE;
395 }
396 if (IsEnvDefined(env, prefix, "NumOffSourcePos", print))
397 {
398 SetNumOffSourcePos(GetEnvValue(env, prefix, "NumOffSourcePos", (Int_t)fNumOffSourcePos));
399 rc = kTRUE;
400 }
401 if (IsEnvDefined(env, prefix, "DoOffCut", print))
402 {
403 SetDoOffCut(GetEnvValue(env, prefix, "DoOffCut", fDoOffCut));
404 rc = kTRUE;
405 }
406 if (IsEnvDefined(env, prefix, "SignificanceCutLevel", print))
407 {
408 SetSignificanceCutLevel(GetEnvValue(env, prefix, "SignificanceCutLevel", fSignificanceCutLevel));
409 rc = kTRUE;
410 }
411 return rc;
412}
Note: See TracBrowser for help on using the repository browser.