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

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