source: trunk/Mars/mhflux/MHThetaSqN.cc@ 13408

Last change on this file since 13408 was 10166, checked in by tbretz, 14 years ago
Removed the old obsolete cvs header line.
File size: 11.9 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 fGeom = (MGeomCam*)pl->FindObject("MGeomCam");
168 if (!fGeom)
169 {
170 *fLog << err << "MGeomCam not found... abort." << endl;
171 return kFALSE;
172 }
173
174 if (fFit.GetScaleMode()==MAlphaFitter::kNone)
175 fFit.SetScaleUser(1./fNumOffSourcePos);
176
177 fThetaSqCut = fSignificanceCutLevel*fFit.GetSignalIntegralMax()/1.7;
178
179 fCounter.Reset();
180
181 return kTRUE;
182}
183
184// --------------------------------------------------------------------------
185//
186// Return the value from fParemeter which should be filled into the plots
187//
188Double_t MHThetaSqN::GetVal() const
189{
190 return static_cast<const MParameterD*>(fParameter)->GetVal();
191}
192
193// --------------------------------------------------------------------------
194//
195// Abbreviation to set the value used by MHAlpha to fill the histogram
196//
197void MHThetaSqN::SetVal(Double_t val)
198{
199 static_cast<MParameterD*>(fParameter)->SetVal(val);
200}
201
202// --------------------------------------------------------------------------
203//
204// Abbreviation to set the value used by MHAlpha to fill the histogram
205//
206TVector2 MHThetaSqN::GetVec(const TVector2 &v, Int_t n1) const
207{
208 if (!fMatrix)
209 return v;
210
211 return TVector2( (*fMatrix)[fMap[n1]], (*fMatrix)[fMap[n1+1]] );
212}
213
214Int_t MHThetaSqN::Fill(const MParContainer *par, const Stat_t weight)
215{
216 const Double_t fMm2Deg = fGeom->GetConvMm2Deg();
217
218 const TVector2 norm(GetVec(fHillas->GetNormAxis(), 6));
219 const TVector2 mean(GetVec(fHillas->GetMean(), 8));
220
221 const TVector2 org = mean*fMm2Deg + norm*fDisp->GetVal();
222
223 // In the case we are filling off-data src0 contains the anti-source position
224 TVector2 src0(GetVec(fSrcPosCam->GetXY(), 10)*fMm2Deg);
225 if (!fOffData)
226 src0 *= -1;
227
228 const UInt_t n = fNumOffSourcePos+1;
229 const Float_t rad = TMath::TwoPi()/n;
230
231 // Calculate distance (theta-sq) to all (off-)source regions
232 TArrayD dist(n);
233 for (UInt_t i=0; i<n; i++)
234 {
235 const TVector2 src = const_cast<TVector2&>(src0).Rotate(i*rad);
236 dist[i] = (src-org).Mod2();
237 }
238
239 // Processing off-data
240 // Check if event's origin is in the on-regions
241 if (!fOffData && fDoOffCut && dist[0]<fThetaSqCut)
242 return kTRUE;
243
244 for (UInt_t i=0; i<n; i++)
245 {
246 // off: is in src region on: is in off regions
247 /// if (!fOffData && i==0) || (fOffData && i!=0)
248 if ((bool)fOffData ^ (i==0) )
249 continue;
250
251 Stat_t w = weight;
252
253 // Processing on-data
254 if (fOffData && fDoOffCut)
255 {
256 /*
257 static int cnt=0;
258 if (dist[1+(cnt++%fNumOffSourcePos)]<fFit.GetSignalIntegralMax())
259 continue;
260 */
261
262 // Check if event's origin is in one of the off-regions
263 for (UInt_t j=1; j<n; j++)
264 if (dist[j]<fThetaSqCut)
265 {
266 w *= (float)(fNumOffSourcePos-1)/fNumOffSourcePos;
267 break;
268 }
269 }
270
271 SetVal(dist[i]);
272
273 const Int_t rc = MHAlpha::Fill(NULL, w);
274 if (rc!=kTRUE)
275 return rc;
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 {
355 fMap[1] = type==0 ? fMatrix->AddColumn("MEnergyEst.fVal") : -1;
356 fMap[2] = type==0 ? -1 : fMatrix->AddColumn("MHillas.fSize");
357 }
358
359 if (!fSkipHistTheta)
360 fMap[3] = fMatrix->AddColumn("MPointingPos.fZd");
361
362 fMap[5] = -1;
363 fMap[6] = fMatrix->AddColumn("MHillas.fCosDelta");
364 fMap[7] = fMatrix->AddColumn("MHillas.fSinDelta");
365 fMap[8] = fMatrix->AddColumn("MHillas.fMeanX");
366 fMap[9] = fMatrix->AddColumn("MHillas.fMeanY");
367 fMap[10] = fMatrix->AddColumn("MSrcPosCam.fX");
368 fMap[11] = fMatrix->AddColumn("MSrcPosCam.fY");
369}
370
371Int_t MHThetaSqN::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
372{
373 Int_t rc = MHAlpha::ReadEnv(env, prefix, print);
374 if (rc==kERROR)
375 return kERROR;
376
377 if (IsEnvDefined(env, prefix, "NumBinsSignal", print))
378 {
379 SetNumBinsSignal(GetEnvValue(env, prefix, "NumBinsSignal", (Int_t)fNumBinsSignal));
380 rc = kTRUE;
381 }
382 if (IsEnvDefined(env, prefix, "NumBinsTotal", print))
383 {
384 SetNumBinsTotal(GetEnvValue(env, prefix, "NumBinsTotal", (Int_t)fNumBinsTotal));
385 rc = kTRUE;
386 }
387 if (IsEnvDefined(env, prefix, "NumOffSourcePos", print))
388 {
389 SetNumOffSourcePos(GetEnvValue(env, prefix, "NumOffSourcePos", (Int_t)fNumOffSourcePos));
390 rc = kTRUE;
391 }
392 if (IsEnvDefined(env, prefix, "DoOffCut", print))
393 {
394 SetDoOffCut(GetEnvValue(env, prefix, "DoOffCut", fDoOffCut));
395 rc = kTRUE;
396 }
397 if (IsEnvDefined(env, prefix, "SignificanceCutLevel", print))
398 {
399 SetSignificanceCutLevel(GetEnvValue(env, prefix, "SignificanceCutLevel", fSignificanceCutLevel));
400 rc = kTRUE;
401 }
402 return rc;
403}
Note: See TracBrowser for help on using the repository browser.