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

Last change on this file since 9344 was 9341, 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.11 2009-02-15 15:43:30 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 fGeom = (MGeomCam*)pl->FindObject("MGeomCam");
170 if (!fGeom)
171 {
172 *fLog << err << "MGeomCam not found... abort." << endl;
173 return kFALSE;
174 }
175
176 if (fFit.GetScaleMode()==MAlphaFitter::kNone)
177 fFit.SetScaleUser(1./fNumOffSourcePos);
178
179 fThetaSqCut = fSignificanceCutLevel*fFit.GetSignalIntegralMax()/1.7;
180
181 fCounter.Reset();
182
183 return kTRUE;
184}
185
186// --------------------------------------------------------------------------
187//
188// Return the value from fParemeter which should be filled into the plots
189//
190Double_t MHThetaSqN::GetVal() const
191{
192 return static_cast<const MParameterD*>(fParameter)->GetVal();
193}
194
195// --------------------------------------------------------------------------
196//
197// Abbreviation to set the value used by MHAlpha to fill the histogram
198//
199void MHThetaSqN::SetVal(Double_t val)
200{
201 static_cast<MParameterD*>(fParameter)->SetVal(val);
202}
203
204// --------------------------------------------------------------------------
205//
206// Abbreviation to set the value used by MHAlpha to fill the histogram
207//
208TVector2 MHThetaSqN::GetVec(const TVector2 &v, Int_t n1) const
209{
210 if (!fMatrix)
211 return v;
212
213 return TVector2( (*fMatrix)[fMap[n1]], (*fMatrix)[fMap[n1+1]] );
214}
215
216Int_t MHThetaSqN::Fill(const MParContainer *par, const Stat_t weight)
217{
218 const Double_t fMm2Deg = fGeom->GetConvMm2Deg();
219
220 const TVector2 norm(GetVec(fHillas->GetNormAxis(), 6));
221 const TVector2 mean(GetVec(fHillas->GetMean(), 8));
222
223 const TVector2 org = mean*fMm2Deg + norm*fDisp->GetVal();
224
225 // In the case we are filling off-data src0 contains the anti-source position
226 TVector2 src0(GetVec(fSrcPosCam->GetXY(), 10)*fMm2Deg);
227 if (!fOffData)
228 src0 *= -1;
229
230 const UInt_t n = fNumOffSourcePos+1;
231 const Float_t rad = TMath::TwoPi()/n;
232
233 // Calculate distance (theta-sq) to all (off-)source regions
234 TArrayD dist(n);
235 for (UInt_t i=0; i<n; i++)
236 {
237 const TVector2 src = const_cast<TVector2&>(src0).Rotate(i*rad);
238 dist[i] = (src-org).Mod2();
239 }
240
241 // Processing off-data
242 // Check if event's origin is in the on-regions
243 if (!fOffData && fDoOffCut && dist[0]<fThetaSqCut)
244 return kTRUE;
245
246 for (UInt_t i=0; i<n; i++)
247 {
248 // off: is in src region on: is in off regions
249 /// if (!fOffData && i==0) || (fOffData && i!=0)
250 if ((bool)fOffData ^ (i==0) )
251 continue;
252
253 Stat_t w = weight;
254
255 // Processing on-data
256 if (fOffData && fDoOffCut)
257 {
258 /*
259 static int cnt=0;
260 if (dist[1+(cnt++%fNumOffSourcePos)]<fFit.GetSignalIntegralMax())
261 continue;
262 */
263
264 // Check if event's origin is in one of the off-regions
265 for (UInt_t j=1; j<n; j++)
266 if (dist[j]<fThetaSqCut)
267 {
268 w *= (float)(fNumOffSourcePos-1)/fNumOffSourcePos;
269 break;
270 }
271 }
272
273 SetVal(dist[i]);
274
275 const Int_t rc = MHAlpha::Fill(NULL, w);
276 if (rc!=kTRUE)
277 return rc;
278 }
279
280 if (!fOffData)
281 {
282 // off 0.4deg 0.3deg
283 // 5: 0.200 ~86% 0.150
284 // 4: 0.235 ~92% 0.176
285 // 3: 0.283 ~96% 0.212
286 // 2: 0.346 ~98% 0.260
287
288 const Double_t dis = src0.Mod()*TMath::Sin(rad/2);
289 const Double_t cut = TMath::Sqrt(fFit.GetSignalIntegralMax());
290 if (dis<cut)
291 {
292 *fLog << warn << "WARNING - Source regions overlap: distance ";
293 *fLog << dis << " less than theta-sq cut " << cut << "!" << endl;
294 if (dis*1.7<cut*1.0)
295 fCounter[3]++;
296 else
297 if (dis*1.7<cut*1.5)
298 fCounter[2]++;
299 else
300 fCounter[1]++;
301
302 }
303 else
304 fCounter[0]++;
305 }
306
307 return kTRUE;
308}
309
310// --------------------------------------------------------------------------
311//
312// Print execution statistics
313//
314Bool_t MHThetaSqN::Finalize()
315{
316 if (GetNumExecutions()>0)
317 {
318 *fLog << inf << endl;
319 *fLog << GetDescriptor() << " execution statistics:" << endl;
320 PrintSkipped(fCounter[3], "Region distance below 1.0sigma (68.3%)");
321 PrintSkipped(fCounter[2], "Region distance below 1.5sigma (86.6%)");
322 PrintSkipped(fCounter[1], "Region distance below 1.7sigma (91.1%)");
323 *fLog << " " << (int)fCounter[0] << " (" << Form("%5.1f", 100.*fCounter[0]/GetNumExecutions()) << "%) Evts ok!" << endl;
324 //*fLog << " " << (int)fCounter[0] << " (" << Form("%5.1f", 100.*fCounter[0]/GetNumExecutions()) << "%) Evts survived calculation!" << endl;
325 *fLog << endl;
326 }
327
328 return MHAlpha::Finalize();
329}
330
331// --------------------------------------------------------------------------
332//
333// You can use this function if you want to use a MHMatrix instead of
334// MMcEvt. This function adds all necessary columns to the
335// given matrix. Afterward you should fill the matrix with the corresponding
336// data (eg from a file by using MHMatrix::Fill). If you now loop
337// through the matrix (eg using MMatrixLoop) MHHadronness::Fill
338// will take the values from the matrix instead of the containers.
339//
340// It takes fSkipHist* into account!
341//
342void MHThetaSqN::InitMapping(MHMatrix *mat, Int_t type)
343{
344 if (fMatrix)
345 return;
346
347 fMatrix = mat;
348
349 fMap[0] = -1;
350 fMap[1] = -1;
351 fMap[2] = -1;
352 fMap[3] = -1;
353 fMap[4] = -1;
354
355 if (!fSkipHistEnergy)
356 if (type==0)
357 {
358 fMap[1] = fMatrix->AddColumn("MEnergyEst.fVal");
359 fMap[2] = -1;
360 }
361 else
362 {
363 fMap[1] = -1;
364 fMap[2] = fMatrix->AddColumn("MHillas.fSize");
365 }
366
367 if (!fSkipHistTheta)
368 fMap[3] = fMatrix->AddColumn("MPointingPos.fZd");
369
370 fMap[5] = -1;
371 fMap[6] = fMatrix->AddColumn("MHillas.fCosDelta");
372 fMap[7] = fMatrix->AddColumn("MHillas.fSinDelta");
373 fMap[8] = fMatrix->AddColumn("MHillas.fMeanX");
374 fMap[9] = fMatrix->AddColumn("MHillas.fMeanY");
375 fMap[10] = fMatrix->AddColumn("MSrcPosCam.fX");
376 fMap[11] = fMatrix->AddColumn("MSrcPosCam.fY");
377}
378
379Int_t MHThetaSqN::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
380{
381 Int_t rc = MHAlpha::ReadEnv(env, prefix, print);
382 if (rc==kERROR)
383 return kERROR;
384
385 if (IsEnvDefined(env, prefix, "NumBinsSignal", print))
386 {
387 SetNumBinsSignal(GetEnvValue(env, prefix, "NumBinsSignal", (Int_t)fNumBinsSignal));
388 rc = kTRUE;
389 }
390 if (IsEnvDefined(env, prefix, "NumBinsTotal", print))
391 {
392 SetNumBinsTotal(GetEnvValue(env, prefix, "NumBinsTotal", (Int_t)fNumBinsTotal));
393 rc = kTRUE;
394 }
395 if (IsEnvDefined(env, prefix, "NumOffSourcePos", print))
396 {
397 SetNumOffSourcePos(GetEnvValue(env, prefix, "NumOffSourcePos", (Int_t)fNumOffSourcePos));
398 rc = kTRUE;
399 }
400 if (IsEnvDefined(env, prefix, "DoOffCut", print))
401 {
402 SetDoOffCut(GetEnvValue(env, prefix, "DoOffCut", fDoOffCut));
403 rc = kTRUE;
404 }
405 if (IsEnvDefined(env, prefix, "SignificanceCutLevel", print))
406 {
407 SetSignificanceCutLevel(GetEnvValue(env, prefix, "SignificanceCutLevel", fSignificanceCutLevel));
408 rc = kTRUE;
409 }
410 return rc;
411}
Note: See TracBrowser for help on using the repository browser.