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

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