source: trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeCam.cc@ 4661

Last change on this file since 4661 was 4658, checked in by gaug, 20 years ago
*** empty log message ***
File size: 37.3 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): Markus Gaug 02/2004 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24/////////////////////////////////////////////////////////////////////////////
25//
26// MHCalibrationChargeCam
27//
28// Fills the extracted signals of MExtractedSignalCam into the MHGausEvents-classes
29// MHCalibrationChargeHiGainPix and MHCalibrationChargeLoGainPix for every:
30//
31// - Pixel, stored in the TObjArray's MHCalibrationCam::fHiGainArray and
32// MHCalibrationCam::fLoGainArray
33//
34// - Average pixel per AREA index (e.g. inner and outer for the MAGIC camera),
35// stored in the TObjArray's MHCalibrationCam::fAverageHiGainAreas and
36// MHCalibrationCam::fAverageLoGainAreas
37//
38// - Average pixel per camera SECTOR (e.g. sectors 1-6 for the MAGIC camera),
39// stored in the TObjArray's MHCalibrationCam::fAverageHiGainSectors and
40// MHCalibrationCam::fAverageLoGainSectors
41//
42// Every signal is taken from MExtractedSignalCam and filled into a histogram and
43// an array, in order to perform a Fourier analysis (see MHGausEvents).
44// The signals are moreover averaged on an event-by-event basis and written into
45// the corresponding average pixels.
46//
47// Additionally, the (FADC slice) position of the maximum is stored in an Absolute
48// Arrival Time histogram. This histogram serves for a rough cross-check if the
49// signal does not lie at or outside the edges of the extraction window.
50//
51// The Charge 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// +- MHGausEvents::fPickupLimit (default: 5) sigma (see MHGausEvents::RepeatFit())
57// In case this does not make the fit valid, the histogram means and RMS's are
58// taken directly (see MHGausEvents::BypassFit()) and the following flags are set:
59// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted ) or
60// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainNotFitted ) and
61// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun )
62//
63// Outliers of more than MHGausEvents::fPickupLimit (default: 5) sigmas
64// from the mean are counted as Pickup events (stored in MHGausEvents::fPickup)
65//
66// Unless more than fNumHiGainSaturationLimit (default: 1%) of the overall FADC
67// slices show saturation, the following flag is set:
68// - MCalibrationChargePix::SetHiGainSaturation();
69// In that case, the calibration constants are derived from the low-gain results.
70//
71// If more than fNumLoGainSaturationLimit (default: 1%) of the overall
72// low-gain FADC slices saturate, the following flags are set:
73// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainSaturation ) and
74// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnsuitableRun )
75//
76// The class also fills arrays with the signal vs. event number, creates a fourier
77// spectrum and investigates if the projected fourier components follow an exponential
78// distribution. In case that the probability of the exponential fit is less than
79// MHGausEvents::fProbLimit (default: 0.5%), the following flags are set:
80// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainOscillating ) or
81// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainOscillating ) and
82// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun )
83//
84// This same procedure is performed for the average pixels.
85//
86// The following results are written into MCalibrationChargeCam:
87//
88// - MCalibrationPix::SetHiGainSaturation()
89// - MCalibrationPix::SetHiGainMean()
90// - MCalibrationPix::SetHiGainMeanErr()
91// - MCalibrationPix::SetHiGainSigma()
92// - MCalibrationPix::SetHiGainSigmaErr()
93// - MCalibrationPix::SetHiGainProb()
94// - MCalibrationPix::SetHiGainNumPickup()
95//
96// - MCalibrationPix::SetLoGainMean()
97// - MCalibrationPix::SetLoGainMeanErr()
98// - MCalibrationPix::SetLoGainSigma()
99// - MCalibrationPix::SetLoGainSigmaErr()
100// - MCalibrationPix::SetLoGainProb()
101// - MCalibrationPix::SetLoGainNumPickup()
102//
103// - MCalibrationChargePix::SetAbsTimeMean()
104// - MCalibrationChargePix::SetAbsTimeRms()
105//
106// For all averaged areas, the fitted sigma is multiplied with the square root of
107// the number involved pixels in order to be able to compare it to the average of
108// sigmas in the camera.
109//
110/////////////////////////////////////////////////////////////////////////////
111#include "MHCalibrationChargeCam.h"
112#include "MHCalibrationCam.h"
113
114#include "MLog.h"
115#include "MLogManip.h"
116
117#include "MParList.h"
118
119#include "MHCalibrationChargeHiGainPix.h"
120#include "MHCalibrationChargeLoGainPix.h"
121#include "MHCalibrationChargePix.h"
122
123#include "MCalibrationCam.h"
124
125#include "MCalibrationChargeCam.h"
126#include "MCalibrationChargePix.h"
127
128#include "MGeomCam.h"
129#include "MGeomPix.h"
130
131#include "MHGausEvents.h"
132
133#include "MBadPixelsCam.h"
134#include "MBadPixelsPix.h"
135
136#include "MRawEvtData.h"
137#include "MRawEvtPixelIter.h"
138
139#include "MExtractedSignalCam.h"
140#include "MExtractedSignalPix.h"
141
142#include <TPad.h>
143#include <TVirtualPad.h>
144#include <TCanvas.h>
145#include <TStyle.h>
146#include <TF1.h>
147#include <TH2D.h>
148#include <TLine.h>
149#include <TLatex.h>
150#include <TLegend.h>
151
152ClassImp(MHCalibrationChargeCam);
153
154using namespace std;
155
156const Float_t MHCalibrationChargeCam::fgNumHiGainSaturationLimit = 0.01;
157const Float_t MHCalibrationChargeCam::fgNumLoGainSaturationLimit = 0.005;
158const Float_t MHCalibrationChargeCam::fgTimeLowerLimit = 1.;
159const Float_t MHCalibrationChargeCam::fgTimeUpperLimit = 2.;
160// 1Led Green, 1 LED blue, 5 LEDs blue, 10 LEDs blue, 10 LEDs UV, CT1, 5Leds Green
161const Float_t MHCalibrationChargeCam::gkHiGainInnerRefLines[7] = { 245., 323. , 1065., 1467., 180., 211. , 533.5};
162const Float_t MHCalibrationChargeCam::gkHiGainOuterRefLines[7] = { 217., 307.5, 932. , 1405., 167., 183.5, 405.5};
163const Float_t MHCalibrationChargeCam::gkLoGainInnerRefLines[7] = { 20.8, 28.0 , 121. , 200.2, 16.5, 13.5 , 41.7 };
164const Float_t MHCalibrationChargeCam::gkLoGainOuterRefLines[7] = { 18.9, 26.0 , 108.3, 198. , 14.0, 11. , 42. };
165// --------------------------------------------------------------------------
166//
167// Default Constructor.
168//
169// Sets:
170// - all pointers to NULL
171//
172// Initializes:
173// - fNumHiGainSaturationLimit to fgNumHiGainSaturationLimit
174// - fNumLoGainSaturationLimit to fgNumLoGainSaturationLimit
175// - fTimeLowerLimit to fgTimeLowerLimit
176// - fTimeUpperLimit to fgTimeUpperLimit
177//
178MHCalibrationChargeCam::MHCalibrationChargeCam(const char *name, const char *title)
179 : fRawEvt(NULL)
180{
181 fName = name ? name : "MHCalibrationChargeCam";
182 fTitle = title ? title : "Class to fill the calibration histograms ";
183
184 SetNumHiGainSaturationLimit(fgNumHiGainSaturationLimit);
185 SetNumLoGainSaturationLimit(fgNumLoGainSaturationLimit);
186 SetTimeLowerLimit();
187 SetTimeUpperLimit();
188
189
190}
191
192// --------------------------------------------------------------------------
193//
194// Gets the pointers to:
195// - MRawEvtData
196//
197Bool_t MHCalibrationChargeCam::SetupHists(const MParList *pList)
198{
199
200 fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
201 if (!fRawEvt)
202 {
203 *fLog << err << dbginf << "MRawEvtData not found... aborting." << endl;
204 return kFALSE;
205 }
206
207 return kTRUE;
208}
209
210// --------------------------------------------------------------------------
211//
212// Gets or creates the pointers to:
213// - MExtractedSignalCam
214// - MCalibrationChargeCam
215// - MBadPixelsCam
216//
217// Initializes the number of used FADC slices from MExtractedSignalCam
218// into MCalibrationChargeCam and test for changes in that variable
219//
220// Initializes, if empty to MGeomCam::GetNumPixels():
221// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray
222//
223// Initializes, if empty to MGeomCam::GetNumAreas() for:
224// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
225//
226// Initializes, if empty to MGeomCam::GetNumSectors() for:
227// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
228//
229// Calls MHCalibrationCam::InitHists() for every entry in:
230// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray
231// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
232// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
233//
234// Sets Titles and Names for the Charge Histograms:
235// - MHCalibrationCam::fAverageHiGainAreas
236// - MHCalibrationCam::fAverageHiGainSectors
237//
238// Sets number of bins to MHCalibrationCam::fAverageNbins for:
239// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
240// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
241//
242Bool_t MHCalibrationChargeCam::ReInitHists(MParList *pList)
243{
244
245 MExtractedSignalCam *signal = (MExtractedSignalCam*)pList->FindObject(AddSerialNumber("MExtractedSignalCam"));
246 if (!signal)
247 {
248 *fLog << err << "MExtractedSignalCam not found... abort." << endl;
249 return kFALSE;
250 }
251
252 fCam = (MCalibrationCam*)pList->FindObject(AddSerialNumber("MCalibrationChargeCam"));
253 if (!fCam)
254 {
255 fCam = (MCalibrationCam*)pList->FindCreateObj(AddSerialNumber("MCalibrationChargeCam"));
256 if (!fCam)
257 {
258 gLog << err << "Cannot find nor create MCalibrationChargeCam ... abort." << endl;
259 return kFALSE;
260 }
261 else
262 fCam->Init(*fGeom);
263 }
264
265 fFirstHiGain = signal->GetFirstUsedSliceHiGain();
266 fLastHiGain = signal->GetLastUsedSliceHiGain();
267 fFirstLoGain = signal->GetFirstUsedSliceLoGain();
268 fLastLoGain = signal->GetLastUsedSliceLoGain();
269
270 const Float_t numhigain = signal->GetNumUsedHiGainFADCSlices();
271 const Float_t numlogain = signal->GetNumUsedLoGainFADCSlices();
272
273 if (fCam->GetNumHiGainFADCSlices() == 0.)
274 fCam->SetNumHiGainFADCSlices ( numhigain );
275 else if (fCam->GetNumHiGainFADCSlices() != numhigain)
276 {
277 *fLog << err << GetDescriptor()
278 << ": Number of High Gain FADC extraction slices has changed, abort..." << endl;
279 return kFALSE;
280 }
281
282 if (fCam->GetNumLoGainFADCSlices() == 0.)
283 fCam->SetNumLoGainFADCSlices ( numlogain );
284 else if (fCam->GetNumLoGainFADCSlices() != numlogain)
285 {
286 *fLog << err << GetDescriptor()
287 << ": Number of Low Gain FADC extraction slices has changes, abort..." << endl;
288 return kFALSE;
289 }
290
291 const Int_t npixels = fGeom->GetNumPixels();
292 const Int_t nsectors = fGeom->GetNumSectors();
293 const Int_t nareas = fGeom->GetNumAreas();
294
295 if (fHiGainArray->GetEntries()==0)
296 {
297 fHiGainArray->Expand(npixels);
298 for (Int_t i=0; i<npixels; i++)
299 {
300 (*fHiGainArray)[i] = new MHCalibrationChargeHiGainPix;
301 InitHists((*this)[i],(*fBadPixels)[i],i);
302 }
303 }
304
305 if (fLoGainArray->GetEntries()==0)
306 {
307 fLoGainArray->Expand(npixels);
308
309 for (Int_t i=0; i<npixels; i++)
310 {
311 (*fLoGainArray)[i] = new MHCalibrationChargeLoGainPix;
312 InitHists((*this)(i),(*fBadPixels)[i],i);
313 }
314
315 }
316
317 if (fAverageHiGainAreas->GetEntries()==0)
318 {
319 fAverageHiGainAreas->Expand(nareas);
320
321 for (Int_t j=0; j<nareas; j++)
322 {
323 (*fAverageHiGainAreas)[j] =
324 new MHCalibrationChargeHiGainPix("AverageHiGainArea",
325 "Average HiGain FADC sums area idx ");
326
327 MHCalibrationChargePix &hist = (MHCalibrationChargePix&)GetAverageHiGainArea(j);
328
329 hist.SetNbins(fAverageNbins);
330 hist.GetHAbsTime()->SetTitle("Absolute Arrival Time average HiGain Area Idx ");
331
332 if (fGeom->InheritsFrom("MGeomCamMagic"))
333 {
334 hist.GetHGausHist()->SetTitle(Form("%s%s%s","Signal averaged on event-by-event basis ",
335 j==0 ? "Inner Pixels " : "Outer Pixels ","High Gain Runs: "));
336 hist.InitBins();
337 hist.SetEventFrequency(fPulserFrequency);
338 }
339 else
340 {
341 hist.GetHGausHist()->SetTitle("Signal averaged on event-by-event basis High Gain Area Idx ");
342 InitHists(hist,fCam->GetAverageBadArea(j),j);
343 }
344 }
345 }
346
347
348 if (fAverageLoGainAreas->GetEntries()==0)
349 {
350 fAverageLoGainAreas->Expand(nareas);
351
352 for (Int_t j=0; j<nareas; j++)
353 {
354 (*fAverageLoGainAreas)[j] =
355 new MHCalibrationChargeLoGainPix("AverageLoGainArea",
356 "Average LoGain FADC sums of pixel area idx ");
357
358 MHCalibrationChargePix &hist = (MHCalibrationChargePix&)GetAverageLoGainArea(j);
359
360 hist.SetNbins(fAverageNbins);
361 hist.GetHAbsTime()->SetTitle("Absolute Arrival Time average LoGain Area Idx ");
362
363 if (fGeom->InheritsFrom("MGeomCamMagic"))
364 {
365 hist.GetHGausHist()->SetTitle(Form("%s%s%s","Signal averaged on event-by-event basis ",
366 j==0 ? "Inner Pixels " : "Outer Pixels ","Low Gain Runs: "));
367 hist.InitBins();
368 hist.SetEventFrequency(fPulserFrequency);
369 }
370 else
371 {
372 hist.GetHGausHist()->SetTitle("Signal averaged on event-by-event basis High Gain Area Idx ");
373 InitHists(hist,fCam->GetAverageBadArea(j),j);
374 }
375 }
376 }
377
378 if (fAverageHiGainSectors->GetEntries()==0)
379 {
380 fAverageHiGainSectors->Expand(nsectors);
381
382 for (Int_t j=0; j<nsectors; j++)
383 {
384 (*fAverageHiGainSectors)[j] =
385 new MHCalibrationChargeHiGainPix("AverageHiGainSector",
386 "Average HiGain FADC sums of pixel sector ");
387
388 MHCalibrationChargePix &hist = (MHCalibrationChargePix&)GetAverageHiGainSector(j);
389
390 hist.GetHGausHist()->SetTitle("Summed FADC slices average HiGain Sector ");
391 hist.SetNbins(fAverageNbins);
392 hist.GetHAbsTime()->SetTitle("Absolute Arrival Time average HiGain Sector ");
393
394 InitHists(hist,fCam->GetAverageBadSector(j),j);
395
396 }
397 }
398
399 if (fAverageLoGainSectors->GetEntries()==0)
400 {
401 fAverageLoGainSectors->Expand(nsectors);
402
403 for (Int_t j=0; j<nsectors; j++)
404 {
405 (*fAverageLoGainSectors)[j] =
406 new MHCalibrationChargeLoGainPix("AverageLoGainSector",
407 "Average LoGain FADC sums of pixel sector ");
408
409 MHCalibrationChargePix &hist = (MHCalibrationChargePix&)GetAverageLoGainSector(j);
410
411 hist.GetHGausHist()->SetTitle("Summed FADC slices average LoGain Sector ");
412 hist.SetNbins(fAverageNbins);
413 hist.GetHAbsTime()->SetTitle("Absolute Arrival Time average LoGain Sector ");
414
415 InitHists(hist,fCam->GetAverageBadSector(j),j);
416
417 }
418 }
419
420 return kTRUE;
421}
422
423
424// --------------------------------------------------------------------------
425//
426// Retrieves from MExtractedSignalCam:
427// - first used LoGain FADC slice
428//
429// Retrieves from MGeomCam:
430// - number of pixels
431// - number of pixel areas
432// - number of sectors
433//
434// For all TObjArray's (including the averaged ones), the following steps are performed:
435//
436// 1) Fill Charges histograms (MHGausEvents::FillHistAndArray()) with:
437// - MExtractedSignalPix::GetExtractedSignalHiGain();
438// - MExtractedSignalPix::GetExtractedSignalLoGain();
439//
440// 2) Set number of saturated slices (MHCalibrationChargePix::SetSaturated()) with:
441// - MExtractedSignalPix::GetNumHiGainSaturated();
442// - MExtractedSignalPix::GetNumLoGainSaturated();
443//
444// 3) Fill AbsTime histograms (MHCalibrationChargePix::FillAbsTime()) with:
445// - MRawEvtPixelIter::GetIdxMaxHiGainSample();
446// - MRawEvtPixelIter::GetIdxMaxLoGainSample(first slice);
447//
448Bool_t MHCalibrationChargeCam::FillHists(const MParContainer *par, const Stat_t w)
449{
450
451 MExtractedSignalCam *signal = (MExtractedSignalCam*)par;
452 if (!signal)
453 {
454 *fLog << err << "No argument in MExtractedSignalCam::Fill... abort." << endl;
455 return kFALSE;
456 }
457
458 const UInt_t npixels = fGeom->GetNumPixels();
459 const UInt_t nareas = fGeom->GetNumAreas();
460 const UInt_t nsectors = fGeom->GetNumSectors();
461 const UInt_t lofirst = signal->GetFirstUsedSliceLoGain();
462
463 Float_t sumhiarea [nareas], sumloarea [nareas], timehiarea [nareas], timeloarea [nareas];
464 Float_t sumhisector[nsectors], sumlosector[nsectors], timehisector[nsectors], timelosector[nsectors];
465 Int_t sathiarea [nareas], satloarea [nareas];
466 Int_t sathisector[nsectors], satlosector[nsectors];
467
468 memset(sumhiarea, 0, nareas * sizeof(Float_t));
469 memset(sumloarea, 0, nareas * sizeof(Float_t));
470 memset(timehiarea, 0, nareas * sizeof(Float_t));
471 memset(timeloarea, 0, nareas * sizeof(Float_t));
472 memset(sathiarea, 0, nareas * sizeof(Int_t ));
473 memset(satloarea, 0, nareas * sizeof(Int_t ));
474 memset(sumhisector, 0, nsectors*sizeof(Float_t));
475 memset(sumlosector, 0, nsectors*sizeof(Float_t));
476 memset(timehisector,0, nsectors*sizeof(Float_t));
477 memset(timelosector,0, nsectors*sizeof(Float_t));
478 memset(sathisector, 0, nsectors*sizeof(Int_t ));
479 memset(satlosector, 0, nsectors*sizeof(Int_t ));
480
481 for (UInt_t i=0; i<npixels; i++)
482 {
483
484 MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)(*this)[i];
485 MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)(*this)(i);
486
487 if (histhi.IsExcluded())
488 continue;
489
490 const MExtractedSignalPix &pix = (*signal)[i];
491
492 const Float_t sumhi = pix.GetExtractedSignalHiGain();
493 const Float_t sumlo = pix.GetExtractedSignalLoGain();
494
495 if (!histhi.FillHistAndArray(sumhi))
496 fHiGainOverFlow++;
497 if (!histlo.FillHistAndArray(sumlo))
498 fLoGainOverFlow++;
499
500 const Int_t sathi = (Int_t)pix.GetNumHiGainSaturated();
501 const Int_t satlo = (Int_t)pix.GetNumLoGainSaturated();
502
503 histhi.SetSaturated(sathi);
504 histlo.SetSaturated(satlo);
505
506 const Int_t aidx = (*fGeom)[i].GetAidx();
507 const Int_t sector = (*fGeom)[i].GetSector();
508
509 sumhiarea[aidx] += sumhi;
510 sumloarea[aidx] += sumlo;
511 sathiarea[aidx] += sathi;
512 satloarea[aidx] += satlo;
513
514 sumhisector[sector] += sumhi;
515 sumlosector[sector] += sumlo;
516 sathisector[sector] += sathi;
517 satlosector[sector] += satlo;
518 }
519
520 MRawEvtPixelIter pixel(fRawEvt);
521 while (pixel.Next())
522 {
523
524 const UInt_t pixid = pixel.GetPixelId();
525
526 MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)(*this)[pixid];
527 MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)(*this)(pixid);
528
529 if (histhi.IsExcluded())
530 continue;
531
532 const Float_t timehi = (Float_t)pixel.GetIdxMaxHiGainSample();
533 const Float_t timelo = (Float_t)pixel.GetIdxMaxLoGainSample(lofirst);
534
535 histhi.FillAbsTime(timehi);
536 histlo.FillAbsTime(timelo);
537
538 const Int_t aidx = (*fGeom)[pixid].GetAidx();
539 const Int_t sector = (*fGeom)[pixid].GetSector();
540
541 timehiarea[aidx] += timehi;
542 timeloarea[aidx] += timelo;
543
544 timehisector[sector] += timehi;
545 timelosector[sector] += timelo;
546 }
547
548 for (UInt_t j=0; j<nareas; j++)
549 {
550
551 const Int_t npix = fAverageAreaNum[j];
552
553 if (npix == 0)
554 continue;
555
556 MHCalibrationChargePix &hipix = (MHCalibrationChargePix&)GetAverageHiGainArea(j);
557 MHCalibrationChargePix &lopix = (MHCalibrationChargePix&)GetAverageLoGainArea(j);
558
559 hipix.FillHistAndArray(sumhiarea[j]/npix);
560 lopix.FillHistAndArray(sumloarea[j]/npix);
561
562 hipix.SetSaturated(sathiarea[j]/npix);
563 lopix.SetSaturated(satloarea[j]/npix);
564
565 hipix.FillAbsTime(timehiarea[j]/npix);
566 lopix.FillAbsTime(timeloarea[j]/npix);
567
568 }
569
570 for (UInt_t j=0; j<nsectors; j++)
571 {
572
573 const Int_t npix = fAverageSectorNum[j];
574
575 if (npix == 0)
576 continue;
577
578 MHCalibrationChargePix &hipix = (MHCalibrationChargePix&)GetAverageHiGainSector(j);
579 MHCalibrationChargePix &lopix = (MHCalibrationChargePix&)GetAverageLoGainSector(j);
580
581 hipix.FillHistAndArray(sumhisector[j]/npix);
582 lopix.FillHistAndArray(sumlosector[j]/npix);
583
584 hipix.SetSaturated(sathisector[j]/npix);
585 lopix.SetSaturated(satlosector[j]/npix);
586
587 hipix.FillAbsTime(timehisector[j]/npix);
588 lopix.FillAbsTime(timelosector[j]/npix);
589
590 }
591
592 return kTRUE;
593}
594
595// --------------------------------------------------------------------------
596//
597// For all TObjArray's (including the averaged ones), the following steps are performed:
598//
599// 1) Returns if the pixel is excluded.
600// 2) Tests saturation. In case yes, set the flag: MCalibrationPix::SetHiGainSaturation()
601// or the flag: MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainSaturated )
602// 3) Store the absolute arrival times in the MCalibrationChargePix's. If flag
603// MCalibrationPix::IsHiGainSaturation() is set, the Low-Gain arrival times are stored,
604// otherwise the Hi-Gain ones.
605// 4) Calls to MHCalibrationCam::FitHiGainArrays() and MCalibrationCam::FitLoGainArrays()
606// with the flags:
607// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted )
608// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainNotFitted )
609// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainOscillating )
610// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainOscillating )
611//
612Bool_t MHCalibrationChargeCam::FinalizeHists()
613{
614
615 *fLog << endl;
616
617 if (fHiGainOverFlow)
618 *fLog << warn << GetDescriptor()
619 << ": Histogram Overflow has occurred " << fHiGainOverFlow << " in the High-Gain! " << endl;
620 if (fLoGainOverFlow)
621 *fLog << warn << GetDescriptor()
622 << ": Histogram Overflow has occurred " << fLoGainOverFlow << " in the Low-Gain! " << endl;
623
624 for (Int_t i=0; i<fHiGainArray->GetSize(); i++)
625 {
626
627 MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)(*this)[i];
628 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
629 MBadPixelsPix &bad = (*fBadPixels)[i];
630
631 if (histhi.IsExcluded())
632 continue;
633
634 if (histhi.GetSaturated() > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries())
635 {
636 pix.SetHiGainSaturation();
637 histhi.CreateFourierSpectrum();
638 continue;
639 }
640
641 FinalizeAbsTimes(histhi, pix, bad, fFirstHiGain, fLastHiGain);
642 }
643
644 for (Int_t i=0; i<fLoGainArray->GetSize(); i++)
645 {
646
647 MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)(*this)(i);
648 MBadPixelsPix &bad = (*fBadPixels)[i];
649
650 if (histlo.IsExcluded())
651 continue;
652
653 if (histlo.GetSaturated() > fNumLoGainSaturationLimit*histlo.GetHGausHist()->GetEntries())
654 {
655 *fLog << warn << "Saturated Lo Gain histogram in pixel: " << i << endl;
656 bad.SetUncalibrated( MBadPixelsPix::kLoGainSaturation );
657 histlo.CreateFourierSpectrum();
658 continue;
659 }
660
661 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
662
663 if (pix.IsHiGainSaturation())
664 FinalizeAbsTimes(histlo, pix, bad, fFirstLoGain, fLastLoGain);
665 }
666
667 for (Int_t j=0; j<fAverageHiGainAreas->GetSize(); j++)
668 {
669
670 MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)GetAverageHiGainArea(j);
671 MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageArea(j);
672 MBadPixelsPix &bad = fCam->GetAverageBadArea(j);
673
674 if (histhi.GetSaturated() > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries())
675 {
676 pix.SetHiGainSaturation();
677 histhi.CreateFourierSpectrum();
678 continue;
679 }
680
681 FinalizeAbsTimes(histhi, pix, bad, fFirstHiGain, fLastHiGain);
682 }
683
684 for (Int_t j=0; j<fAverageLoGainAreas->GetSize(); j++)
685 {
686
687 MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)GetAverageLoGainArea(j);
688 MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageArea(j);
689 MBadPixelsPix &bad = fCam->GetAverageBadArea(j);
690
691 if (histlo.GetSaturated() > fNumLoGainSaturationLimit*histlo.GetHGausHist()->GetEntries())
692 {
693 *fLog << warn << "Saturated Lo Gain histogram in area idx: " << j << endl;
694 histlo.CreateFourierSpectrum();
695 continue;
696 }
697
698 if (pix.IsHiGainSaturation())
699 FinalizeAbsTimes(histlo, pix, bad, fFirstLoGain, fLastLoGain);
700 }
701
702 for (Int_t j=0; j<fAverageHiGainSectors->GetSize(); j++)
703 {
704
705 MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)GetAverageHiGainSector(j);
706 MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageSector(j);
707 MBadPixelsPix &bad = fCam->GetAverageBadSector(j);
708
709 if (histhi.GetSaturated() > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries())
710 {
711 pix.SetHiGainSaturation();
712 histhi.CreateFourierSpectrum();
713 continue;
714 }
715
716 FinalizeAbsTimes(histhi, pix, bad, fFirstHiGain, fLastHiGain);
717 }
718
719 for (Int_t j=0; j<fAverageLoGainSectors->GetSize(); j++)
720 {
721
722 MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)GetAverageLoGainSector(j);
723 MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageSector(j);
724 MBadPixelsPix &bad = fCam->GetAverageBadSector(j);
725
726 if (histlo.GetSaturated() > fNumLoGainSaturationLimit*histlo.GetHGausHist()->GetEntries())
727 {
728 *fLog << warn << "Saturated Lo Gain histogram in sector: " << j << endl;
729 bad.SetUncalibrated( MBadPixelsPix::kLoGainSaturation );
730 histlo.CreateFourierSpectrum();
731 continue;
732 }
733
734 if (pix.IsHiGainSaturation())
735 FinalizeAbsTimes(histlo, pix, bad, fFirstLoGain, fLastLoGain);
736 }
737
738 //
739 // Perform the fitting for the High Gain (done in MHCalibrationCam)
740 //
741 FitHiGainArrays((MCalibrationCam&)(*fCam),(*fBadPixels),
742 MBadPixelsPix::kHiGainNotFitted,
743 MBadPixelsPix::kHiGainOscillating);
744 //
745 // Perform the fitting for the Low Gain (done in MHCalibrationCam)
746 //
747 FitLoGainArrays((MCalibrationCam&)(*fCam),(*fBadPixels),
748 MBadPixelsPix::kLoGainNotFitted,
749 MBadPixelsPix::kLoGainOscillating);
750
751 return kTRUE;
752}
753
754// --------------------------------------------------------------------------------
755//
756// Fill the absolute time results into MCalibrationChargePix
757//
758// Check absolute time validity:
759// - Mean arrival time is at least fTimeLowerLimit slices from the lower edge
760// - Mean arrival time is at least fUpperLimit slices from the upper edge
761//
762void MHCalibrationChargeCam::FinalizeAbsTimes(MHCalibrationChargePix &hist, MCalibrationChargePix &pix, MBadPixelsPix &bad,
763 Byte_t first, Byte_t last)
764{
765
766 const Float_t mean = hist.GetAbsTimeMean();
767 const Float_t rms = hist.GetAbsTimeRms();
768
769 pix.SetAbsTimeMean ( mean );
770 pix.SetAbsTimeRms ( rms );
771
772 const Float_t lowerlimit = (Float_t)first + fTimeLowerLimit;
773 const Float_t upperlimit = (Float_t)last + fTimeUpperLimit;
774
775 if ( mean < lowerlimit)
776 {
777 *fLog << warn << GetDescriptor()
778 << Form("%s%3.1f%s%2.1f%s%4i",": Mean ArrivalTime: ",mean," smaller than ",fTimeLowerLimit,
779 " FADC slices from lower edge in pixel ",hist.GetPixId()) << endl;
780 bad.SetUncalibrated( MBadPixelsPix::kMeanTimeInFirstBin );
781 }
782
783 if ( mean > upperlimit )
784 {
785 *fLog << warn << GetDescriptor()
786 << Form("%s%3.1f%s%2.1f%s%4i",": Mean ArrivalTime: ",mean," greater than ",fTimeUpperLimit,
787 " FADC slices from upper edge in pixel ",hist.GetPixId()) << endl;
788 bad.SetUncalibrated( MBadPixelsPix::kMeanTimeInLast2Bins );
789 }
790}
791
792// --------------------------------------------------------------------------
793//
794// Sets all pixels to MBadPixelsPix::kUnsuitableRun, if following flags are set:
795// - MBadPixelsPix::kLoGainSaturation
796//
797// Sets all pixels to MBadPixelsPix::kUnreliableRun, if following flags are set:
798// - if MBadPixelsPix::kHiGainNotFitted and !MCalibrationPix::IsHiGainSaturation()
799// - if MBadPixelsPix::kHiGainOscillating and !MCalibrationPix::IsHiGainSaturation()
800// - if MBadPixelsPix::kLoGainNotFitted and MCalibrationPix::IsLoGainSaturation()
801// - if MBadPixelsPix::kLoGainOscillating and MCalibrationPix::IsLoGainSaturation()
802//
803void MHCalibrationChargeCam::FinalizeBadPixels()
804{
805
806 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
807 {
808
809 MBadPixelsPix &bad = (*fBadPixels)[i];
810 MCalibrationPix &pix = (*fCam)[i];
811
812 if (bad.IsUncalibrated( MBadPixelsPix::kHiGainNotFitted ))
813 if (!pix.IsHiGainSaturation())
814 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun );
815
816 if (bad.IsUncalibrated( MBadPixelsPix::kHiGainOscillating ))
817 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun );
818
819 if (bad.IsUncalibrated( MBadPixelsPix::kLoGainNotFitted ))
820 if (pix.IsHiGainSaturation())
821 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun );
822
823 if (bad.IsUncalibrated( MBadPixelsPix::kLoGainOscillating ))
824 if (pix.IsHiGainSaturation())
825 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun );
826
827 if (bad.IsUncalibrated( MBadPixelsPix::kLoGainSaturation ))
828 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );
829 }
830}
831
832// --------------------------------------------------------------------------
833//
834// Dummy, needed by MCamEvent
835//
836Bool_t MHCalibrationChargeCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
837{
838 return kTRUE;
839}
840
841// --------------------------------------------------------------------------
842//
843// Calls MHGausEvents::DrawClone() for pixel idx
844//
845void MHCalibrationChargeCam::DrawPixelContent(Int_t idx) const
846{
847 (*this)[idx].DrawClone();
848}
849
850
851// -----------------------------------------------------------------------------
852//
853// Default draw:
854//
855// Displays the averaged areas, both High Gain and Low Gain
856//
857// Calls the Draw of the fAverageHiGainAreas and fAverageLoGainAreas objects with options
858//
859void MHCalibrationChargeCam::Draw(const Option_t *opt)
860{
861
862 const Int_t nareas = fAverageHiGainAreas->GetEntries();
863 if (nareas == 0)
864 return;
865
866 TString option(opt);
867 option.ToLower();
868
869 if (!option.Contains("datacheck"))
870 {
871 MHCalibrationCam::Draw(opt);
872 return;
873 }
874
875 //
876 // From here on , the datacheck - Draw
877 //
878 TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas(this);
879 pad->SetBorderMode(0);
880 pad->Divide(1,nareas);
881
882 //
883 // Loop over inner and outer pixels
884 //
885 for (Int_t i=0; i<nareas;i++)
886 {
887
888 pad->cd(i+1);
889
890 MHCalibrationChargePix &hipix = (MHCalibrationChargePix&)GetAverageHiGainArea(i);
891 MHCalibrationChargePix &lopix = (MHCalibrationChargePix&)GetAverageLoGainArea(i);
892 //
893 // Ask for Hi-Gain saturation
894 //
895 if (hipix.GetSaturated() > fNumHiGainSaturationLimit*hipix.GetHGausHist()->GetEntries())
896 DrawDataCheckPixel(lopix,i ? gkLoGainOuterRefLines : gkLoGainInnerRefLines);
897 else
898 DrawDataCheckPixel(hipix,i ? gkHiGainOuterRefLines : gkHiGainInnerRefLines);
899 }
900}
901
902
903// --------------------------------------------------------------------------
904//
905// Our own clone function is necessary since root 3.01/06 or Mars 0.4
906// I don't know the reason.
907//
908// Creates new MHCalibrationCam
909//
910TObject *MHCalibrationChargeCam::Clone(const char *name) const
911{
912
913 const Int_t navhi = fAverageHiGainAreas->GetEntries();
914 const Int_t navlo = fAverageLoGainAreas->GetEntries();
915 const Int_t nsehi = fAverageHiGainSectors->GetEntries();
916 const Int_t nselo = fAverageLoGainSectors->GetEntries();
917
918 //
919 // FIXME, this might be done faster and more elegant, by direct copy.
920 //
921 MHCalibrationChargeCam *cam = new MHCalibrationChargeCam();
922
923 cam->fAverageHiGainAreas->Expand(navhi);
924 cam->fAverageLoGainAreas->Expand(navlo);
925 cam->fAverageHiGainSectors->Expand(nsehi);
926 cam->fAverageLoGainSectors->Expand(nselo);
927
928 cam->fAverageHiGainAreas->Expand(navhi);
929 cam->fAverageLoGainAreas->Expand(navlo);
930 cam->fAverageHiGainSectors->Expand(nsehi);
931 cam->fAverageLoGainSectors->Expand(nselo);
932
933 for (int i=0; i<navhi; i++)
934 (*cam->fAverageHiGainAreas) [i] = (*fAverageHiGainAreas) [i]->Clone();
935 for (int i=0; i<navlo; i++)
936 (*cam->fAverageLoGainAreas) [i] = (*fAverageLoGainAreas) [i]->Clone();
937 for (int i=0; i<nsehi; i++)
938 (*cam->fAverageHiGainSectors)[i] = (*fAverageHiGainSectors)[i]->Clone();
939 for (int i=0; i<nselo; i++)
940 (*cam->fAverageLoGainSectors)[i] = (*fAverageLoGainSectors)[i]->Clone();
941
942 cam->fAverageAreaNum = fAverageAreaNum;
943 cam->fAverageAreaSat = fAverageAreaSat;
944 cam->fAverageAreaSigma = fAverageAreaSigma;
945 cam->fAverageAreaSigmaVar = fAverageAreaSigmaVar;
946 cam->fAverageAreaRelSigma = fAverageAreaRelSigma;
947 cam->fAverageAreaRelSigmaVar = fAverageAreaRelSigmaVar;
948 cam->fAverageSectorNum = fAverageSectorNum;
949 cam->fRunNumbers = fRunNumbers;
950
951 cam->fColor = fColor;
952 cam->fPulserFrequency = fPulserFrequency;
953 cam->fAverageNbins = fAverageNbins;
954
955 return cam;
956
957}
958
959void MHCalibrationChargeCam::DrawDataCheckPixel(MHCalibrationChargePix &pix, const Float_t refline[])
960{
961
962 TVirtualPad *newpad = gPad;
963 newpad->Divide(1,2);
964 newpad->cd(1);
965
966 gPad->SetTicks();
967 if (!pix.IsEmpty() && !pix.IsOnlyOverflow())
968 gPad->SetLogy();
969
970 gStyle->SetOptStat(0);
971
972 TH1F *hist = pix.GetHGausHist();
973
974 TH2D *null = new TH2D("Null",hist->GetTitle(),100,pix.GetFirst() > -1. ? 0. : 100.,pix.GetLast()/2.,
975 100,0.,hist->GetEntries()/10.);
976
977 null->SetDirectory(NULL);
978 null->SetBit(kCanDelete);
979 //
980 // set the labels bigger
981 //
982 TAxis *xaxe = null->GetXaxis();
983 TAxis *yaxe = null->GetYaxis();
984 xaxe->CenterTitle();
985 yaxe->CenterTitle();
986 xaxe->SetTitleSize(0.07);
987 yaxe->SetTitleSize(0.07);
988 xaxe->SetTitleOffset(0.7);
989 yaxe->SetTitleOffset(0.55);
990 xaxe->SetLabelSize(0.06);
991 yaxe->SetLabelSize(0.06);
992
993 xaxe->SetTitle(hist->GetXaxis()->GetTitle());
994 yaxe->SetTitle(hist->GetYaxis()->GetTitle());
995 null->Draw();
996 hist->Draw("same");
997
998 gStyle->SetOptFit();
999
1000 if (pix.GetFGausFit())
1001 {
1002 switch ( fColor )
1003 {
1004 case MCalibrationCam::kGREEN:
1005 pix.GetFGausFit()->SetLineColor(kGreen);
1006 break;
1007 case MCalibrationCam::kBLUE:
1008 pix.GetFGausFit()->SetLineColor(kBlue);
1009 break;
1010 case MCalibrationCam::kUV:
1011 pix.GetFGausFit()->SetLineColor(106);
1012 break;
1013 case MCalibrationCam::kCT1:
1014 pix.GetFGausFit()->SetLineColor(006);
1015 break;
1016 default:
1017 pix.GetFGausFit()->SetLineColor(kRed);
1018 }
1019 pix.GetFGausFit()->Draw("same");
1020 }
1021
1022 DisplayRefLines(null,refline);
1023
1024 newpad->cd(2);
1025 gPad->SetTicks();
1026
1027 pix.DrawEvents();
1028 return;
1029
1030}
1031
1032
1033void MHCalibrationChargeCam::DisplayRefLines(const TH2D *hist, const Float_t refline[]) const
1034{
1035
1036 TLine *green1 = new TLine(refline[0],0.,refline[0],hist->GetYaxis()->GetXmax());
1037 green1->SetBit(kCanDelete);
1038 green1->SetLineColor(kGreen);
1039 green1->SetLineStyle(2);
1040 green1->SetLineWidth(3);
1041 green1->Draw();
1042
1043 TLine *green5 = new TLine(refline[6],0.,refline[6],hist->GetYaxis()->GetXmax());
1044 green5->SetBit(kCanDelete);
1045 green5->SetLineColor(8);
1046 green5->SetLineStyle(2);
1047 green5->SetLineWidth(3);
1048 green5->Draw();
1049
1050 TLine *blue1 = new TLine(refline[1],0.,refline[1],hist->GetYaxis()->GetXmax());
1051 blue1->SetBit(kCanDelete);
1052 blue1->SetLineColor(227);
1053 blue1->SetLineStyle(2);
1054 blue1->SetLineWidth(3);
1055 blue1->Draw();
1056
1057 TLine *blue5 = new TLine(refline[2],0.,refline[2],hist->GetYaxis()->GetXmax());
1058 blue5->SetBit(kCanDelete);
1059 blue5->SetLineColor(68);
1060 blue5->SetLineStyle(2);
1061 blue5->SetLineWidth(3);
1062 blue5->Draw();
1063
1064 TLine *blue10 = new TLine(refline[3],0.,refline[3],hist->GetYaxis()->GetXmax());
1065 blue10->SetBit(kCanDelete);
1066 blue10->SetLineColor(4);
1067 blue10->SetLineStyle(2);
1068 blue10->SetLineWidth(3);
1069 blue10->Draw();
1070
1071 TLine *uv10 = new TLine(refline[4],0.,refline[4],hist->GetYaxis()->GetXmax());
1072 uv10->SetBit(kCanDelete);
1073 uv10->SetLineColor(106);
1074 uv10->SetLineStyle(2);
1075 uv10->SetLineWidth(3);
1076 uv10->Draw();
1077
1078 TLine *ct1 = new TLine(refline[5],0.,refline[5],hist->GetYaxis()->GetXmax());
1079 ct1->SetBit(kCanDelete);
1080 ct1->SetLineColor(6);
1081 ct1->SetLineStyle(2);
1082 ct1->SetLineWidth(3);
1083 ct1->Draw();
1084
1085 TLegend *leg = new TLegend(0.7,0.35,0.9,0.99);
1086 leg->SetBit(kCanDelete);
1087 leg->AddEntry(green1,"1 Led GREEN","l");
1088 leg->AddEntry(green5,"5 Leds GREEN","l");
1089 leg->AddEntry(blue1,"1 Led BLUE","l");
1090 leg->AddEntry(blue5,"5 Leds BLUE","l");
1091 leg->AddEntry(blue10,"10 Leds BLUE","l");
1092 leg->AddEntry(uv10,"10 Leds UV","l");
1093 leg->AddEntry(ct1,"CT1-Pulser","l");
1094
1095 leg->Draw();
1096}
Note: See TracBrowser for help on using the repository browser.