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

Last change on this file since 8177 was 8132, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 33.3 KB
Line 
1/* ======================================================================== *\
2! $Name: not supported by cvs2svn $:$Id: MHCalibrationHiLoCam.cc,v 1.20 2006-10-19 14:01:49 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) : 0.;
741
742 // The number of used slices are just a mean value
743 // the real number might change from event to event.
744 // (up to 50%!)
745 const Float_t offset = fit->GetParameter(0)/fUsedLoGainSlices;
746
747 fHivsLoResults.Fill(offset,gainr);
748
749 MCalibrationHiLoPix &pix = (MCalibrationHiLoPix&)(*fCam)[i];
750 pix.SetOffsetPerSlice(offset);
751 pix.SetGainRatio (gainr );
752
753 }
754
755 //
756 // Check histogram overflow
757 //
758 if (IsAverageing())
759 {
760 for (Int_t j=0; j<nareas; j++)
761 CheckOverflow(GetAverageHiGainArea(j));
762
763 for (Int_t j=0; j<fAverageHiGainSectors->GetSize(); j++)
764 CheckOverflow(GetAverageHiGainSector(j));
765 }
766
767
768 FitHiGainArrays(*hilocam,*badcam,
769 MBadPixelsPix::kHiLoNotFitted,
770 MBadPixelsPix::kHiLoOscillating);
771
772 if (!IsLoGain())
773 return kTRUE;
774
775 for (Int_t i=0; i<fLoGainArray->GetSize(); i++)
776 {
777
778 MHCalibrationPix &hist = (*this)(i);
779
780 if (hist.IsExcluded())
781 continue;
782
783 CheckOverflow(hist);
784 }
785
786 if (IsAverageing())
787 {
788 for (Int_t j=0; j<nareas; j++)
789 {
790
791 MHCalibrationHiLoPix &hist = (MHCalibrationHiLoPix&)GetAverageHiGainArea(j);
792 //
793 // Check histogram overflow
794 //
795 CheckOverflow(hist);
796
797 TProfile *h = hist.GetHivsLo();
798 h->Fit("pol1","RQ","",fLowerFitLimitProfile,fUpperFitLimitProfile);
799
800 TF1 *fit = h->GetFunction("pol1");
801
802 const Float_t gainr = fit->GetParameter(1) > 0.001
803 ? 1./fit->GetParameter(1) : 0.;
804
805 // The number of used slices are just a mean value
806 // the real number might change from event to event.
807 // (up to 50%!)
808 const Float_t offset = fit->GetParameter(0)/fUsedLoGainSlices;
809
810 MCalibrationHiLoPix &pix = (MCalibrationHiLoPix&)fCam->GetAverageArea(0);
811 pix.SetOffsetPerSlice(offset);
812 pix.SetGainRatio (gainr );
813
814 }
815
816 for (Int_t j=0; j<fAverageHiGainSectors->GetSize(); j++)
817 {
818
819 MHCalibrationHiLoPix &hist = (MHCalibrationHiLoPix&)GetAverageHiGainSector(j);
820 //
821 // Check histogram overflow
822 //
823 CheckOverflow(hist);
824
825 TProfile *h = hist.GetHivsLo();
826 h->Fit("pol1","RQ","",fLowerFitLimitProfile,fUpperFitLimitProfile);
827
828 TF1 *fit = h->GetFunction("pol1");
829
830 const Float_t gainr = fit->GetParameter(1) > 0.001
831 ? 1./fit->GetParameter(1) : 0.;
832
833 // The number of used slices are just a mean value
834 // the real number might change from event to event.
835 // (up to 50%!)
836 const Float_t offset = fit->GetParameter(0)/fUsedLoGainSlices;
837
838 MCalibrationHiLoPix &pix = (MCalibrationHiLoPix&)fCam->GetAverageSector(0);
839 pix.SetOffsetPerSlice(offset);
840 pix.SetGainRatio (gainr );
841
842 }
843 }
844
845 FitLoGainArrays(*hilocam,*badcam,
846 MBadPixelsPix::kHiLoNotFitted,
847 MBadPixelsPix::kHiLoOscillating);
848
849 return kTRUE;
850}
851
852// --------------------------------------------------------------------------
853//
854// Sets all pixels to MBadPixelsPix::kUnreliableRun, if following flags are set:
855// - MBadPixelsPix::kHiLoNotFitted
856// - MBadPixelsPix::kHiLoOscillating
857//
858void MHCalibrationHiLoCam::FinalizeBadPixels()
859{
860
861 MBadPixelsCam *badcam = fIntensBad ? fIntensBad->GetCam() : fBadPixels;
862
863 for (Int_t i=0; i<badcam->GetSize(); i++)
864 {
865 MBadPixelsPix &bad = (*badcam)[i];
866
867 if (bad.IsUncalibrated( MBadPixelsPix::kHiLoNotFitted ))
868 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun );
869
870 if (bad.IsUncalibrated( MBadPixelsPix::kHiLoOscillating))
871 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun );
872
873 }
874}
875
876// --------------------------------------------------------------------------
877//
878// The types are as follows:
879//
880// Fitted values:
881// ==============
882//
883// 0: Fitted Mean High-Gain Low-Gain Charge Ratio in FADC slices (MHGausEvents::GetMean()
884// 1: Error Mean High-Gain Low-Gain Charge Ratio in FADC slices (MHGausEvents::GetMeanErr()
885// 2: Sigma fitted High-Gain Low-Gain Charge Ratio in FADC slices (MHGausEvents::GetSigma()
886// 3: Error Sigma High-Gain Low-Gain Charge Ratio in FADC slices (MHGausEvents::GetSigmaErr()
887//
888// Useful variables derived from the fit results:
889// =============================================
890//
891// 4: Returned probability of Gauss fit (calls: MHGausEvents::GetProb())
892//
893// Localized defects:
894// ==================
895//
896// 5: Gaus fit not OK (calls: MHGausEvents::IsGausFitOK())
897// 6: Fourier spectrum not OK (calls: MHGausEvents::IsFourierSpectrumOK())
898//
899Bool_t MHCalibrationHiLoCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
900{
901
902 if (fHiGainArray->GetSize() <= idx)
903 return kFALSE;
904
905 const MHCalibrationPix &pixhi = (*this)[idx];
906 const MHCalibrationPix &pixlo = (*this)(idx);
907
908 switch (type)
909 {
910 case 0:
911 val = pixhi.GetMean();
912 break;
913 case 1:
914 val = pixhi.GetMeanErr();
915 break;
916 case 2:
917 val = pixhi.GetSigma();
918 break;
919 case 3:
920 val = pixhi.GetSigmaErr();
921 break;
922 case 4:
923 val = pixhi.GetProb();
924 break;
925 case 5:
926 if (!pixhi.IsGausFitOK())
927 val = 1.;
928 break;
929 case 6:
930 if (!pixhi.IsFourierSpectrumOK())
931 val = 1.;
932 break;
933 case 7:
934 if (!IsLoGain())
935 break;
936 val = pixlo.GetMean();
937 break;
938 case 8:
939 if (!IsLoGain())
940 break;
941 val = pixlo.GetMeanErr();
942 break;
943 case 9:
944 if (!IsLoGain())
945 break;
946 val = pixlo.GetSigma();
947 break;
948 case 10:
949 if (!IsLoGain())
950 break;
951 val = pixlo.GetSigmaErr();
952 break;
953 case 11:
954 if (!IsLoGain())
955 break;
956 val = pixlo.GetProb();
957 break;
958 case 12:
959 if (!IsLoGain())
960 break;
961 if (!pixlo.IsGausFitOK())
962 val = 1.;
963 break;
964 case 13:
965 if (!IsLoGain())
966 break;
967 if (!pixlo.IsFourierSpectrumOK())
968 val = 1.;
969 break;
970 default:
971 return kFALSE;
972 }
973 return kTRUE;
974}
975
976// --------------------------------------------------------------------------
977//
978// Calls MHCalibrationPix::DrawClone() for pixel idx
979//
980void MHCalibrationHiLoCam::DrawPixelContent(Int_t idx) const
981{
982 (*this)[idx].DrawClone();
983}
984
985void MHCalibrationHiLoCam::CheckOverflow( MHCalibrationPix &pix )
986{
987
988 if (pix.IsExcluded())
989 return;
990
991 TH1F *hist = pix.GetHGausHist();
992
993 Stat_t overflow = hist->GetBinContent(hist->GetNbinsX()+1);
994 if (overflow > fOverflowLimit*hist->GetEntries())
995 {
996 *fLog << warn << "Hist-overflow " << overflow
997 << " times in " << pix.GetName() << endl;
998 }
999
1000 overflow = hist->GetBinContent(0);
1001 if (overflow > fOverflowLimit*hist->GetEntries())
1002 {
1003 *fLog << warn << "Hist-underflow " << overflow
1004 << " times in " << pix.GetName() << endl;
1005 }
1006}
1007
1008// -----------------------------------------------------------------------------
1009//
1010// Default draw:
1011//
1012// Displays the averaged areas, both amplification ratio as time difference
1013//
1014void MHCalibrationHiLoCam::Draw(const Option_t *opt)
1015{
1016
1017 if (!IsAverageing())
1018 return;
1019
1020 const Int_t nareas = fAverageHiGainAreas->GetSize();
1021 if (nareas == 0)
1022 return;
1023
1024 TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas(this);
1025 pad->SetBorderMode(0);
1026
1027 pad->Divide(IsLoGain() ? 2 : 1,nareas);
1028
1029 for (Int_t i=0; i<nareas;i++)
1030 {
1031
1032 pad->cd(IsLoGain() ? 2*(i+1)-1 : i+1);
1033
1034 GetAverageHiGainArea(i).Draw(opt);
1035
1036 if (IsLoGain())
1037 {
1038 pad->cd(2*(i+1));
1039
1040 TH1F *hist = GetAverageLoGainArea(i).GetHGausHist();
1041 hist->SetXTitle("Extracted Time Difference [FADC sl.]");
1042 GetAverageLoGainArea(i).Draw(opt);
1043 }
1044 }
1045}
1046
1047Int_t MHCalibrationHiLoCam::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
1048{
1049
1050 Bool_t rc = kFALSE;
1051
1052 if (MHCalibrationCam::ReadEnv(env,prefix,print))
1053 rc = kTRUE;
1054
1055 if (IsEnvDefined(env, prefix, "LowerFitLimitProfile", print))
1056 {
1057 SetLowerFitLimitProfile(GetEnvValue(env, prefix, "LowerFitLimitProfile", fLowerFitLimitProfile));
1058 rc = kTRUE;
1059 }
1060
1061 if (IsEnvDefined(env, prefix, "UpperFitLimitProfile", print))
1062 {
1063 SetUpperFitLimitProfile(GetEnvValue(env, prefix, "UpperFitLimitProfile", fUpperFitLimitProfile));
1064 rc = kTRUE;
1065 }
1066
1067 if (IsEnvDefined(env, prefix, "HivsLoNbins", print))
1068 {
1069 SetHivsLoNbins(GetEnvValue(env, prefix, "HivsLoNbins", fHivsLoNbins));
1070 rc = kTRUE;
1071 }
1072
1073 if (IsEnvDefined(env, prefix, "HivsLoFirst", print))
1074 {
1075 SetHivsLoFirst(GetEnvValue(env, prefix, "HivsLoFirst", fHivsLoFirst));
1076 rc = kTRUE;
1077 }
1078
1079 if (IsEnvDefined(env, prefix, "HivsLoLast", print))
1080 {
1081 SetHivsLoLast(GetEnvValue(env, prefix, "HivsLoLast", fHivsLoLast));
1082 rc = kTRUE;
1083 }
1084
1085 return rc;
1086}
Note: See TracBrowser for help on using the repository browser.