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

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