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

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