source: trunk/MagicSoft/Mars/mimage/MHillasExt.cc@ 4054

Last change on this file since 4054 was 3682, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 7.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 12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
19! Author(s): Rudolf Bock 10/2001 <mailto:Rudolf.Bock@cern.ch>
20! Author(s): Wolfgang Wittek 06/2002 <mailto:wittek@mppmu.mpg.de>
21!
22! Copyright: MAGIC Software Development, 2000-2004
23!
24!
25\* ======================================================================== */
26
27/////////////////////////////////////////////////////////////////////////////
28//
29// MHillasExt
30//
31// Storage Container for extended image parameters
32//
33// extended image parameters
34//
35//
36// Version 1:
37// ----------
38// fConc ratio of sum of two highest pixels over fSize
39// fConc1 ratio of highest pixel over fSize
40// fAsym distance from highest pixel to center, projected onto major axis
41// fM3Long third moment along major axis
42// fM3Trans third moment along minor axis
43//
44// Version 2:
45// ----------
46// fConc removed
47// fConc1 removed
48//
49// Version 3:
50// ----------
51// fMaxDist added. Distance between center and most distant used pixel
52//
53//
54// WARNING: Before you can use fAsym, fM3Long and fM3Trans you must
55// multiply by the sign of MHillasSrc::fCosDeltaAlpha
56//
57////////////////////////////////////////////////////////////////////////////
58/*
59 // fAsymna d/(d na) of ( sum(x*q^na)/sum(q^na), sum(y*q^na)/sum(q^na) )
60 // projected onto the major axis
61 // fAsym0 (F-B)/(F+B) along the major axis
62 */
63#include "MHillasExt.h"
64
65#include <fstream>
66#include <TArrayF.h>
67
68#include "MGeomPix.h"
69#include "MGeomCam.h"
70
71#include "MCerPhotPix.h"
72#include "MCerPhotEvt.h"
73
74#include "MLog.h"
75#include "MLogManip.h"
76
77#include "MHillas.h"
78
79ClassImp(MHillasExt);
80
81using namespace std;
82
83// -------------------------------------------------------------------------
84//
85// Default constructor.
86//
87MHillasExt::MHillasExt(const char *name, const char *title)
88{
89 fName = name ? name : "MHillasExt";
90 fTitle = title ? title : "Storage container for extended parameter set of one event";
91
92 Reset();
93}
94
95// -------------------------------------------------------------------------
96//
97void MHillasExt::Reset()
98{
99 fAsym = 0;
100 fM3Long = 0;
101 fM3Trans = 0;
102
103 fMaxDist = 0;
104}
105
106// -------------------------------------------------------------------------
107//
108// Print contents of MHillasExt to *fLog.
109//
110void MHillasExt::Print(Option_t *) const
111{
112 *fLog << all;
113 *fLog << "Extended Image Parameters (" << GetName() << ")" << endl;
114 *fLog << " - Asymmetry = " << fAsym << endl;
115 *fLog << " - 3.Moment Long [mm] = " << fM3Long << endl;
116 *fLog << " - 3.Moment Trans [mm] = " << fM3Trans << endl;
117 *fLog << " - Max.Dist [mm] = " << fMaxDist << endl;
118}
119
120// -------------------------------------------------------------------------
121//
122// Print contents of MHillasExt to *fLog, depending on the geometry in
123// units of deg.
124//
125void MHillasExt::Print(const MGeomCam &geom) const
126{
127 *fLog << all;
128 *fLog << "Extended Image Parameters (" << GetName() << ")" << endl;
129 *fLog << " - Asymmetry = " << fAsym *geom.GetConvMm2Deg() << endl;
130 *fLog << " - 3.Moment Long [deg] = " << fM3Long *geom.GetConvMm2Deg() << endl;
131 *fLog << " - 3.Moment Trans [deg] = " << fM3Trans*geom.GetConvMm2Deg() << endl;
132 *fLog << " - Max.Dist [deg] = " << fMaxDist*geom.GetConvMm2Deg() << endl;
133}
134
135// -------------------------------------------------------------------------
136//
137// calculation of additional parameters based on the camera geometry
138// and the cerenkov photon event
139//
140Int_t MHillasExt::Calc(const MGeomCam &geom, const MCerPhotEvt &evt, const MHillas &hil)
141{
142 //
143 // calculate the additional image parameters
144 // --------------------------------------------
145 //
146 // loop to get third moments along ellipse axes and two max pixels
147 //
148 // For the moments double precision is used to make sure, that
149 // the complex matrix multiplication and sum is evaluated correctly.
150 //
151 Double_t m3x = 0;
152 Double_t m3y = 0;
153
154 const UInt_t npixevt = evt.GetNumPixels();
155
156 Int_t maxpixid = 0;
157 Float_t maxpix = 0;
158 Float_t maxdist = 0;
159
160 for (UInt_t i=0; i<npixevt; i++)
161 {
162 const MCerPhotPix &pix = evt[i];
163 if (!pix.IsPixelUsed())
164 continue;
165
166 const Int_t pixid = pix.GetPixId();
167
168 const MGeomPix &gpix = geom[pixid];
169 const Double_t dx = gpix.GetX() - hil.GetMeanX(); // [mm]
170 const Double_t dy = gpix.GetY() - hil.GetMeanY(); // [mm]
171
172 Double_t nphot = pix.GetNumPhotons(); // [1]
173
174 const Double_t dzx = hil.GetCosDelta()*dx + hil.GetSinDelta()*dy; // [mm]
175 const Double_t dzy = -hil.GetSinDelta()*dx + hil.GetCosDelta()*dy; // [mm]
176
177 const Double_t dist = dx*dx+dy*dy;
178 if (TMath::Abs(dist)>TMath::Abs(maxdist))
179 maxdist = dzx<0 ? -dist : dist; // [mm^2]
180
181 m3x += nphot * dzx*dzx*dzx; // [mm^3]
182 m3y += nphot * dzy*dzy*dzy; // [mm^3]
183
184 //
185 // Now we are working on absolute values of nphot, which
186 // must take pixel size into account
187 //
188 nphot *= geom.GetPixRatio(pixid);
189
190 if (nphot>maxpix)
191 {
192 maxpix = nphot; // [1]
193 maxpixid = pixid;
194 continue; // [1]
195 }
196 }
197
198 const MGeomPix &maxp = geom[maxpixid];
199
200 fAsym = (hil.GetMeanX()-maxp.GetX())*hil.GetCosDelta() +
201 (hil.GetMeanY()-maxp.GetY())*hil.GetSinDelta(); // [mm]
202
203 //
204 // Third moments along axes get normalized
205 //
206 m3x /= hil.GetSize();
207 m3y /= hil.GetSize();
208
209 fM3Long = m3x<0 ? -pow(-m3x, 1./3) : pow(m3x, 1./3); // [mm]
210 fM3Trans = m3y<0 ? -pow(-m3y, 1./3) : pow(m3y, 1./3); // [mm]
211
212 const Double_t md = TMath::Sqrt(TMath::Abs(maxdist));
213 fMaxDist = maxdist<0 ? -md : md; // [mm]
214
215 SetReadyToSave();
216
217 return 0;
218}
219
220// --------------------------------------------------------------------------
221//
222// This function is ment for special usage, please never try to set
223// values via this function
224//
225void MHillasExt::Set(const TArrayF &arr)
226{
227 if (arr.GetSize() != 4)
228 return;
229
230 fAsym = arr.At(0);
231 fM3Long = arr.At(1);
232 fM3Trans = arr.At(2);
233 fMaxDist = arr.At(3);
234}
Note: See TracBrowser for help on using the repository browser.