source: trunk/MagicSoft/Mars/manalysis/MSigmabar.cc@ 3881

Last change on this file since 3881 was 3530, checked in by wittek, 21 years ago
*** empty log message ***
File size: 7.5 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): Robert Wagner, 10/2002 <mailto:magicsoft@rwagner.de>
19! Author(s): Wolfgang Wittek, 01/2003 <mailto:wittek@mppmu.mpg.de>
20! Author(s): Thomas Bretz, 04/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
21!
22! Copyright: MAGIC Software Development, 2003
23!
24!
25\* ======================================================================== */
26
27/////////////////////////////////////////////////////////////////////////////
28// //
29// MSigmabar //
30// //
31// This is the storage container to hold information about //
32// "average" pedestal sigmas //
33// //
34// In calculating averages all sigmas are normalized to the area of pixel 0//
35// //
36/////////////////////////////////////////////////////////////////////////////
37#include "MSigmabar.h"
38
39#include <TMath.h>
40
41#include "MLog.h"
42#include "MLogManip.h"
43
44#include "MParList.h"
45
46#include "MGeomCam.h"
47#include "MGeomPix.h"
48
49#include "MCerPhotEvt.h"
50#include "MCerPhotPix.h"
51
52#include "MPedPhotCam.h"
53#include "MPedPhotPix.h"
54
55ClassImp(MSigmabar);
56
57using namespace std;
58
59// --------------------------------------------------------------------------
60//
61MSigmabar::MSigmabar(const char *name, const char *title)
62{
63 fName = name ? name : "MSigmabar";
64 fTitle = title ? title : "Storage container for Sigmabar";
65}
66
67// --------------------------------------------------------------------------
68//
69MSigmabar::~MSigmabar()
70{
71 // do nothing special.
72}
73
74void MSigmabar::Reset()
75{
76 fSigmabar = -1;
77 fInnerPixels = -1;
78 fOuterPixels = -1;
79 fSigmabarInner = -1;
80 fSigmabarOuter = -1;
81
82 memset(fSigmabarSector, 0, sizeof(fSigmabarSector));
83 memset(fSigmabarInnerSector, 0, sizeof(fSigmabarInnerSector));
84 memset(fSigmabarOuterSector, 0, sizeof(fSigmabarOuterSector));
85}
86
87// --------------------------------------------------------------------------
88//
89// Actual calculation of sigmabar. This is done for each of the six sectors
90// separately due to their possibly different HV behavior. Also inner and
91// outer pixels are treated separately.
92//
93// Preliminary! Works for CT1 test, for real MAGIC crosschecks still have
94// to be done. Also implementation details will be updated, like
95// determination of sector to which a respective pixel belongs
96//
97Float_t MSigmabar::Calc(const MGeomCam &geom, const MPedPhotCam &ped,
98 const MCerPhotEvt &evt)
99{
100 Int_t innerPixels[6];
101 Int_t outerPixels[6];
102 Float_t innerSum[6];
103 Float_t outerSum[6];
104
105 memset(innerPixels, 0, sizeof(innerPixels));
106 memset(outerPixels, 0, sizeof(outerPixels));
107 memset(innerSum, 0, sizeof(innerSum));
108 memset(outerSum, 0, sizeof(outerSum));
109
110 //
111 // sum up sigma/sqrt(area) for each sector,
112 // separately for inner and outer region;
113 //
114 // consider all pixels with Cherenkov photon information
115 // and require "Used"
116 //
117
118 const UInt_t npix = evt.GetNumPixels();
119
120 //*fLog << "MSigmabar : npix = " << npix << endl;
121
122 for (UInt_t i=0; i<npix; i++)
123 {
124 MCerPhotPix &cerpix = evt.operator[](i);
125 if (!cerpix.IsPixelUsed())
126 continue;
127
128 /*
129 if ( cerpix.GetNumPhotons() == 0 )
130 {
131 *fLog << "MSigmabar::Calc(); no.of photons is 0 for used pixel"
132 << endl;
133 continue;
134 }
135 */
136
137 const Int_t idx = cerpix.GetPixId();
138 if (idx == 0)
139 {
140 //*fLog << "MSigmabar : id = 0; pixel '0' is used, ignore it"
141 // << endl;
142 continue;
143 }
144
145 // ratio is the area of pixel 0
146 // divided by the area of the current pixel
147 const Double_t ratio = geom.GetPixRatio(idx);
148
149 const MGeomPix &gpix = geom[idx];
150
151 // This is wrong : rounding has to be done on positive values
152 //Int_t sector = (Int_t)(atan2(gpix.GetY(),gpix.GetX())*6
153 // / (TMath::Pi()*2));
154 //if (sector<0)
155 // sector+=6;
156
157 Float_t sectorf = atan2(gpix.GetY(),gpix.GetX()) * 6.0
158 / (TMath::Pi()*2);
159 if (sectorf < 0.0)
160 sectorf += 6.0;
161 Int_t sector = (Int_t)sectorf;
162
163 // count only those pixels which have a sigma != 0.0
164 const Float_t sigma = ped[idx].GetRms();
165
166 if ( sigma <= 0 )
167 continue;
168
169 //if (area < 1.5)
170 if (ratio > 0.5)
171 {
172 innerPixels[sector]++;
173 innerSum[sector]+= sigma * sqrt(ratio);
174 }
175 else
176 {
177 outerPixels[sector]++;
178 outerSum[sector]+= sigma * sqrt(ratio);
179 }
180 }
181
182 fInnerPixels = 0;
183 fOuterPixels = 0;
184 Double_t fSumInner = 0;
185 Double_t fSumOuter = 0;
186 for (UInt_t i=0; i<6; i++)
187 {
188 fSumInner += innerSum[i];
189 fInnerPixels += innerPixels[i];
190 fSumOuter += outerSum[i];
191 fOuterPixels += outerPixels[i];
192 }
193
194 if (fInnerPixels > 0) fSigmabarInner = fSumInner / fInnerPixels;
195 if (fOuterPixels > 0) fSigmabarOuter = fSumOuter / fOuterPixels;
196
197 //
198 // this is the average sigma/sqrt(area) (over all pixels)
199 //
200 fSigmabar = (fInnerPixels+fOuterPixels)<=0 ? 0:
201 (fSumInner+fSumOuter)/(fInnerPixels+fOuterPixels);
202
203 for (UInt_t i=0; i<6; i++)
204 {
205 const Double_t ip = innerPixels[i];
206 const Double_t op = outerPixels[i];
207 const Double_t iss = innerSum[i];
208 const Double_t oss = outerSum[i];
209
210 const Double_t sum = ip + op;
211 fSigmabarInnerSector[i] = ip <=0 ? 0 : iss/ip;
212 fSigmabarOuterSector[i] = op <=0 ? 0 : oss/op;
213 fSigmabarSector[i] = sum<=0 ? 0 : (iss+oss)/sum;
214 }
215
216 //TString opt = "";
217 //Print(opt);
218
219 return fSigmabar;
220}
221
222// --------------------------------------------------------------------------
223//
224void MSigmabar::Print(Option_t *) const
225{
226 *fLog << all << endl;
227 *fLog << "Total number of inner pixels is " << fInnerPixels << endl;
228 *fLog << "Total number of outer pixels is " << fOuterPixels << endl;
229 *fLog << endl;
230
231 *fLog << "Sigmabar Overall : " << fSigmabar << " ";
232 *fLog << " Sectors: ";
233 for (Int_t i=0;i<6;i++)
234 *fLog << fSigmabarSector[i] << ", ";
235 *fLog << endl;
236
237
238 *fLog << "Sigmabar Inner : " << fSigmabarInner << " ";
239 *fLog << " Sectors: ";
240 for (Int_t i=0;i<6;i++)
241 *fLog << fSigmabarInnerSector[i] << ", ";
242 *fLog << endl;
243
244
245 *fLog << "Sigmabar Outer : " << fSigmabarOuter << " ";
246 *fLog << " Sectors: ";
247 for (Int_t i=0;i<6;i++)
248 *fLog << fSigmabarOuterSector[i] << ", ";
249 *fLog << endl;
250
251}
252
253
254
255
256
257
258
259
260
Note: See TracBrowser for help on using the repository browser.