source: branches/MarsMoreSimulationTruth/mhcalib/MHCalibrationHiLoCam.cc

Last change on this file was 10166, checked in by tbretz, 14 years ago
Removed the old obsolete cvs header line.
File size: 32.7 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// MHCalibrationHiLoCam
27//
28// Fills the extracted high-gain low-gain charge ratios of MArrivalTimeCam into
29// the MHCalibrationPix-classes MHCalibrationPix for every:
30//
31// - Pixel, stored in the TObjArray's MHCalibrationCam::fHiGainArray
32// or MHCalibrationCam::fHiGainArray, respectively, depending if
33// MArrivalTimePix::IsLoGainUsed() is set.
34//
35// - Average pixel per AREA index (e.g. inner and outer for the MAGIC camera),
36// stored in the TObjArray's MHCalibrationCam::fAverageHiGainAreas and
37// MHCalibrationCam::fAverageHiGainAreas
38//
39// - Average pixel per camera SECTOR (e.g. sectors 1-6 for the MAGIC camera),
40// stored in the TObjArray's MHCalibrationCam::fAverageHiGainSectors
41// and MHCalibrationCam::fAverageHiGainSectors
42//
43// The histograms are fitted to a Gaussian, mean and sigma with its errors
44// and the fit probability are extracted. If none of these values are NaN's and
45// if the probability is bigger than MHGausEvents::fProbLimit (default: 0.5%),
46// the fit is declared valid.
47// Otherwise, the fit is repeated within ranges of the previous mean
48// +- MHCalibrationPix::fPickupLimit (default: 5) sigma (see MHCalibrationPix::RepeatFit())
49// In case this does not make the fit valid, the histogram means and RMS's are
50// taken directly (see MHCalibrationPix::BypassFit()) and the following flags are set:
51// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiLoNotFitted ) and
52// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun )
53//
54// Outliers of more than MHCalibrationPix::fPickupLimit (default: 5) sigmas
55// from the mean are counted as Pickup events (stored in MHCalibrationPix::fPickup)
56//
57// The class also fills arrays with the signal vs. event number, creates a fourier
58// spectrum (see MHGausEvents::CreateFourierSpectrum()) and investigates if the
59// projected fourier components follow an exponential distribution.
60// In case that the probability of the exponential fit is less than
61// MHGausEvents::fProbLimit (default: 0.5%), the following flags are set:
62// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiLoOscillating ) and
63// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun )
64//
65// This same procedure is performed for the average pixels.
66//
67// The following results are written into MCalibrationHiLoCam:
68//
69// - MCalibrationPix::SetMean()
70// - MCalibrationPix::SetMeanErr()
71// - MCalibrationPix::SetSigma()
72// - MCalibrationPix::SetSigmaErr()
73// - MCalibrationPix::SetProb()
74// - MCalibrationPix::SetNumPickup()
75//
76// For all averaged areas, the fitted sigma is multiplied with the square root of
77// the number involved pixels in order to be able to compare it to the average of
78// sigmas in the camera.
79//
80/////////////////////////////////////////////////////////////////////////////
81#include "MHCalibrationHiLoCam.h"
82
83#include <TOrdCollection.h>
84#include <TPad.h>
85#include <TVirtualPad.h>
86#include <TCanvas.h>
87#include <TStyle.h>
88#include <TF1.h>
89#include <TLine.h>
90#include <TLatex.h>
91#include <TLegend.h>
92#include <TGraph.h>
93#include <TProfile.h>
94
95#include "MHCalibrationHiLoPix.h"
96
97#include "MLog.h"
98#include "MLogManip.h"
99
100#include "MParList.h"
101
102#include "MCalibrationHiLoCam.h"
103#include "MCalibrationHiLoPix.h"
104#include "MCalibrationCam.h"
105#include "MCalibrationPix.h"
106
107#include "MExtractedSignalCam.h"
108#include "MExtractedSignalPix.h"
109#include "MArrivalTimeCam.h"
110#include "MArrivalTimePix.h"
111
112#include "MGeomCam.h"
113#include "MGeomPix.h"
114
115#include "MBadPixelsCam.h"
116#include "MBadPixelsPix.h"
117
118ClassImp(MHCalibrationHiLoCam);
119
120using namespace std;
121
122const Int_t MHCalibrationHiLoCam::fgNbins = 175;
123const Axis_t MHCalibrationHiLoCam::fgFirst = -5.1;
124const Axis_t MHCalibrationHiLoCam::fgLast = 29.9;
125const Float_t MHCalibrationHiLoCam::fgProbLimit = 0.;
126const Int_t MHCalibrationHiLoCam::fgHivsLoNbins = 90;
127const Axis_t MHCalibrationHiLoCam::fgHivsLoFirst = 95.;
128const Axis_t MHCalibrationHiLoCam::fgHivsLoLast = 995.;
129const Axis_t MHCalibrationHiLoCam::fgLowerFitLimitProfile = 480.;
130const Axis_t MHCalibrationHiLoCam::fgUpperFitLimitProfile = 680.;
131const TString MHCalibrationHiLoCam::gsHistName = "HiLo";
132const TString MHCalibrationHiLoCam::gsHistTitle = "HiGain vs. LoGain";
133const TString MHCalibrationHiLoCam::gsHistXTitle = "Amplification Ratio [1]";
134const TString MHCalibrationHiLoCam::gsHistYTitle = "Nr. events";
135const TString MHCalibrationHiLoCam::gsHivsLoHistName = "HivsLo";
136const TString MHCalibrationHiLoCam::gsHivsLoHistTitle = "High-gain vs. Low-gain Charge";
137const TString MHCalibrationHiLoCam::gsHivsLoHistXTitle = "Q High-Gain [FADC counts]";
138const TString MHCalibrationHiLoCam::gsHivsLoHistYTitle = "Q Low-Gain [FADC counts]";
139
140// --------------------------------------------------------------------------
141//
142// Default Constructor.
143//
144// Sets:
145// - fNbins to fgNbins
146// - fFirst to fgFirst
147// - fLast to fgLast
148//
149// - fHistName to gsHistName
150// - fHistTitle to gsHistTitle
151// - fHistXTitle to gsHistXTitle
152// - fHistYTitle to gsHistYTitle
153//
154// - fLowerLimt to fgLowerLim
155// - fUpperLimt to fgUpperLim
156//
157MHCalibrationHiLoCam::MHCalibrationHiLoCam(const char *name, const char *title)
158 : fArrTimes(NULL), fHivsLoResults("Results","Fit Results high-gain vs. low-gain",
159 200,-10.,10.,200,0.,20.),
160 fUsedLoGainSlices(0)
161{
162
163 fName = name ? name : "MHCalibrationHiLoCam";
164 fTitle = title ? title : "Histogram class for the high-gain vs. low-gain amplification ratio calibration";
165
166 SetNbins(fgNbins);
167 SetFirst(fgFirst);
168 SetLast (fgLast );
169
170 SetProbLimit(fgProbLimit);
171
172 SetHistName (gsHistName .Data());
173 SetHistTitle (gsHistTitle .Data());
174 SetHistXTitle(gsHistXTitle.Data());
175 SetHistYTitle(gsHistYTitle.Data());
176
177 SetHivsLoNbins(fgHivsLoNbins);
178 SetHivsLoFirst(fgHivsLoFirst);
179 SetHivsLoLast (fgHivsLoLast );
180
181 SetLowerFitLimitProfile();
182 SetUpperFitLimitProfile();
183
184 SetHivsLoHistName (gsHivsLoHistName .Data());
185 SetHivsLoHistTitle (gsHivsLoHistTitle .Data());
186 SetHivsLoHistXTitle(gsHivsLoHistXTitle.Data());
187 SetHivsLoHistYTitle(gsHivsLoHistYTitle.Data());
188
189 SetOscillations(kFALSE);
190
191 fHivsLoResults.GetXaxis()->SetTitle("Offset per FADC slices [FADC cnts]");
192 fHivsLoResults.GetYaxis()->SetTitle("Gains ratio [1]");
193 fHivsLoResults.SetDirectory(0);
194
195}
196
197// --------------------------------------------------------------------------
198//
199// Creates new MHCalibrationHiLoCam only with the averaged areas:
200// the rest has to be retrieved directly, e.g. via:
201// MHCalibrationHiLoCam *cam = MParList::FindObject("MHCalibrationHiLoCam");
202// - cam->GetAverageSector(5).DrawClone();
203// - (*cam)[100].DrawClone()
204//
205TObject *MHCalibrationHiLoCam::Clone(const char *) const
206{
207
208 MHCalibrationHiLoCam *cam = new MHCalibrationHiLoCam();
209
210 //
211 // Copy the data members
212 //
213 cam->fColor = fColor;
214 cam->fRunNumbers = fRunNumbers;
215 cam->fPulserFrequency = fPulserFrequency;
216 cam->fFlags = fFlags;
217 cam->fNbins = fNbins;
218 cam->fFirst = fFirst;
219 cam->fLast = fLast;
220
221 //
222 // Copy the MArrays
223 //
224 cam->fAverageAreaRelSigma = fAverageAreaRelSigma;
225 cam->fAverageAreaRelSigmaVar = fAverageAreaRelSigmaVar;
226 cam->fAverageAreaSat = fAverageAreaSat;
227 cam->fAverageAreaSigma = fAverageAreaSigma;
228 cam->fAverageAreaSigmaVar = fAverageAreaSigmaVar;
229 cam->fAverageAreaNum = fAverageAreaNum;
230 cam->fAverageSectorNum = fAverageSectorNum;
231
232 if (!IsAverageing())
233 return cam;
234
235 const Int_t navhi = fAverageHiGainAreas->GetSize();
236 const Int_t navlo = fAverageLoGainAreas->GetSize();
237
238 for (int i=0; i<navhi; i++)
239 cam->fAverageHiGainAreas->AddAt(GetAverageHiGainArea(i).Clone(),i);
240
241 for (int i=0; i<navlo; i++)
242 cam->fAverageLoGainAreas->AddAt(GetAverageLoGainArea(i).Clone(),i);
243
244 return cam;
245}
246
247// --------------------------------------------------------------------------
248//
249// Gets or creates the pointers to:
250// - MCalibrationHiLoCam
251//
252// Searches pointer to:
253// - MExtractedSignalCam
254// - MArrivalTimeCam
255//
256// Calls:
257// - MHCalibrationCam::InitHiGainArrays()
258// - MHCalibrationCam::InitLoGainArrays()
259//
260// Sets:
261// - fSumarea to nareas
262// - fSumsector to nareas
263// - fNumarea to nareas
264// - fNumsector to nareas
265//
266Bool_t MHCalibrationHiLoCam::ReInitHists(MParList *pList)
267{
268
269 fCam = (MCalibrationCam*)pList->FindObject(AddSerialNumber("MCalibrationHiLoCam"));
270 if (!fCam)
271 {
272 fCam = (MCalibrationCam*)pList->FindCreateObj(AddSerialNumber("MCalibrationHiLoCam"));
273 if (!fCam)
274 return kFALSE;
275 fCam->Init(*fGeom);
276 }
277
278 MExtractedSignalCam *signal = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
279 if (!signal)
280 {
281 *fLog << err << "MExtractedSignalCam not found... abort." << endl;
282 return kFALSE;
283 }
284
285 fUsedLoGainSlices = signal->GetNumUsedLoGainFADCSlices();
286
287 fArrTimes = (MArrivalTimeCam*)pList->FindObject("MArrivalTimeCam");
288 if (!fArrTimes)
289 {
290 *fLog << warn << "MArrivalTimeCam not found... cannot calibrated arrival times between "
291 <<"high and low-gain" << endl;
292 SetLoGain(kFALSE);
293 }
294
295 const Int_t npixels = fGeom->GetNumPixels();
296 const Int_t nsectors = fGeom->GetNumSectors();
297 const Int_t nareas = fGeom->GetNumAreas();
298
299 InitHiGainArrays(npixels,nareas,nsectors);
300 InitLoGainArrays(npixels,nareas,nsectors);
301
302 fSumareahi .Set(nareas);
303 fSumsectorhi.Set(nsectors);
304 fNumareahi .Set(nareas);
305 fNumsectorhi.Set(nsectors);
306 if (IsLoGain())
307 {
308 fSumarealo .Set(nareas);
309 fSumsectorlo.Set(nsectors);
310 fNumarealo .Set(nareas);
311 fNumsectorlo.Set(nsectors);
312 }
313 return kTRUE;
314}
315
316// --------------------------------------------------------------------------
317//
318// Retrieve:
319// - fRunHeader->GetNumSamplesHiGain();
320//
321// Initializes the High Gain Arrays:
322//
323// - For every entry in the expanded arrays:
324// * Initialize an MHCalibrationHiLoPix
325// * Set Binning from fNbins, fFirst and fLast
326// * Set Binning of Abs Times histogram from fAbsNbins, fAbsFirst and fAbsLast
327// * Set Histgram names and titles from fHistName and fHistTitle
328// * Set Abs Times Histgram names and titles from fAbsHistName and fAbsHistTitle
329// * Set X-axis and Y-axis titles from fHistXTitle and fHistYTitle
330// * Set X-axis and Y-axis titles of Abs Times Histogram from fAbsHistXTitle and fAbsHistYTitle
331// * Call InitHists
332//
333//
334void MHCalibrationHiLoCam::InitHiGainArrays(const Int_t npixels, const Int_t nareas, const Int_t nsectors)
335{
336
337 TProfile *h;
338
339 if (fHiGainArray->GetSize()==0)
340 {
341 for (Int_t i=0; i<npixels; i++)
342 {
343 fHiGainArray->AddAt(new MHCalibrationHiLoPix(Form("%sHiGainPix%04d",fHistName.Data(),i),
344 Form("%s High Gain Pixel%04d",fHistTitle.Data(),i)),i);
345
346 MHCalibrationHiLoPix &pix = (MHCalibrationHiLoPix&)(*this)[i];
347
348 pix.SetNbins(fNbins);
349 pix.SetFirst(fFirst);
350 pix.SetLast (fLast);
351
352 pix.SetProbLimit(fProbLimit);
353
354 pix.SetHivsLoNbins(fHivsLoNbins);
355 pix.SetHivsLoFirst(fHivsLoFirst);
356 pix.SetHivsLoLast (fHivsLoLast);
357
358 InitHists(pix,(*fBadPixels)[i],i);
359
360 if (fCam)
361 (*fCam)[i].SetPixId(i);
362
363 h = pix.GetHivsLo();
364
365 h->SetName (Form("H%sHiGainPix%04d",fHivsLoHistName.Data(),i));
366 h->SetTitle(Form("%s High Gain Pixel %04d",fHivsLoHistTitle.Data(),i));
367 h->SetXTitle(fHivsLoHistXTitle.Data());
368 h->SetYTitle(fHivsLoHistYTitle.Data());
369 h->SetDirectory(0);
370 }
371 }
372
373
374 if (fAverageHiGainAreas->GetSize()==0)
375 {
376 for (Int_t j=0; j<nareas; j++)
377 {
378 fAverageHiGainAreas->AddAt(new MHCalibrationHiLoPix(Form("%sHiGainArea%d",fHistName.Data(),j),
379 Form("%s High Gain Area Idx %d",fHistTitle.Data(),j)),j);
380
381 MHCalibrationHiLoPix &pix = (MHCalibrationHiLoPix&)GetAverageHiGainArea(j);
382
383 pix.SetNbins(fNbins*(Int_t)TMath::Sqrt((Float_t)npixels/nareas));
384 pix.SetFirst(fFirst);
385 pix.SetLast (fLast);
386
387 pix.SetHivsLoNbins(fHivsLoNbins);
388 pix.SetHivsLoFirst(fHivsLoFirst);
389 pix.SetHivsLoLast (fHivsLoLast);
390
391 InitHists(pix, fCam->GetAverageBadArea(j), j);
392
393 if (fCam)
394 fCam->GetAverageArea(j).SetPixId(j);
395
396 h = pix.GetHivsLo();
397
398 h->SetName (Form("H%sHiGainArea%d",fHivsLoHistName.Data(),j));
399 h->SetTitle(Form("%s averaged on event-by-event basis High Gain Area Idx %d",
400 fHivsLoHistTitle.Data(), j));
401 h->SetXTitle(fHivsLoHistXTitle.Data());
402 h->SetYTitle(fHivsLoHistYTitle.Data());
403 h->SetDirectory(0);
404 }
405 }
406
407 if (fAverageHiGainSectors->GetSize()==0)
408 {
409 for (Int_t j=0; j<nsectors; j++)
410 {
411 fAverageHiGainSectors->AddAt(new MHCalibrationHiLoPix(Form("%sHiGainSector%02d",fHistName.Data(),j),
412 Form("%s High Gain Sector %02d",fHistTitle.Data(),j)),j);
413
414 MHCalibrationHiLoPix &pix = (MHCalibrationHiLoPix&)GetAverageHiGainSector(j);
415
416 pix.SetNbins(fNbins*(Int_t)TMath::Sqrt((Float_t)npixels/nareas));
417 pix.SetFirst(fFirst);
418 pix.SetLast (fLast);
419
420 pix.SetHivsLoNbins(fHivsLoNbins);
421 pix.SetHivsLoFirst(fHivsLoFirst);
422 pix.SetHivsLoLast (fHivsLoLast);
423
424 InitHists(pix, fCam->GetAverageBadSector(j), j);
425
426 if (fCam)
427 fCam->GetAverageSector(j).SetPixId(j);
428
429 h = pix.GetHivsLo();
430
431 h->SetName (Form("H%sHiGainSector%02d",fHivsLoHistName.Data(),j));
432 h->SetTitle(Form("%s averaged on event-by-event basis High Gain Area Sector %02d",
433 fHivsLoHistTitle.Data(),j));
434 h->SetXTitle(fHivsLoHistXTitle.Data());
435 h->SetYTitle(fHivsLoHistYTitle.Data());
436 h->SetDirectory(0);
437 }
438 }
439}
440
441//--------------------------------------------------------------------------------------
442//
443// Return, if IsLoGain() is kFALSE
444//
445// Retrieve:
446// - fRunHeader->GetNumSamplesHiGain();
447//
448// Initializes the Low Gain Arrays:
449//
450// - For every entry in the expanded arrays:
451// * Initialize an MHCalibrationHiLoPix
452// * Set Binning from fNbins, fFirst and fLast
453// * Set Binning of HivsLo Times histogram from fHivsLoNbins, fHivsLoFirst and fHivsLoLast
454// * Set Histgram names and titles from fHistName and fHistTitle
455// * Set HivsLo Times Histgram names and titles from fHivsLoHistName and fHivsLoHistTitle
456// * Set X-axis and Y-axis titles from fHistXTitle and fHistYTitle
457// * Set X-axis and Y-axis titles of HivsLo Times Histogram from fHivsLoHistXTitle and fHivsLoHistYTitle
458// * Call InitHists
459//
460void MHCalibrationHiLoCam::InitLoGainArrays(const Int_t npixels, const Int_t nareas, const Int_t nsectors)
461{
462 if (!IsLoGain())
463 return;
464
465 TProfile *h;
466
467 if (fLoGainArray->GetSize()==0 )
468 {
469 for (Int_t i=0; i<npixels; i++)
470 {
471 fLoGainArray->AddAt(new MHCalibrationHiLoPix(Form("%sLoGainPix%04d",fHistName.Data(),i),
472 Form("%s Low Gain Pixel %04d",fHistTitle.Data(),i)),i);
473
474 MHCalibrationHiLoPix &pix = (MHCalibrationHiLoPix&)(*this)(i);
475
476 pix.SetNbins(fNbins);
477 pix.SetFirst(fFirst);
478 pix.SetLast (fLast);
479
480 pix.SetProbLimit(fProbLimit);
481
482 pix.SetHivsLoNbins(fHivsLoNbins);
483 pix.SetHivsLoFirst(fHivsLoFirst);
484 pix.SetHivsLoLast (fHivsLoLast );
485
486 InitHists(pix,(*fBadPixels)[i],i);
487
488 h = pix.GetHivsLo();
489
490 h->SetName (Form("H%sLoGainPix%04d",fHivsLoHistName.Data(),i));
491 h->SetTitle(Form("%s Low Gain Pixel %04d",fHivsLoHistTitle.Data(),i));
492 h->SetXTitle(fHivsLoHistXTitle.Data());
493 h->SetYTitle(fHivsLoHistYTitle.Data());
494 h->SetDirectory(0);
495 }
496 }
497
498 if (fAverageLoGainAreas->GetSize()==0)
499 {
500 for (Int_t j=0; j<nareas; j++)
501 {
502 fAverageLoGainAreas->AddAt(new MHCalibrationHiLoPix(Form("%sLoGainArea%d",fHistName.Data(),j),
503 Form("%s Low Gain Area Idx %d",fHistTitle.Data(),j)),j);
504
505 MHCalibrationHiLoPix &pix = (MHCalibrationHiLoPix&)GetAverageLoGainArea(j);
506
507 pix.SetNbins(fNbins*(Int_t)TMath::Sqrt((Float_t)npixels/nareas));
508 pix.SetFirst(fFirst);
509 pix.SetLast (fLast);
510
511 pix.SetHivsLoNbins(fHivsLoNbins);
512 pix.SetHivsLoFirst(fHivsLoFirst);
513 pix.SetHivsLoLast (fHivsLoLast );
514
515 InitHists(pix, fCam->GetAverageBadArea(j), j);
516
517 h = pix.GetHivsLo();
518
519 h->SetName (Form("H%sLoGainArea%02d",fHivsLoHistName.Data(),j));
520 h->SetTitle(Form("%s%s%02d",fHivsLoHistTitle.Data(),
521 " averaged on event-by-event basis Low Gain Area Idx ",j));
522 h->SetXTitle(fHivsLoHistXTitle.Data());
523 h->SetYTitle(fHivsLoHistYTitle.Data());
524 h->SetDirectory(0);
525 }
526 }
527
528
529 if (fAverageLoGainSectors->GetSize()==0 && IsLoGain())
530 {
531 for (Int_t j=0; j<nsectors; j++)
532 {
533 fAverageLoGainSectors->AddAt(new MHCalibrationHiLoPix(Form("%sLoGainSector%02d",fHistName.Data(),j),
534 Form("%s Low Gain Sector %02d",fHistTitle.Data(),j)),j);
535
536 MHCalibrationHiLoPix &pix = (MHCalibrationHiLoPix&)GetAverageLoGainSector(j);
537
538 pix.SetNbins(fNbins*(Int_t)TMath::Sqrt((Float_t)npixels/nareas));
539 pix.SetFirst(fFirst);
540 pix.SetLast (fLast);
541
542 pix.SetHivsLoNbins(fHivsLoNbins);
543 pix.SetHivsLoFirst(fHivsLoFirst);
544 pix.SetHivsLoLast (fHivsLoLast);
545
546 InitHists(pix, fCam->GetAverageBadSector(j), j);
547
548 h = pix.GetHivsLo();
549
550 h->SetName (Form("H%sLoGainSector%02d",fHivsLoHistName.Data(),j));
551 h->SetTitle(Form("%s%s%02d",fHivsLoHistTitle.Data(),
552 " averaged on event-by-event basis Low Gain Area Sector ",j));
553 h->SetXTitle(fHivsLoHistXTitle.Data());
554 h->SetYTitle(fHivsLoHistYTitle.Data());
555 h->SetDirectory(0);
556 }
557 }
558}
559
560
561
562// -------------------------------------------------------------------------------
563//
564// Retrieves pointer to MExtractedSignalCam:
565//
566// Retrieves from MGeomCam:
567// - number of pixels
568// - number of pixel areas
569// - number of sectors
570//
571// Fills histograms (MHGausEvents::FillHistAndArray()) with:
572// - MExtractedSignalPix::GetExtractedSignalHiGain(pixid) / MExtractedSignalPix::GetExtractedSignalLoGain;
573// if the high-gain signal does not show high-gain saturation, but the low-gain
574// has been extracted.
575// - MArrivalTimePix::GetArrivalTimeHiGain(pixid) / MArrivalTimePix::GetArrivalTimeLoGain;
576// if the high-gain signal does not show high-gain saturation, but the low-gain
577// has been extracted.
578//
579Bool_t MHCalibrationHiLoCam::FillHists(const MParContainer *par, const Stat_t w)
580{
581
582 MExtractedSignalCam *signal = (MExtractedSignalCam*)par;
583 if (!signal)
584 {
585 gLog << err << "No argument in MExtractedSignal::Fill... abort." << endl;
586 return kFALSE;
587 }
588
589 const Int_t npixels = fGeom->GetNumPixels();
590 const Int_t nareas = fGeom->GetNumAreas();
591 const Int_t nsectors = fGeom->GetNumSectors();
592
593 fSumareahi .Reset();
594 fSumsectorhi.Reset();
595 fNumareahi .Reset();
596 fNumsectorhi.Reset();
597 fSumarealo .Reset();
598 fSumsectorlo.Reset();
599 fNumarealo .Reset();
600 fNumsectorlo.Reset();
601
602 for (Int_t i=0; i<npixels; i++)
603 {
604 const MExtractedSignalPix &pix = (*signal)[i];
605 const Int_t aidx = (*fGeom)[i].GetAidx();
606 const Int_t sector = (*fGeom)[i].GetSector();
607
608 const Float_t siglo = pix.GetExtractedSignalLoGain();
609
610 //
611 // Skip all pixels with:
612 // - Saturated high-gain
613 // - Not extracted low-gain
614 // (see MExtractTimeAndCharge::fLoGainSwitch for setting the criteria)
615 //
616 if (siglo < 0.5 || pix.IsHiGainSaturated())
617 continue;
618
619 const Float_t sighi = pix.GetExtractedSignalHiGain();
620
621 // *fLog << err << sighi << " " << siglo << endl;
622 const Float_t ratio = sighi / siglo;
623
624 MHCalibrationHiLoPix &histhi = (MHCalibrationHiLoPix&)(*this)[i];
625
626 histhi.FillHist(ratio);
627 histhi.FillHivsLo(sighi,siglo);
628
629 if (IsAverageing())
630 {
631 MHCalibrationHiLoPix &histhi2 = (MHCalibrationHiLoPix&)GetAverageHiGainArea(aidx);
632 histhi2.FillHivsLo(sighi,siglo);
633 }
634
635 fSumareahi [aidx] += ratio;
636 fNumareahi [aidx] ++;
637 fSumsectorhi[sector] += ratio;
638 fNumsectorhi[sector] ++;
639
640 if (IsLoGain())
641 {
642 const MArrivalTimePix &tix = (*fArrTimes)[i];
643 MHCalibrationPix &histlo = (*this)(i);
644
645 const Float_t diff = tix.GetArrivalTimeLoGain() - tix.GetArrivalTimeHiGain();
646
647 histlo.FillHist(diff);
648 fSumarealo [aidx] += diff;
649 fNumarealo [aidx] ++;
650 fSumsectorlo[sector] += diff;
651 fNumsectorlo[sector] ++;
652 }
653 }
654
655 if (!IsAverageing())
656 return kTRUE;
657
658 for (Int_t j=0; j<nareas; j++)
659 {
660
661 MHCalibrationHiLoPix &histhi = (MHCalibrationHiLoPix&)GetAverageHiGainArea(j);
662
663 if (IsOscillations())
664 histhi.FillHistAndArray(fNumareahi[j] == 0 ? 0. : fSumareahi[j]/fNumareahi[j]);
665 else
666 histhi.FillHist(fNumareahi[j] == 0 ? 0. : fSumareahi[j]/fNumareahi[j]);
667
668 if (IsLoGain())
669 {
670 MHCalibrationPix &histlo = GetAverageLoGainArea(j);
671 if (IsOscillations())
672 histlo.FillHistAndArray(fNumarealo[j] == 0 ? 0. : fSumarealo[j]/fNumarealo[j]);
673 else
674 histlo.FillHist(fNumarealo[j] == 0 ? 0. : fSumarealo[j]/fNumarealo[j]);
675 }
676 }
677
678 for (Int_t j=0; j<nsectors; j++)
679 {
680 MHCalibrationPix &hist = GetAverageHiGainSector(j);
681
682 if (IsOscillations())
683 hist.FillHistAndArray(fNumsectorhi[j] == 0 ? 0. : fSumsectorhi[j]/fNumsectorhi[j]);
684 else
685 hist.FillHist(fNumsectorhi[j] == 0 ? 0. : fSumsectorhi[j]/fNumsectorhi[j]);
686
687 if (IsLoGain())
688 {
689
690 MHCalibrationPix &histlo = GetAverageLoGainSector(j);
691
692 if (IsOscillations())
693 histlo.FillHistAndArray(fNumsectorlo[j] == 0 ? 0. : fSumsectorlo[j]/fNumsectorlo[j]);
694 else
695 histlo.FillHist(fNumsectorlo[j] == 0 ? 0. : fSumsectorlo[j]/fNumsectorlo[j]);
696 }
697 }
698
699 return kTRUE;
700}
701
702// --------------------------------------------------------------------------
703//
704// Calls:
705// - MHCalibrationCam::FitHiGainArrays() with flags:
706// MBadPixelsPix::kHiLoNotFitted and MBadPixelsPix::kHiLoOscillating
707// - MHCalibrationCam::FitLoGainArrays() with flags:
708// MBadPixelsPix::kHiLoNotFitted and MBadPixelsPix::kHiLoOscillating
709//
710Bool_t MHCalibrationHiLoCam::FinalizeHists()
711{
712
713 *fLog << endl;
714
715 MCalibrationCam *hilocam = fCam;
716
717 const Int_t nareas = fAverageHiGainAreas->GetSize();
718
719 for (Int_t i=0; i<fHiGainArray->GetSize(); i++)
720 {
721
722 MHCalibrationHiLoPix &hist = (MHCalibrationHiLoPix&)(*this)[i];
723
724 if (hist.IsExcluded())
725 continue;
726
727 CheckOverflow(hist);
728
729 TProfile *h = hist.GetHivsLo();
730 h->Fit("pol1","RQ","",fLowerFitLimitProfile,fUpperFitLimitProfile);
731
732 TF1 *fit = h->GetFunction("pol1");
733
734 const Float_t gainr = fit->GetParameter(1) > 0.001
735 ? 1./fit->GetParameter(1) : 0.;
736
737 // The number of used slices are just a mean value
738 // the real number might change from event to event.
739 // (up to 50%!)
740 const Float_t offset = fit->GetParameter(0)/fUsedLoGainSlices;
741
742 fHivsLoResults.Fill(offset,gainr);
743
744 MCalibrationHiLoPix &pix = (MCalibrationHiLoPix&)(*fCam)[i];
745 pix.SetOffsetPerSlice(offset);
746 pix.SetGainRatio (gainr );
747
748 }
749
750 //
751 // Check histogram overflow
752 //
753 if (IsAverageing())
754 {
755 for (Int_t j=0; j<nareas; j++)
756 CheckOverflow(GetAverageHiGainArea(j));
757
758 for (Int_t j=0; j<fAverageHiGainSectors->GetSize(); j++)
759 CheckOverflow(GetAverageHiGainSector(j));
760 }
761
762
763 FitHiGainArrays(*hilocam, *fBadPixels,
764 MBadPixelsPix::kHiLoNotFitted,
765 MBadPixelsPix::kHiLoOscillating);
766
767 if (!IsLoGain())
768 return kTRUE;
769
770 for (Int_t i=0; i<fLoGainArray->GetSize(); i++)
771 {
772
773 MHCalibrationPix &hist = (*this)(i);
774
775 if (hist.IsExcluded())
776 continue;
777
778 CheckOverflow(hist);
779 }
780
781 if (IsAverageing())
782 {
783 for (Int_t j=0; j<nareas; j++)
784 {
785
786 MHCalibrationHiLoPix &hist = (MHCalibrationHiLoPix&)GetAverageHiGainArea(j);
787 //
788 // Check histogram overflow
789 //
790 CheckOverflow(hist);
791
792 TProfile *h = hist.GetHivsLo();
793 h->Fit("pol1","RQ","",fLowerFitLimitProfile,fUpperFitLimitProfile);
794
795 TF1 *fit = h->GetFunction("pol1");
796
797 const Float_t gainr = fit->GetParameter(1) > 0.001
798 ? 1./fit->GetParameter(1) : 0.;
799
800 // The number of used slices are just a mean value
801 // the real number might change from event to event.
802 // (up to 50%!)
803 const Float_t offset = fit->GetParameter(0)/fUsedLoGainSlices;
804
805 MCalibrationHiLoPix &pix = (MCalibrationHiLoPix&)fCam->GetAverageArea(0);
806 pix.SetOffsetPerSlice(offset);
807 pix.SetGainRatio (gainr );
808
809 }
810
811 for (Int_t j=0; j<fAverageHiGainSectors->GetSize(); j++)
812 {
813
814 MHCalibrationHiLoPix &hist = (MHCalibrationHiLoPix&)GetAverageHiGainSector(j);
815 //
816 // Check histogram overflow
817 //
818 CheckOverflow(hist);
819
820 TProfile *h = hist.GetHivsLo();
821 h->Fit("pol1","RQ","",fLowerFitLimitProfile,fUpperFitLimitProfile);
822
823 TF1 *fit = h->GetFunction("pol1");
824
825 const Float_t gainr = fit->GetParameter(1) > 0.001
826 ? 1./fit->GetParameter(1) : 0.;
827
828 // The number of used slices are just a mean value
829 // the real number might change from event to event.
830 // (up to 50%!)
831 const Float_t offset = fit->GetParameter(0)/fUsedLoGainSlices;
832
833 MCalibrationHiLoPix &pix = (MCalibrationHiLoPix&)fCam->GetAverageSector(0);
834 pix.SetOffsetPerSlice(offset);
835 pix.SetGainRatio (gainr );
836
837 }
838 }
839
840 FitLoGainArrays(*hilocam, *fBadPixels,
841 MBadPixelsPix::kHiLoNotFitted,
842 MBadPixelsPix::kHiLoOscillating);
843
844 return kTRUE;
845}
846
847// --------------------------------------------------------------------------
848//
849// Sets all pixels to MBadPixelsPix::kUnreliableRun, if following flags are set:
850// - MBadPixelsPix::kHiLoNotFitted
851// - MBadPixelsPix::kHiLoOscillating
852//
853void MHCalibrationHiLoCam::FinalizeBadPixels()
854{
855
856 for (Int_t i=0; i<fBadPixels->GetSize(); i++)
857 {
858 MBadPixelsPix &bad = (*fBadPixels)[i];
859
860 if (bad.IsUncalibrated( MBadPixelsPix::kHiLoNotFitted ))
861 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun );
862
863 if (bad.IsUncalibrated( MBadPixelsPix::kHiLoOscillating))
864 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun );
865
866 }
867}
868
869// --------------------------------------------------------------------------
870//
871// The types are as follows:
872//
873// Fitted values:
874// ==============
875//
876// 0: Fitted Mean High-Gain Low-Gain Charge Ratio in FADC slices (MHGausEvents::GetMean()
877// 1: Error Mean High-Gain Low-Gain Charge Ratio in FADC slices (MHGausEvents::GetMeanErr()
878// 2: Sigma fitted High-Gain Low-Gain Charge Ratio in FADC slices (MHGausEvents::GetSigma()
879// 3: Error Sigma High-Gain Low-Gain Charge Ratio in FADC slices (MHGausEvents::GetSigmaErr()
880//
881// Useful variables derived from the fit results:
882// =============================================
883//
884// 4: Returned probability of Gauss fit (calls: MHGausEvents::GetProb())
885//
886// Localized defects:
887// ==================
888//
889// 5: Gaus fit not OK (calls: MHGausEvents::IsGausFitOK())
890// 6: Fourier spectrum not OK (calls: MHGausEvents::IsFourierSpectrumOK())
891//
892Bool_t MHCalibrationHiLoCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
893{
894
895 if (fHiGainArray->GetSize() <= idx)
896 return kFALSE;
897
898 const MHCalibrationPix &pixhi = (*this)[idx];
899 const MHCalibrationPix &pixlo = (*this)(idx);
900
901 switch (type)
902 {
903 case 0:
904 val = pixhi.GetMean();
905 break;
906 case 1:
907 val = pixhi.GetMeanErr();
908 break;
909 case 2:
910 val = pixhi.GetSigma();
911 break;
912 case 3:
913 val = pixhi.GetSigmaErr();
914 break;
915 case 4:
916 val = pixhi.GetProb();
917 break;
918 case 5:
919 if (!pixhi.IsGausFitOK())
920 val = 1.;
921 break;
922 case 6:
923 if (!pixhi.IsFourierSpectrumOK())
924 val = 1.;
925 break;
926 case 7:
927 if (!IsLoGain())
928 break;
929 val = pixlo.GetMean();
930 break;
931 case 8:
932 if (!IsLoGain())
933 break;
934 val = pixlo.GetMeanErr();
935 break;
936 case 9:
937 if (!IsLoGain())
938 break;
939 val = pixlo.GetSigma();
940 break;
941 case 10:
942 if (!IsLoGain())
943 break;
944 val = pixlo.GetSigmaErr();
945 break;
946 case 11:
947 if (!IsLoGain())
948 break;
949 val = pixlo.GetProb();
950 break;
951 case 12:
952 if (!IsLoGain())
953 break;
954 if (!pixlo.IsGausFitOK())
955 val = 1.;
956 break;
957 case 13:
958 if (!IsLoGain())
959 break;
960 if (!pixlo.IsFourierSpectrumOK())
961 val = 1.;
962 break;
963 default:
964 return kFALSE;
965 }
966 return kTRUE;
967}
968
969// --------------------------------------------------------------------------
970//
971// Calls MHCalibrationPix::DrawClone() for pixel idx
972//
973void MHCalibrationHiLoCam::DrawPixelContent(Int_t idx) const
974{
975 (*this)[idx].DrawClone();
976}
977
978void MHCalibrationHiLoCam::CheckOverflow( MHCalibrationPix &pix )
979{
980
981 if (pix.IsExcluded())
982 return;
983
984 TH1F *hist = pix.GetHGausHist();
985
986 Stat_t overflow = hist->GetBinContent(hist->GetNbinsX()+1);
987 if (overflow > fOverflowLimit*hist->GetEntries())
988 {
989 *fLog << warn << "Hist-overflow " << overflow
990 << " times in " << pix.GetName() << endl;
991 }
992
993 overflow = hist->GetBinContent(0);
994 if (overflow > fOverflowLimit*hist->GetEntries())
995 {
996 *fLog << warn << "Hist-underflow " << overflow
997 << " times in " << pix.GetName() << endl;
998 }
999}
1000
1001// -----------------------------------------------------------------------------
1002//
1003// Default draw:
1004//
1005// Displays the averaged areas, both amplification ratio as time difference
1006//
1007void MHCalibrationHiLoCam::Draw(const Option_t *opt)
1008{
1009
1010 if (!IsAverageing())
1011 return;
1012
1013 const Int_t nareas = fAverageHiGainAreas->GetSize();
1014 if (nareas == 0)
1015 return;
1016
1017 TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas(this);
1018 pad->SetBorderMode(0);
1019
1020 pad->Divide(IsLoGain() ? 2 : 1,nareas);
1021
1022 for (Int_t i=0; i<nareas;i++)
1023 {
1024
1025 pad->cd(IsLoGain() ? 2*(i+1)-1 : i+1);
1026
1027 GetAverageHiGainArea(i).Draw(opt);
1028
1029 if (IsLoGain())
1030 {
1031 pad->cd(2*(i+1));
1032
1033 TH1F *hist = GetAverageLoGainArea(i).GetHGausHist();
1034 hist->SetXTitle("Extracted Time Difference [FADC sl.]");
1035 GetAverageLoGainArea(i).Draw(opt);
1036 }
1037 }
1038}
1039
1040Int_t MHCalibrationHiLoCam::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
1041{
1042
1043 Bool_t rc = kFALSE;
1044
1045 if (MHCalibrationCam::ReadEnv(env,prefix,print))
1046 rc = kTRUE;
1047
1048 if (IsEnvDefined(env, prefix, "LowerFitLimitProfile", print))
1049 {
1050 SetLowerFitLimitProfile(GetEnvValue(env, prefix, "LowerFitLimitProfile", fLowerFitLimitProfile));
1051 rc = kTRUE;
1052 }
1053
1054 if (IsEnvDefined(env, prefix, "UpperFitLimitProfile", print))
1055 {
1056 SetUpperFitLimitProfile(GetEnvValue(env, prefix, "UpperFitLimitProfile", fUpperFitLimitProfile));
1057 rc = kTRUE;
1058 }
1059
1060 if (IsEnvDefined(env, prefix, "HivsLoNbins", print))
1061 {
1062 SetHivsLoNbins(GetEnvValue(env, prefix, "HivsLoNbins", fHivsLoNbins));
1063 rc = kTRUE;
1064 }
1065
1066 if (IsEnvDefined(env, prefix, "HivsLoFirst", print))
1067 {
1068 SetHivsLoFirst(GetEnvValue(env, prefix, "HivsLoFirst", fHivsLoFirst));
1069 rc = kTRUE;
1070 }
1071
1072 if (IsEnvDefined(env, prefix, "HivsLoLast", print))
1073 {
1074 SetHivsLoLast(GetEnvValue(env, prefix, "HivsLoLast", fHivsLoLast));
1075 rc = kTRUE;
1076 }
1077
1078 return rc;
1079}
Note: See TracBrowser for help on using the repository browser.