source: trunk/MagicSoft/Mars/mhcalib/MHCalibrationHiLoCam.cc@ 8106

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