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@atsro.uni-wuerzburh.de>
|
---|
19 | ! Author(s): Markus Gaug 02/2004 <mailto:markus@ifae.es>
|
---|
20 | !
|
---|
21 | ! Copyright: MAGIC Software Development, 2000-2004
|
---|
22 | !
|
---|
23 | !
|
---|
24 | \* ======================================================================== */
|
---|
25 |
|
---|
26 | /////////////////////////////////////////////////////////////////////////////
|
---|
27 | //
|
---|
28 | // MHPedestalCam
|
---|
29 | //
|
---|
30 | // Fills the extracted pedestals of MExtractedSignalCam into
|
---|
31 | // the MHCalibrationPix-classes MHPedestalPix for every:
|
---|
32 | //
|
---|
33 | // - Pixel, stored in the TObjArray's MHCalibrationCam::fHiGainArray
|
---|
34 | // or MHCalibrationCam::fHiGainArray, respectively
|
---|
35 | //
|
---|
36 | // - Average pixel per AREA index (e.g. inner and outer for the MAGIC camera),
|
---|
37 | // stored in the TObjArray's MHCalibrationCam::fAverageHiGainAreas and
|
---|
38 | // MHCalibrationCam::fAverageHiGainAreas
|
---|
39 | //
|
---|
40 | // - Average pixel per camera SECTOR (e.g. sectors 1-6 for the MAGIC camera),
|
---|
41 | // stored in the TObjArray's MHCalibrationCam::fAverageHiGainSectors
|
---|
42 | // and MHCalibrationCam::fAverageHiGainSectors
|
---|
43 | //
|
---|
44 | // Every pedestal is directly taken from MExtractedSignalCam, filled by the
|
---|
45 | // appropriate extractor.
|
---|
46 | //
|
---|
47 | // The pedestals are filled into a histogram and an array, in order to perform
|
---|
48 | // a Fourier analysis (see MHGausEvents). The signals are moreover averaged on an
|
---|
49 | // event-by-event basis and written into the corresponding average pixels.
|
---|
50 | //
|
---|
51 | // The histograms are fitted to a Gaussian, mean and sigma with its errors
|
---|
52 | // and the fit probability are extracted. If none of these values are NaN's and
|
---|
53 | // if the probability is bigger than MHGausEvents::fProbLimit (default: 0.5%),
|
---|
54 | // the fit is declared valid.
|
---|
55 | // Otherwise, the fit is repeated within ranges of the previous mean
|
---|
56 | // +- MHCalibrationPix::fPickupLimit (default: 5) sigma (see MHCalibrationPix::RepeatFit())
|
---|
57 | // In case this does not make the fit valid, the histogram means and RMS's are
|
---|
58 | // taken directly (see MHCalibrationPix::BypassFit()).
|
---|
59 | //
|
---|
60 | // Outliers of more than MHCalibrationPix::fPickupLimit (default: 5) sigmas
|
---|
61 | // from the mean are counted as Pickup events (stored in MHCalibrationPix::fPickup)
|
---|
62 | //
|
---|
63 | // The number of summed FADC slices is taken to re-normalize
|
---|
64 | // the result to a pedestal per pixel with the formulas (see MHPedestalPix::Renorm()):
|
---|
65 | // - Mean Pedestal / slice = Mean Pedestal / Number slices
|
---|
66 | // - Mean Pedestal Error / slice = Mean Pedestal Error / Number slices
|
---|
67 | // - Sigma Pedestal / slice = Sigma Pedestal / Sqrt (Number slices)
|
---|
68 | // - Sigma Pedestal Error / slice = Sigma Pedestal Error / Sqrt (Number slices)
|
---|
69 | //
|
---|
70 | // The class also fills arrays with the signal vs. event number, creates a fourier
|
---|
71 | // spectrum (see MHGausEvents::CreateFourierSpectrum()) and investigates if the
|
---|
72 | // projected fourier components follow an exponential distribution.
|
---|
73 | //
|
---|
74 | // This same procedure is performed for the average pixels.
|
---|
75 | //
|
---|
76 | // The following results are written into MPedestalCam:
|
---|
77 | //
|
---|
78 | // - MCalibrationPix::SetMean()
|
---|
79 | // - MCalibrationPix::SetMeanErr()
|
---|
80 | // - MCalibrationPix::SetSigma()
|
---|
81 | // - MCalibrationPix::SetSigmaErr()
|
---|
82 | // - MCalibrationPix::SetProb()
|
---|
83 | // - MCalibrationPix::SetNumPickup()
|
---|
84 | //
|
---|
85 | // For all averaged areas, the fitted sigma is multiplied with the square root of
|
---|
86 | // the number involved pixels in order to be able to compare it to the average of
|
---|
87 | // sigmas in the camera.
|
---|
88 | //
|
---|
89 | /////////////////////////////////////////////////////////////////////////////
|
---|
90 | #include "MHPedestalCam.h"
|
---|
91 | #include "MHPedestalPix.h"
|
---|
92 |
|
---|
93 | #include "MBadPixelsCam.h"
|
---|
94 | #include "MBadPixelsPix.h"
|
---|
95 |
|
---|
96 | #include "MLog.h"
|
---|
97 | #include "MLogManip.h"
|
---|
98 |
|
---|
99 | #include "MParList.h"
|
---|
100 |
|
---|
101 | #include "MExtractedSignalCam.h"
|
---|
102 | #include "MExtractedSignalPix.h"
|
---|
103 |
|
---|
104 | #include "MPedestalCam.h"
|
---|
105 | #include "MPedestalPix.h"
|
---|
106 | #include "MCalibrationPix.h"
|
---|
107 |
|
---|
108 | #include "MGeomCam.h"
|
---|
109 | #include "MGeomPix.h"
|
---|
110 |
|
---|
111 | #include "MCalibrationPedCam.h"
|
---|
112 |
|
---|
113 | #include "TH1.h"
|
---|
114 | #include "TH1F.h"
|
---|
115 |
|
---|
116 | #include <TOrdCollection.h>
|
---|
117 |
|
---|
118 | ClassImp(MHPedestalCam);
|
---|
119 |
|
---|
120 | using namespace std;
|
---|
121 |
|
---|
122 | const Int_t MHPedestalCam::fgNbins = 100;
|
---|
123 | const Axis_t MHPedestalCam::fgFirst = -50.;
|
---|
124 | const Axis_t MHPedestalCam::fgLast = 50.;
|
---|
125 | const TString MHPedestalCam::gsHistName = "Pedestal";
|
---|
126 | const TString MHPedestalCam::gsHistTitle = "Pedestal";
|
---|
127 | const TString MHPedestalCam::gsHistXTitle = "Charge [FADC slices]";
|
---|
128 | const TString MHPedestalCam::gsHistYTitle = "Nr. events";
|
---|
129 | // --------------------------------------------------------------------------
|
---|
130 | //
|
---|
131 | // Default constructor.
|
---|
132 | //
|
---|
133 | // Initializes:
|
---|
134 | // - fExtractHiGainSlices to 0.
|
---|
135 | // - the event frequency to 1200 Hz.
|
---|
136 | // - the fRenorm flag to kFALSE
|
---|
137 | //
|
---|
138 | // - fNbins to fgNbins
|
---|
139 | // - fFirst to fgFirst
|
---|
140 | // - fLast to fgLast
|
---|
141 | //
|
---|
142 | // - fHistName to gsHistName
|
---|
143 | // - fHistTitle to gsHistTitle
|
---|
144 | // - fHistXTitle to gsHistXTitle
|
---|
145 | // - fHistYTitle to gsHistYTitle
|
---|
146 | //
|
---|
147 | MHPedestalCam::MHPedestalCam(const char *name, const char *title)
|
---|
148 | : fNumEvents(0), fExtractHiGainSlices(0.)
|
---|
149 | {
|
---|
150 |
|
---|
151 | fName = name ? name : "MHPedestalCam";
|
---|
152 | fTitle = title ? title : "Class to fill the pedestal histograms";
|
---|
153 |
|
---|
154 | SetPulserFrequency(1200);
|
---|
155 | SetRenorm(kTRUE);
|
---|
156 | SetLoGain(kFALSE);
|
---|
157 |
|
---|
158 | SetNbins(fgNbins);
|
---|
159 | SetFirst(fgFirst);
|
---|
160 | SetLast (fgLast );
|
---|
161 |
|
---|
162 | SetHistName (gsHistName .Data());
|
---|
163 | SetHistTitle (gsHistTitle .Data());
|
---|
164 | SetHistXTitle(gsHistXTitle.Data());
|
---|
165 | SetHistYTitle(gsHistYTitle.Data());
|
---|
166 | }
|
---|
167 |
|
---|
168 | // --------------------------------------------------------------------------
|
---|
169 | //
|
---|
170 | // Calls MHCalibrationCam::ResetHists()
|
---|
171 | //
|
---|
172 | // Resets:
|
---|
173 | // - fSum
|
---|
174 | // - fSum2
|
---|
175 | // - fAreaSum
|
---|
176 | // - fAreaSum2
|
---|
177 | // - fAreaVal
|
---|
178 | // - fSectorSum
|
---|
179 | // - fSectorSum2
|
---|
180 | // - fSectorVal
|
---|
181 | //
|
---|
182 | void MHPedestalCam::ResetHists()
|
---|
183 | {
|
---|
184 |
|
---|
185 | MHCalibrationCam::ResetHists();
|
---|
186 |
|
---|
187 | fNumEvents = 0;
|
---|
188 |
|
---|
189 | // If the size is yet set, set the size
|
---|
190 | if (fSum.GetSize()>0)
|
---|
191 | {
|
---|
192 | // Reset contents of arrays.
|
---|
193 | fSum.Reset();
|
---|
194 | fSum2.Reset();
|
---|
195 |
|
---|
196 | fAreaSum. Reset();
|
---|
197 | fAreaSum2.Reset();
|
---|
198 | fAreaVal.Reset();
|
---|
199 |
|
---|
200 | fSectorSum. Reset();
|
---|
201 | fSectorSum2.Reset();
|
---|
202 | fSectorVal.Reset();
|
---|
203 | }
|
---|
204 | }
|
---|
205 |
|
---|
206 | // --------------------------------------------------------------------------
|
---|
207 | //
|
---|
208 | // Creates new MHCalibrationChargeCam only with the averaged areas:
|
---|
209 | // the rest has to be retrieved directly, e.g. via:
|
---|
210 | // MHPedestalCam *cam = MParList::FindObject("MHPedestalCam");
|
---|
211 | // - cam->GetAverageSector(5).DrawClone();
|
---|
212 | // - (*cam)[100].DrawClone()
|
---|
213 | //
|
---|
214 | TObject *MHPedestalCam::Clone(const char *) const
|
---|
215 | {
|
---|
216 |
|
---|
217 | MHPedestalCam *cam = new MHPedestalCam();
|
---|
218 |
|
---|
219 | //
|
---|
220 | // Copy the data members
|
---|
221 | //
|
---|
222 | cam->fRunNumbers = fRunNumbers;
|
---|
223 | cam->fPulserFrequency = fPulserFrequency;
|
---|
224 | cam->fFlags = fFlags;
|
---|
225 | cam->fNbins = fNbins;
|
---|
226 | cam->fFirst = fFirst;
|
---|
227 | cam->fLast = fLast;
|
---|
228 |
|
---|
229 | if (!IsAverageing())
|
---|
230 | return cam;
|
---|
231 |
|
---|
232 | const Int_t navhi = fAverageHiGainAreas->GetSize();
|
---|
233 |
|
---|
234 | for (int i=0; i<navhi; i++)
|
---|
235 | cam->fAverageHiGainAreas->AddAt(GetAverageHiGainArea(i).Clone(),i);
|
---|
236 |
|
---|
237 | return cam;
|
---|
238 | }
|
---|
239 |
|
---|
240 | // --------------------------------------------------------------------------
|
---|
241 | //
|
---|
242 | // Searches pointer to:
|
---|
243 | // - MPedestalCam
|
---|
244 | // - MExtractedSignalCam
|
---|
245 | //
|
---|
246 | // Searches or creates:
|
---|
247 | // - MCalibrationPedCam
|
---|
248 | //
|
---|
249 | // Retrieves from MExtractedSignalCam:
|
---|
250 | // - number of used High Gain FADC slices (MExtractedSignalCam::GetNumUsedHiGainFADCSlices())
|
---|
251 | //
|
---|
252 | // Initializes, if empty to MGeomCam::GetNumPixels():
|
---|
253 | // - MHCalibrationCam::fHiGainArray
|
---|
254 | //
|
---|
255 | // Initializes, if empty to MGeomCam::GetNumAreas() for:
|
---|
256 | // - MHCalibrationCam::fAverageHiGainAreas
|
---|
257 | //
|
---|
258 | // Initializes, if empty to MGeomCam::GetNumSectors() for:
|
---|
259 | // - MHCalibrationCam::fAverageHiGainSectors
|
---|
260 | //
|
---|
261 | // Calls MHCalibrationCam::InitPedHists() for every entry in:
|
---|
262 | // - MHCalibrationCam::fHiGainArray
|
---|
263 | // - MHCalibrationCam::fAverageHiGainAreas
|
---|
264 | // - MHCalibrationCam::fAverageHiGainSectors
|
---|
265 | //
|
---|
266 | // Sets Titles and Names for the Histograms
|
---|
267 | // - MHCalibrationCam::fAverageHiGainAreas
|
---|
268 | //
|
---|
269 | Bool_t MHPedestalCam::ReInitHists(MParList *pList)
|
---|
270 | {
|
---|
271 |
|
---|
272 | MExtractedSignalCam *signal = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
|
---|
273 | if (!signal)
|
---|
274 | {
|
---|
275 | *fLog << err << "Cannot find MExtractedSignalCam... abort." << endl;
|
---|
276 | return kFALSE;
|
---|
277 | }
|
---|
278 |
|
---|
279 |
|
---|
280 | fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
|
---|
281 |
|
---|
282 | if (!fPedestals)
|
---|
283 | *fLog << warn << "Cannot find MPedestalCam ..." << endl;
|
---|
284 |
|
---|
285 | const Int_t npixels = fGeom->GetNumPixels();
|
---|
286 | const Int_t nsectors = fGeom->GetNumSectors();
|
---|
287 | const Int_t nareas = fGeom->GetNumAreas();
|
---|
288 |
|
---|
289 | fCam = (MCalibrationCam*)pList->FindObject("MCalibrationPedCam");
|
---|
290 | if (!fCam)
|
---|
291 | {
|
---|
292 | fCam = (MCalibrationCam*)pList->FindCreateObj(AddSerialNumber("MCalibrationPedCam"));
|
---|
293 | if (!fCam)
|
---|
294 | return kFALSE;
|
---|
295 |
|
---|
296 | fCam->Init(*fGeom);
|
---|
297 | }
|
---|
298 |
|
---|
299 | Float_t sliceshi = signal->GetNumUsedHiGainFADCSlices();
|
---|
300 |
|
---|
301 | if (sliceshi == 0.)
|
---|
302 | {
|
---|
303 | *fLog << err << "Number of used signal slices in MExtractedSignalCam is zero ... abort."
|
---|
304 | << endl;
|
---|
305 | return kFALSE;
|
---|
306 | }
|
---|
307 |
|
---|
308 | if (fExtractHiGainSlices != 0. && sliceshi != fExtractHiGainSlices )
|
---|
309 | {
|
---|
310 | *fLog << err << "Number of used High Gain signal slices changed in MExtractedSignalCam ... abort."
|
---|
311 | << endl;
|
---|
312 | return kFALSE;
|
---|
313 | }
|
---|
314 |
|
---|
315 | fExtractHiGainSlices = sliceshi;
|
---|
316 |
|
---|
317 | *fLog << endl;
|
---|
318 | *fLog << inf << GetDescriptor()
|
---|
319 | << ": Number of found High Gain signal slices: " << fExtractHiGainSlices << endl;
|
---|
320 |
|
---|
321 | InitHiGainArrays(npixels,nareas,nsectors);
|
---|
322 |
|
---|
323 | if (fSum.GetSize()==0)
|
---|
324 | {
|
---|
325 | fSum. Set(npixels);
|
---|
326 | fSum2. Set(npixels);
|
---|
327 | fAreaSum. Set(nareas);
|
---|
328 | fAreaSum2. Set(nareas);
|
---|
329 | fAreaVal. Set(nareas);
|
---|
330 | fSectorSum. Set(nsectors);
|
---|
331 | fSectorSum2.Set(nsectors);
|
---|
332 | fSectorVal.Set(nsectors);
|
---|
333 | }
|
---|
334 |
|
---|
335 | fNumEvents = 0;
|
---|
336 |
|
---|
337 | return kTRUE;
|
---|
338 | }
|
---|
339 |
|
---|
340 |
|
---|
341 | // -------------------------------------------------------------------------------
|
---|
342 | //
|
---|
343 | // Retrieves pointer to MExtractedSignalCam:
|
---|
344 | //
|
---|
345 | // Retrieves from MGeomCam:
|
---|
346 | // - number of pixels
|
---|
347 | // - number of pixel areas
|
---|
348 | // - number of sectors
|
---|
349 | //
|
---|
350 | // Fills HiGain histograms (MHGausEvents::FillHistAndArray()), respectively
|
---|
351 | // with the signals MExtractedSignalCam::GetExtractedSignalHiGain().
|
---|
352 | //
|
---|
353 | Bool_t MHPedestalCam::FillHists(const MParContainer *par, const Stat_t w)
|
---|
354 | {
|
---|
355 |
|
---|
356 | MExtractedSignalCam *signal = (MExtractedSignalCam*)par;
|
---|
357 | if (!signal)
|
---|
358 | {
|
---|
359 | *fLog << err << "No argument in MExtractedSignalCam::Fill... abort." << endl;
|
---|
360 | return kFALSE;
|
---|
361 | }
|
---|
362 |
|
---|
363 | const UInt_t npixels = fGeom->GetNumPixels();
|
---|
364 | const UInt_t nareas = fGeom->GetNumAreas();
|
---|
365 | const UInt_t nsectors = fGeom->GetNumSectors();
|
---|
366 |
|
---|
367 | MArrayF sumareahi(nareas);
|
---|
368 | MArrayF sumsectorhi(nsectors);
|
---|
369 | MArrayI numareahi(nareas);
|
---|
370 | MArrayI numsectorhi(nsectors);
|
---|
371 |
|
---|
372 | for (UInt_t i=0; i<npixels; i++)
|
---|
373 | {
|
---|
374 |
|
---|
375 | MHCalibrationPix &histhi = (*this)[i];
|
---|
376 |
|
---|
377 | if (histhi.IsExcluded())
|
---|
378 | continue;
|
---|
379 |
|
---|
380 | const MExtractedSignalPix &pix = (*signal)[i];
|
---|
381 |
|
---|
382 | const Float_t pedhi = pix.GetExtractedSignalHiGain();
|
---|
383 |
|
---|
384 | const Int_t aidx = (*fGeom)[i].GetAidx();
|
---|
385 | const Int_t sector = (*fGeom)[i].GetSector();
|
---|
386 |
|
---|
387 | histhi.FillHistAndArray(pedhi) ;
|
---|
388 |
|
---|
389 | fSum[i] += pedhi;
|
---|
390 | fAreaSum[aidx] += pedhi;
|
---|
391 | fSectorSum[sector] += pedhi;
|
---|
392 |
|
---|
393 | const Float_t sqrsum = pedhi*pedhi;
|
---|
394 | fSum2[i] += sqrsum;
|
---|
395 | fAreaSum2[aidx] += sqrsum;
|
---|
396 | fSectorSum2[sector] += sqrsum;
|
---|
397 |
|
---|
398 | sumareahi [aidx] += pedhi;
|
---|
399 | numareahi [aidx] ++;
|
---|
400 | sumsectorhi[sector] += pedhi;
|
---|
401 | numsectorhi[sector] ++;
|
---|
402 |
|
---|
403 | }
|
---|
404 |
|
---|
405 | for (UInt_t j=0; j<nareas; j++)
|
---|
406 | {
|
---|
407 | MHCalibrationPix &histhi = GetAverageHiGainArea(j);
|
---|
408 | histhi.FillHistAndArray(numareahi[j] == 0 ? 0. : sumareahi[j]/numareahi[j]);
|
---|
409 | }
|
---|
410 |
|
---|
411 | for (UInt_t j=0; j<nsectors; j++)
|
---|
412 | {
|
---|
413 | MHCalibrationPix &histhi = GetAverageHiGainSector(j);
|
---|
414 | histhi.FillHistAndArray(numsectorhi[j] == 0 ? 0. : sumsectorhi[j]/numsectorhi[j]);
|
---|
415 |
|
---|
416 | }
|
---|
417 |
|
---|
418 | fNumEvents++;
|
---|
419 |
|
---|
420 | return kTRUE;
|
---|
421 | }
|
---|
422 |
|
---|
423 | // --------------------------------------------------------------------------
|
---|
424 | //
|
---|
425 | // Calls:
|
---|
426 | // - MHCalibrationCam::FitHiGainArrays() with Bad Pixels flags 0
|
---|
427 | //
|
---|
428 | Bool_t MHPedestalCam::FinalizeHists()
|
---|
429 | {
|
---|
430 |
|
---|
431 | if (fNumEvents <= 1)
|
---|
432 | {
|
---|
433 | *fLog << err << GetDescriptor()
|
---|
434 | << ": Number of processed events smaller or equal to 1" << endl;
|
---|
435 | return kFALSE;
|
---|
436 | }
|
---|
437 |
|
---|
438 | FitHists();
|
---|
439 |
|
---|
440 | if (fRenorm)
|
---|
441 | RenormResults();
|
---|
442 |
|
---|
443 | return kTRUE;
|
---|
444 |
|
---|
445 | }
|
---|
446 |
|
---|
447 | void MHPedestalCam::FitHists()
|
---|
448 | {
|
---|
449 |
|
---|
450 | for (Int_t i=0; i<fHiGainArray->GetSize(); i++)
|
---|
451 | {
|
---|
452 |
|
---|
453 | MHCalibrationPix &hist = (*this)[i];
|
---|
454 | MCalibrationPix &pix = (*fCam)[i];
|
---|
455 |
|
---|
456 | //
|
---|
457 | // 1) Store calculated means and variances in the low-gain slices
|
---|
458 | //
|
---|
459 | pix.SetLoGainMean ( fSum[i] / fNumEvents );
|
---|
460 | const Double_t diff = fSum2[i] - fSum[i]*fSum[i]/fNumEvents;
|
---|
461 | pix.SetLoGainSigma ( diff > 0. ? TMath::Sqrt( diff / (fNumEvents-1) ) : 0.);
|
---|
462 | pix.SetLoGainMeanVar ( pix.GetLoGainSigma() * pix.GetLoGainSigma() / fNumEvents );
|
---|
463 | pix.SetLoGainSigmaVar( pix.GetLoGainMeanVar() / 4. );
|
---|
464 |
|
---|
465 | if (hist.IsEmpty() || hist.IsOnlyOverflow() || hist.IsOnlyUnderflow())
|
---|
466 | continue;
|
---|
467 |
|
---|
468 | //
|
---|
469 | // 2) Fit the Hi Gain histograms with a Gaussian
|
---|
470 | //
|
---|
471 | hist.FitGaus();
|
---|
472 | //
|
---|
473 | // 3) Take histogram means and RMS
|
---|
474 | //
|
---|
475 | hist.BypassFit();
|
---|
476 | //
|
---|
477 | // 4) Check for oscillations
|
---|
478 | //
|
---|
479 | if (IsOscillations())
|
---|
480 | hist.CreateFourierSpectrum();
|
---|
481 | //
|
---|
482 | // 5) Retrieve the results and store them in this class
|
---|
483 | //
|
---|
484 | pix.SetHiGainMean ( hist.GetMean() );
|
---|
485 | pix.SetHiGainMeanVar ( hist.GetMeanErr() * hist.GetMeanErr() );
|
---|
486 | pix.SetHiGainRms ( hist.GetHistRms() );
|
---|
487 | pix.SetHiGainSigma ( hist.GetSigma() );
|
---|
488 | pix.SetHiGainSigmaVar ( hist.GetSigmaErr()* hist.GetSigmaErr() );
|
---|
489 | pix.SetHiGainProb ( hist.GetProb() );
|
---|
490 | pix.SetHiGainNumBlackout( hist.GetBlackout() );
|
---|
491 | pix.SetHiGainNumPickup ( hist.GetPickup() );
|
---|
492 | }
|
---|
493 |
|
---|
494 | if (!IsAverageing())
|
---|
495 | return;
|
---|
496 |
|
---|
497 | for (Int_t j=0; j<fAverageHiGainAreas->GetSize(); j++)
|
---|
498 | {
|
---|
499 |
|
---|
500 | MHCalibrationPix &hist = GetAverageHiGainArea(j);
|
---|
501 | MCalibrationPix &pix = fCam->GetAverageArea(j);
|
---|
502 |
|
---|
503 | if (hist.IsEmpty() || hist.IsOnlyOverflow() || hist.IsOnlyUnderflow())
|
---|
504 | continue;
|
---|
505 |
|
---|
506 | //
|
---|
507 | // 2) Fit the Hi Gain histograms with a Gaussian
|
---|
508 | //
|
---|
509 | hist.FitGaus();
|
---|
510 | //
|
---|
511 | // 3) Take histogram means and RMS
|
---|
512 | //
|
---|
513 | hist.BypassFit();
|
---|
514 | //
|
---|
515 | // 4) Check for oscillations
|
---|
516 | //
|
---|
517 | if (IsOscillations())
|
---|
518 | hist.CreateFourierSpectrum();
|
---|
519 | //
|
---|
520 | // 5) Retrieve the results and store them in this class
|
---|
521 | //
|
---|
522 | pix.SetHiGainMean ( hist.GetMean() );
|
---|
523 | pix.SetHiGainMeanVar ( hist.GetMeanErr() * hist.GetMeanErr() );
|
---|
524 | pix.SetHiGainRms ( hist.GetHistRms() );
|
---|
525 | pix.SetHiGainSigma ( hist.GetSigma() );
|
---|
526 | pix.SetHiGainSigmaVar ( hist.GetSigmaErr()* hist.GetSigmaErr() );
|
---|
527 | pix.SetHiGainProb ( hist.GetProb() );
|
---|
528 | pix.SetHiGainNumBlackout( hist.GetBlackout() );
|
---|
529 | pix.SetHiGainNumPickup ( hist.GetPickup() );
|
---|
530 | //
|
---|
531 | // 6) Store calculated means and variances in the low-gain slices
|
---|
532 | //
|
---|
533 | const ULong_t aevts = fNumEvents * fAreaVal[j];
|
---|
534 | if (aevts <= 1)
|
---|
535 | continue;
|
---|
536 |
|
---|
537 | pix.SetLoGainMean ( fAreaSum[j] / aevts );
|
---|
538 | const Double_t diff = fAreaSum2[j] - fAreaSum[j]*fAreaSum[j]/aevts ;
|
---|
539 | pix.SetLoGainSigma( diff > 0. ? TMath::Sqrt( diff / (aevts-1) ) : 0.);
|
---|
540 | }
|
---|
541 |
|
---|
542 | for (Int_t j=0; j<fAverageHiGainSectors->GetSize(); j++)
|
---|
543 | {
|
---|
544 | MHCalibrationPix &hist = GetAverageHiGainSector(j);
|
---|
545 | MCalibrationPix &pix = fCam->GetAverageSector(j);
|
---|
546 |
|
---|
547 | if (hist.IsEmpty() || hist.IsOnlyOverflow() || hist.IsOnlyUnderflow())
|
---|
548 | continue;
|
---|
549 |
|
---|
550 | //
|
---|
551 | // 2) Fit the Hi Gain histograms with a Gaussian
|
---|
552 | //
|
---|
553 | hist.FitGaus();
|
---|
554 | //
|
---|
555 | // 3) Take histogram means and RMS
|
---|
556 | //
|
---|
557 | hist.BypassFit();
|
---|
558 | //
|
---|
559 | // 4) Check for oscillations
|
---|
560 | //
|
---|
561 | if (IsOscillations())
|
---|
562 | hist.CreateFourierSpectrum();
|
---|
563 | //
|
---|
564 | // 5) Retrieve the results and store them in this class
|
---|
565 | //
|
---|
566 | pix.SetHiGainMean ( hist.GetMean() );
|
---|
567 | pix.SetHiGainMeanVar ( hist.GetMeanErr() * hist.GetMeanErr() );
|
---|
568 | pix.SetHiGainRms ( hist.GetHistRms() );
|
---|
569 | pix.SetHiGainSigma ( hist.GetSigma() );
|
---|
570 | pix.SetHiGainSigmaVar ( hist.GetSigmaErr()* hist.GetSigmaErr() );
|
---|
571 | pix.SetHiGainProb ( hist.GetProb() );
|
---|
572 | pix.SetHiGainNumBlackout( hist.GetBlackout() );
|
---|
573 | pix.SetHiGainNumPickup ( hist.GetPickup() );
|
---|
574 | //
|
---|
575 | // 6) Store calculated means and variances in the low-gain slices
|
---|
576 | //
|
---|
577 | const ULong_t sevts = fNumEvents * fSectorVal[j];
|
---|
578 | if (sevts <= 1)
|
---|
579 | continue;
|
---|
580 |
|
---|
581 | pix.SetLoGainMean ( fSectorSum[j] / sevts );
|
---|
582 | const Double_t diff = fSectorSum2[j] - fSectorSum[j]*fSectorSum[j]/sevts ;
|
---|
583 | pix.SetLoGainSigma( diff > 0. ? TMath::Sqrt( diff / (sevts-1) ) : 0.);
|
---|
584 | }
|
---|
585 | }
|
---|
586 |
|
---|
587 | // --------------------------------------------------------------------------
|
---|
588 | //
|
---|
589 | // Renormalizes the pedestal fit results by the following formulae:
|
---|
590 | //
|
---|
591 | // - Mean Pedestal / slice -> Mean Pedestal / Number slices
|
---|
592 | // - Mean Pedestal Error / slice -> Mean Pedestal Error / Number slices
|
---|
593 | // - Sigma Pedestal / slice -> Sigma Pedestal / Sqrt (Number slices)
|
---|
594 | // - Sigma Pedestal Error / slice -> Sigma Pedestal Error / Sqrt (Number slices)
|
---|
595 | //
|
---|
596 | void MHPedestalCam::RenormResults()
|
---|
597 | {
|
---|
598 |
|
---|
599 | //
|
---|
600 | // One never knows...
|
---|
601 | //
|
---|
602 | if (fExtractHiGainSlices <= 0)
|
---|
603 | return;
|
---|
604 |
|
---|
605 | const Float_t sqslices = TMath::Sqrt(fExtractHiGainSlices);
|
---|
606 |
|
---|
607 | for (Int_t i=0; i<fHiGainArray->GetSize(); i++)
|
---|
608 | {
|
---|
609 |
|
---|
610 | MCalibrationPix &pix = (*fCam)[i];
|
---|
611 | pix.SetHiGainMean ( pix.GetHiGainMean() / fExtractHiGainSlices );
|
---|
612 | pix.SetLoGainMean ( pix.GetLoGainMean() / fExtractHiGainSlices );
|
---|
613 | //
|
---|
614 | // Mean error goes with PedestalRMS/Sqrt(entries) -> scale with sqrt(slices)
|
---|
615 | //
|
---|
616 | pix.SetHiGainMeanVar ( pix.GetHiGainMeanVar() / fExtractHiGainSlices );
|
---|
617 | pix.SetLoGainMeanVar ( pix.GetHiGainMeanVar() / fExtractHiGainSlices );
|
---|
618 | //
|
---|
619 | // Sigma goes like PedestalRMS -> scale with sqrt(slices)
|
---|
620 | //
|
---|
621 | pix.SetHiGainSigma ( pix.GetHiGainSigma() / sqslices );
|
---|
622 | pix.SetLoGainSigma ( pix.GetLoGainSigma() / sqslices );
|
---|
623 | //
|
---|
624 | // Sigma error goes like PedestalRMS/2.(entries) -> scale with sqrt(slices)
|
---|
625 | //
|
---|
626 | pix.SetHiGainSigmaVar ( pix.GetHiGainSigmaVar() / fExtractHiGainSlices );
|
---|
627 | pix.SetLoGainSigmaVar ( pix.GetLoGainSigmaVar() / fExtractHiGainSlices );
|
---|
628 | }
|
---|
629 |
|
---|
630 | if (!IsAverageing())
|
---|
631 | return;
|
---|
632 |
|
---|
633 | for (Int_t j=0; j<fAverageHiGainAreas->GetSize(); j++)
|
---|
634 | {
|
---|
635 |
|
---|
636 | MCalibrationPix &pix = fCam->GetAverageArea(j);
|
---|
637 | pix.SetHiGainMean ( pix.GetHiGainMean() / fExtractHiGainSlices );
|
---|
638 | pix.SetLoGainMean ( pix.GetLoGainMean() / fExtractHiGainSlices );
|
---|
639 | pix.SetHiGainMeanVar ( pix.GetHiGainMeanVar() / fExtractHiGainSlices );
|
---|
640 | pix.SetHiGainSigma ( pix.GetHiGainSigma() / sqslices );
|
---|
641 | pix.SetLoGainSigma ( pix.GetLoGainSigma() / sqslices );
|
---|
642 | pix.SetHiGainSigmaVar ( pix.GetHiGainSigmaVar() / fExtractHiGainSlices );
|
---|
643 | }
|
---|
644 |
|
---|
645 | for (Int_t j=0; j<fAverageHiGainSectors->GetSize(); j++)
|
---|
646 | {
|
---|
647 | MCalibrationPix &pix = fCam->GetAverageSector(j);
|
---|
648 | pix.SetHiGainMean ( pix.GetHiGainMean() / fExtractHiGainSlices );
|
---|
649 | pix.SetLoGainMean ( pix.GetLoGainMean() / fExtractHiGainSlices );
|
---|
650 | pix.SetHiGainMeanVar ( pix.GetHiGainMeanVar() / fExtractHiGainSlices );
|
---|
651 | pix.SetHiGainSigma ( pix.GetHiGainSigma() / sqslices );
|
---|
652 | pix.SetLoGainSigma ( pix.GetLoGainSigma() / sqslices );
|
---|
653 | pix.SetHiGainSigmaVar ( pix.GetHiGainSigmaVar() / fExtractHiGainSlices );
|
---|
654 | }
|
---|
655 |
|
---|
656 | }
|
---|
657 |
|
---|
658 |
|
---|
659 | // ------------------------------------------------------------------
|
---|
660 | //
|
---|
661 | // The types are as follows:
|
---|
662 | //
|
---|
663 | // Fitted values:
|
---|
664 | // ==============
|
---|
665 | //
|
---|
666 | // 0: Fitted Charge
|
---|
667 | // 1: Error of fitted Charge
|
---|
668 | // 2: Sigma of fitted Charge
|
---|
669 | // 3: Error of Sigma of fitted Charge
|
---|
670 | //
|
---|
671 | //
|
---|
672 | // Useful variables derived from the fit results:
|
---|
673 | // =============================================
|
---|
674 | //
|
---|
675 | // 4: Returned probability of Gauss fit to Charge distribution
|
---|
676 | // 5: Relative differenc of calculated pedestal (per slice) and fitted (per slice)
|
---|
677 | // 6: Error of the Relative differenc of calculated pedestal (per slice) and fitted (per slice)
|
---|
678 | // 7: Relative difference of the error of the mean pedestal (per slice) - calculated and fitted
|
---|
679 | // 8: Relative differenc of calculated pedestal RMS (per slice) and fitted sigma (per slice)
|
---|
680 | // 9: Error of Relative differenc of calculated pedestal RMS (per slice) and fitted sigma (per slice)
|
---|
681 | // 10: Relative difference of the error of the pedestal RMS (per slice) - calculated and fitted
|
---|
682 | //
|
---|
683 | // Localized defects:
|
---|
684 | // ==================
|
---|
685 | //
|
---|
686 | // 11: Gaus fit not OK
|
---|
687 | // 12: Fourier spectrum not OK
|
---|
688 | //
|
---|
689 | Bool_t MHPedestalCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
|
---|
690 | {
|
---|
691 |
|
---|
692 | if (fHiGainArray->GetSize() <= idx)
|
---|
693 | return kFALSE;
|
---|
694 |
|
---|
695 | if ((*this)[idx].IsExcluded())
|
---|
696 | return kFALSE;
|
---|
697 |
|
---|
698 | const Float_t ped = (*fPedestals)[idx].GetPedestal();
|
---|
699 | const Float_t rms = (*fPedestals)[idx].GetPedestalRms();
|
---|
700 |
|
---|
701 | const Float_t entsqr = TMath::Sqrt((Float_t)fPedestals->GetTotalEntries());
|
---|
702 |
|
---|
703 | const Float_t pederr = rms/entsqr;
|
---|
704 | const Float_t rmserr = rms/entsqr/2.;
|
---|
705 |
|
---|
706 | const Float_t mean = (*this)[idx].GetMean();
|
---|
707 | const Float_t meanerr = (*this)[idx].GetMeanErr();
|
---|
708 | const Float_t sigma = (*this)[idx].GetSigma() ;
|
---|
709 | const Float_t sigmaerr = (*this)[idx].GetSigmaErr();
|
---|
710 |
|
---|
711 | switch (type)
|
---|
712 | {
|
---|
713 | case 0:
|
---|
714 | val = mean;
|
---|
715 | break;
|
---|
716 | case 1:
|
---|
717 | val = meanerr;
|
---|
718 | break;
|
---|
719 | case 2:
|
---|
720 | val = sigma;
|
---|
721 | break;
|
---|
722 | case 3:
|
---|
723 | val = sigmaerr;
|
---|
724 | break;
|
---|
725 | case 4:
|
---|
726 | val = (*this)[idx].GetProb();
|
---|
727 | break;
|
---|
728 | case 5:
|
---|
729 | val = 2.*(mean-ped)/(ped+mean);
|
---|
730 | break;
|
---|
731 | case 6:
|
---|
732 | val = TMath::Sqrt((pederr*pederr + meanerr*meanerr) * (ped*ped + mean*mean))
|
---|
733 | *2./(ped+mean)/(ped+mean);
|
---|
734 | break;
|
---|
735 | case 7:
|
---|
736 | val = 2.*(meanerr-pederr)/(pederr + meanerr);
|
---|
737 | break;
|
---|
738 | case 8:
|
---|
739 | val = 2.*(sigma-rms)/(sigma+rms);
|
---|
740 | break;
|
---|
741 | case 9:
|
---|
742 | val = TMath::Sqrt((rmserr*rmserr + sigmaerr*sigmaerr) * (rms*rms + sigma*sigma))
|
---|
743 | *2./(rms+sigma)/(rms+sigma);
|
---|
744 | break;
|
---|
745 | case 10:
|
---|
746 | val = 2.*(sigmaerr - rmserr)/(sigmaerr + rmserr);
|
---|
747 | break;
|
---|
748 | case 11:
|
---|
749 | if (!(*this)[idx].IsGausFitOK())
|
---|
750 | val = 1.;
|
---|
751 | break;
|
---|
752 | case 12:
|
---|
753 | if (!(*this)[idx].IsFourierSpectrumOK())
|
---|
754 | val = 1.;
|
---|
755 | break;
|
---|
756 | default:
|
---|
757 | return kFALSE;
|
---|
758 | }
|
---|
759 | return kTRUE;
|
---|
760 | }
|
---|
761 |
|
---|
762 | // --------------------------------------------------------------------------
|
---|
763 | //
|
---|
764 | // Calls MHGausEvents::DrawClone() for pixel idx
|
---|
765 | //
|
---|
766 | void MHPedestalCam::DrawPixelContent(Int_t idx) const
|
---|
767 | {
|
---|
768 | (*this)[idx].DrawClone();
|
---|
769 | }
|
---|