Context Navigation

source: trunk/MagicSoft/Mars/mcalib/MCalibrationCalc.cc@ 3029

Visit:

Last change on this file since 3029 was 3029, checked in by gaug, 21 years ago
* empty log message *
File size: 26.4 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 09/2003 <mailto:markus@ifae.es>
19	!
20	! Copyright: MAGIC Software Development, 2000-2001
21	!
22	!
23	\* ======================================================================== */
24
25	//////////////////////////////////////////////////////////////////////////////
26	//
27	// MCalibrationCalc
28	//
29	// Task to calculate the calibration conversion factors from the FADC
30	// time slices. The integrated time slices have to be delivered by an
31	// MExtractedSignalCam. The pedestals by an MPedestalCam and possibly
32	// arrival times from MArrivalTime
33	//
34	// The output container MCalibrationCam holds one entry of type MCalibrationPix
35	// for every pixel. It is filled in the following way:
36	//
37	// ProProcess: Search for MPedestalCam, MExtractedSignalCam
38	// Initialize MCalibrationCam
39	// Initialize MArrivalTime if exists
40	// Initialize pulser light wavelength
41	//
42	// ReInit: MCalibrationCam::InitSize(NumPixels) is called which allocates
43	// memory in a TClonesArray of type MCalibrationPix
44	// Initialize number of used FADC slices
45	// Optionally exclude pixels from calibration
46	//
47	// Process: Optionally, a cut on cosmics can be performed
48	//
49	// Every MCalibrationPix holds a histogram class,
50	// MHCalibrationPixel which itself hold histograms of type:
51	// HCharge(npix) (distribution of summed FADC time slice
52	// entries)
53	// HTime(npix) (distribution of position of maximum)
54	// HChargevsN(npix) (distribution of charges vs. event number.
55	//
56	// PostProcess: All histograms HCharge(npix) are fitted to a Gaussian
57	// All histograms HTime(npix) are fitted to a Gaussian
58	// The histogram HBlindPixelCharge (blind pixel) is fitted to
59	// a single PhE fit
60	//
61	// The histograms of the PIN Diode are fitted to Gaussians
62	//
63	// Fits can be excluded via the commands:
64	// MalibrationCam::SkipTimeFits() (skip all time fits)
65	// MalibrationCam::SkipBlindPixelFits() (skip all blind
66	// pixel fits)
67	// MalibrationCam::SkipPinDiodeFits() (skip all PIN Diode
68	// fits)
69	//
70	// Hi-Gain vs. Lo-Gain Calibration (very memory-intensive)
71	// can be skipped with the command:
72	// MalibrationCam::SkipHiLoGainCalibration()
73	//
74	// Input Containers:
75	// MRawEvtData
76	// MPedestalCam
77	// MExtractedSignalCam
78	//
79	// Output Containers:
80	// MCalibrationCam
81	//
82	//////////////////////////////////////////////////////////////////////////////
83	#include "MCalibrationCalc.h"
84
85	// FXIME: Usage of fstream is a preliminary workaround!
86	#include <fstream>
87
88	// FXIME: This has to be removed!!!! (YES, WHEN WE HAVE ACCESS TO THE DATABASE!!!!!)
89	#include "MCalibrationConfig.h"
90
91	#include <TSystem.h>
92
93	#include "MLog.h"
94	#include "MLogManip.h"
95
96	#include "MParList.h"
97
98	#include "MGeomCam.h"
99	#include "MRawRunHeader.h"
100	#include "MRawEvtPixelIter.h"
101
102	#include "MPedestalCam.h"
103	#include "MPedestalPix.h"
104
105	#include "MCalibrationCam.h"
106	#include "MCalibrationPix.h"
107
108	#include "MExtractedSignalCam.h"
109	#include "MExtractedSignalPix.h"
110
111	#include "MCalibrationBlindPix.h"
112	#include "MCalibrationPINDiode.h"
113
114	#include "MArrivalTime.h"
115
116	ClassImp(MCalibrationCalc);
117
118	using namespace std;
119
120	const Int_t MCalibrationCalc::fBlindPixelId = 559;
121	const Int_t MCalibrationCalc::fPINDiodeId = 9999;
122	const Byte_t MCalibrationCalc::fgSaturationLimit = 254;
123	const Byte_t MCalibrationCalc::fgBlindPixelFirst = 3;
124	const Byte_t MCalibrationCalc::fgBlindPixelLast = 12;
125
126	// --------------------------------------------------------------------------
127	//
128	// Default constructor.
129	//
130	MCalibrationCalc::MCalibrationCalc(const char name, const char title)
131	: fPedestals(NULL), fCalibrations(NULL), fSignals(NULL),
132	fRawEvt(NULL), fRunHeader(NULL), fArrivalTime(NULL), fEvtTime(NULL)
133	{
134
135	fName = name ? name : "MCalibrationCalc";
136	fTitle = title ? title : "Task to calculate the calibration constants and MCalibrationCam ";
137
138	AddToBranchList("MRawEvtData.fHiGainPixId");
139	AddToBranchList("MRawEvtData.fLoGainPixId");
140	AddToBranchList("MRawEvtData.fHiGainFadcSamples");
141	AddToBranchList("MRawEvtData.fLoGainFadcSamples");
142
143	Clear();
144	}
145
146	void MCalibrationCalc::Clear(const Option_t *o)
147	{
148
149	SETBIT(fFlags, kUseTimes);
150	SETBIT(fFlags, kUseBlindPixelFit);
151	SETBIT(fFlags, kUseCosmicsRejection);
152	SETBIT(fFlags, kUseQualityChecks);
153	SETBIT(fFlags, kHiLoGainCalibration);
154
155	CLRBIT(fFlags, kBlindPixelOverFlow);
156	CLRBIT(fFlags, kPINDiodeOverFlow);
157	CLRBIT(fFlags, kHiGainOverFlow);
158	CLRBIT(fFlags, kLoGainOverFlow);
159	// As long as we don't have the PIN Diode:
160	CLRBIT(fFlags, kUsePinDiodeFit);
161
162
163
164	fEvents = 0;
165	fCosmics = 0;
166
167	fNumHiGainSamples = 0;
168	fNumLoGainSamples = 0;
169	fConversionHiLo = 0;
170	fNumExcludedPixels = 0;
171
172	fColor = kECT1;
173	}
174
175
176	MCalibrationBlindPix *MCalibrationCalc::GetBlindPixel() const
177	{
178	return fCalibrations->GetBlindPixel();
179	}
180
181	MCalibrationPINDiode *MCalibrationCalc::GetPINDiode() const
182	{
183	return fCalibrations->GetPINDiode();
184	}
185
186	// --------------------------------------------------------------------------
187	//
188	// The PreProcess searches for the following input containers:
189	// - MRawEvtData
190	// - MPedestalCam
191	//
192	// The following output containers are also searched and created if
193	// they were not found:
194	//
195	// - MHCalibrationBlindPixel
196	// - MCalibrationCam
197	//
198	// The following output containers are only searched, but not created
199	//
200	// - MArrivaltime
201	// - MTime
202	//
203	Int_t MCalibrationCalc::PreProcess(MParList *pList)
204	{
205
206	fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
207	if (!fRawEvt)
208	{
209	*fLog << err << dbginf << "MRawEvtData not found... aborting." << endl;
210	return kFALSE;
211	}
212
213	const MRawRunHeader runheader = (MRawRunHeader)pList->FindObject("MRawRunHeader");
214	if (!runheader)
215	*fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl;
216	else
217	if (runheader->GetRunType() == kRTMonteCarlo)
218	{
219	return kTRUE;
220	}
221
222	fCalibrations = (MCalibrationCam*)pList->FindCreateObj("MCalibrationCam");
223	if (!fCalibrations)
224	{
225	*fLog << err << dbginf << "MCalibrationCam could not be created ... aborting." << endl;
226	return kFALSE;
227	}
228
229	fArrivalTime = (MArrivalTime*)pList->FindObject("MArrivalTime");
230
231	fEvtTime = (MTime*)pList->FindObject("MTime");
232
233	switch (fColor)
234	{
235	case kEBlue:
236	fCalibrations->SetColor(MCalibrationCam::kECBlue);
237	break;
238	case kEGreen:
239	fCalibrations->SetColor(MCalibrationCam::kECGreen);
240	break;
241	case kEUV:
242	fCalibrations->SetColor(MCalibrationCam::kECUV);
243	break;
244	case kECT1:
245	fCalibrations->SetColor(MCalibrationCam::kECCT1);
246	break;
247	default:
248	fCalibrations->SetColor(MCalibrationCam::kECCT1);
249	}
250
251	fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
252	if (!fPedestals)
253	{
254	*fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl;
255	return kFALSE;
256	}
257
258
259	fSignals = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
260	if (!fSignals)
261	{
262	*fLog << err << dbginf << "Cannot find MExtractedSignalCam ... aborting" << endl;
263	return kFALSE;
264	}
265
266	return kTRUE;
267	}
268
269
270	// --------------------------------------------------------------------------
271	//
272	// The ReInit searches for the following input containers:
273	// - MRawRunHeader
274	//
275	Bool_t MCalibrationCalc::ReInit(MParList *pList )
276	{
277
278	fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
279	if (!fRunHeader)
280	{
281	*fLog << err << dbginf << ": MRawRunHeader not found... aborting." << endl;
282	return kFALSE;
283	}
284
285
286	MGeomCam cam = (MGeomCam)pList->FindObject("MGeomCam");
287	if (!cam)
288	{
289	*fLog << err << GetDescriptor() << ": No MGeomCam found... aborting." << endl;
290	return kFALSE;
291	}
292
293	fNumHiGainSamples = fSignals->GetNumUsedHiGainFADCSlices();
294	fNumLoGainSamples = fSignals->GetNumUsedLoGainFADCSlices();
295	fSqrtHiGainSamples = TMath::Sqrt((Float_t)fNumHiGainSamples);
296
297	UInt_t npixels = cam->GetNumPixels();
298
299	for (UInt_t i=0; i<npixels; i++)
300	{
301
302	MCalibrationPix &pix = (*fCalibrations)[i];
303	pix.DefinePixId(i);
304
305	pix.SetAbsTimeBordersHiGain(fSignals->GetFirstUsedSliceHiGain(),
306	fSignals->GetLastUsedSliceHiGain());
307	pix.SetAbsTimeBordersLoGain(fSignals->GetFirstUsedSliceLoGain(),
308	fSignals->GetLastUsedSliceLoGain());
309
310	if (!TESTBIT(fFlags,kUseQualityChecks))
311	pix.SetExcludeQualityCheck();
312
313	}
314
315	//
316	// Look for file to exclude pixels from analysis
317	//
318	if (!fExcludedPixelsFile.IsNull())
319	{
320
321	fExcludedPixelsFile = gSystem->ExpandPathName(fExcludedPixelsFile.Data());
322
323	//
324	// Initialize reading the file
325	//
326	ifstream in(fExcludedPixelsFile.Data(),ios::in);
327
328	if (in)
329	{
330	*fLog << inf << "Use excluded pixels from file: '" << fExcludedPixelsFile.Data() << "'" << endl;
331	//
332	// Read the file and count the number of entries
333	//
334	UInt_t pixel = 0;
335
336	while (++fNumExcludedPixels)
337	{
338
339	in >> pixel;
340
341	if (!in.good())
342	break;
343	//
344	// Check for out of range
345	//
346	if (pixel > npixels)
347	{
348	*fLog << warn << "WARNING: To be excluded pixel: " << pixel
349	<< " is out of range " << endl;
350	continue;
351	}
352	//
353	// Exclude pixel
354	//
355	MCalibrationPix &pix = (*fCalibrations)[pixel];
356	pix.SetExcluded();
357
358	*fLog << GetDescriptor() << inf << ": Exclude Pixel: " << pixel << endl;
359
360	}
361
362	if (--fNumExcludedPixels == 0)
363	*fLog << warn << "WARNING: File '" << fExcludedPixelsFile.Data()
364	<< "'" << " is empty " << endl;
365	else
366	fCalibrations->SetNumPixelsExcluded(fNumExcludedPixels);
367
368	}
369	else
370	*fLog << warn << dbginf << "Cannot open file '" << fExcludedPixelsFile.Data() << "'" << endl;
371	}
372
373	return kTRUE;
374	}
375
376
377	// --------------------------------------------------------------------------
378	//
379	// Calculate the integral of the FADC time slices and store them as a new
380	// pixel in the MCerPhotEvt container.
381	//
382	Int_t MCalibrationCalc::Process()
383	{
384
385	//
386	// Initialize pointers to blind pixel, PIN Diode and individual pixels
387	//
388	MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
389	MCalibrationPINDiode &pindiode = *(fCalibrations->GetPINDiode());
390
391	MRawEvtPixelIter pixel(fRawEvt);
392
393	//
394	// Perform cosmics cut
395	//
396	if (TESTBIT(fFlags,kUseCosmicsRejection))
397	{
398
399	Int_t cosmicpix = 0;
400
401	//
402	// Create a first loop to sort out the cosmics ...
403	//
404	// This is a very primitive check for the number of cosmicpixs
405	// The cut will be applied in the fit, but for the blind pixel,
406	// we need to remove this event
407	//
408	// FIXME: In the future need a much more sophisticated one!!!
409	//
410
411	while (pixel.Next())
412	{
413
414	const UInt_t pixid = pixel.GetPixelId();
415
416	MExtractedSignalPix &sig = (*fSignals)[pixid];
417	MPedestalPix &ped = (*fPedestals)[pixid];
418	Float_t pedrms = ped.GetPedestalRms()*fSqrtHiGainSamples;
419	Float_t sumhi = sig.GetExtractedSignalHiGain();
420
421	//
422	// We consider a pixel as presumably due to cosmics
423	// if its sum of FADC slices is lower than 3 pedestal RMS
424	//
425	if (sumhi < 3.*pedrms )
426	cosmicpix++;
427	}
428
429	//
430	// If the camera contains more than 230
431	// (this is the number of outer pixels plus about 50 inner ones)
432	// presumed pixels due to cosmics, then the event is discarted.
433	// This procedure is more or less equivalent to keeping only events
434	// with at least 350 pixels with high signals.
435	//
436	if (cosmicpix > 230.)
437	{
438	fCosmics++;
439	return kCONTINUE;
440	}
441
442	pixel.Reset();
443	}
444
445	fEvents++;
446
447	Float_t referencetime = 0.;
448
449	//
450	// Create a (second) loop to do fill the calibration histograms
451	// Search for: a signal in MExtractedSignalCam
452	// a time in MArrivalTime.
453	// Fill histograms with:
454	// charge
455	// charge vs. event nr.
456	// relative arrival time w.r.t. pixel 1
457	//
458	while (pixel.Next())
459	{
460
461	const UInt_t pixid = pixel.GetPixelId();
462
463	MCalibrationPix &pix = (*fCalibrations)[pixid];
464
465	if (pix.IsExcluded())
466	continue;
467
468	MExtractedSignalPix &sig = (*fSignals)[pixid];
469
470	const Float_t sumhi = sig.GetExtractedSignalHiGain();
471	const Float_t sumlo = sig.GetExtractedSignalLoGain();
472
473	Float_t abstime = 0.;
474	Float_t reltime = 0.;
475
476	if (TESTBIT(fFlags,kUseTimes))
477	{
478
479	//
480	// Have a look in MArrivalTime,
481	// otherwise search the position of maximum bin
482	// in MRawEvtData
483	//
484	if (fArrivalTime)
485	{
486	abstime = (*fArrivalTime)[pixid];
487	reltime = abstime - (*fArrivalTime)[1];
488	}
489
490	else
491	{
492	if (pixid == 1)
493	referencetime = (Float_t)pixel.GetIdxMaxHiGainSample();
494
495	if (sig.IsLoGainUsed())
496	{
497	abstime = (Float_t)pixel.GetIdxMaxLoGainSample();
498	reltime = abstime - referencetime;
499	}
500	else
501	{
502	abstime = (Float_t)pixel.GetIdxMaxHiGainSample();
503	reltime = abstime - referencetime;
504	}
505	}
506	} /* if Use Times */
507
508	switch(pixid)
509	{
510
511	case fBlindPixelId:
512
513	if (TESTBIT(fFlags,kUseBlindPixelFit))
514	{
515
516	Float_t blindpixelsum = 0.;
517	//
518	// We need a dedicated signal extractor for the blind pixel
519	//
520	MPedestalPix &ped = (*fPedestals)[pixid];
521	if (!CalcSignalBlindPixel(pixel.GetHiGainSamples(), blindpixelsum, ped.GetPedestal()))
522	return kFALSE;
523
524	if (!blindpixel.FillCharge(blindpixelsum))
525	*fLog << warn <<
526	"Overflow or Underflow occurred filling Blind Pixel sum = " << blindpixelsum << endl;
527
528	if (TESTBIT(fFlags,kUseTimes))
529	{
530	if (!blindpixel.FillTime(reltime))
531	*fLog << warn <<
532	"Overflow or Underflow occurred filling Blind Pixel time = " << reltime << endl;
533	}
534
535	if (!TESTBIT(fFlags,kBlindPixelOverFlow))
536	if (!blindpixel.FillRChargevsTime(blindpixelsum,fEvents))
537	{
538	*fLog << warn <<
539	"Overflow or Underflow occurred filling Blind Pixel eventnr = " << fEvents << endl;
540	SETBIT(fFlags,kBlindPixelOverFlow);
541	}
542	} /* if use blind pixel */
543
544	break;
545
546	case fPINDiodeId:
547
548	if (TESTBIT(fFlags,kUsePinDiodeFit))
549	{
550
551	if (!pindiode.FillCharge(sumhi))
552	*fLog << warn <<
553	"Overflow or Underflow occurred filling PINDiode: sum = " << sumhi << endl;
554
555	if (TESTBIT(fFlags,kUseTimes))
556	{
557	if (!pindiode.FillAbsTime(abstime))
558	*fLog << warn <<
559	"Overflow or Underflow occurred filling PINDiode abs. time = " << abstime << endl;
560	if (!pindiode.FillRelTime(reltime))
561	*fLog << warn <<
562	"Overflow or Underflow occurred filling PINDiode rel. time = " << reltime << endl;
563	}
564
565	if (!TESTBIT(fFlags,kPINDiodeOverFlow))
566	if (!pindiode.FillRChargevsTime(sumhi,fEvents))
567	{
568	*fLog << warn
569	<< "Overflow or Underflow occurred filling PINDiode: eventnr = "
570	<< fEvents << endl;
571	SETBIT(fFlags,kPINDiodeOverFlow);
572	}
573
574	} /* if use PIN Diode */
575
576	break;
577
578	default:
579
580	if (!TESTBIT(fFlags,kLoGainOverFlow))
581	if (!pix.FillRChargevsTimeLoGain(sumlo,fEvents))
582	{
583	*fLog << warn
584	<< "Overflow filling Histogram ChargevsNLoGain eventnr = "
585	<< fEvents << endl;
586	SETBIT(fFlags,kLoGainOverFlow);
587	SkipHiLoGainCalibration();
588	}
589
590	if (!TESTBIT(fFlags,kHiGainOverFlow))
591	if (!pix.FillRChargevsTimeHiGain(sumhi,fEvents))
592	{
593	*fLog << warn
594	<< "Overflow filling Histogram ChargevsNHiGain eventnr = "
595	<< fEvents << endl;
596	SETBIT(fFlags,kHiGainOverFlow);
597	SkipHiLoGainCalibration();
598	}
599
600
601	if (TESTBIT(fFlags,kHiLoGainCalibration))
602	pix.FillChargesInGraph(sumhi,sumlo);
603
604	if (sig.IsLoGainUsed())
605	{
606
607	if (!pix.FillChargeLoGain(sumlo))
608	*fLog << warn << "Could not fill Lo Gain Charge of pixel: " << pixid
609	<< " signal = " << sumlo << endl;
610
611	if (TESTBIT(fFlags,kUseTimes))
612	{
613	if (!pix.FillAbsTimeLoGain(abstime))
614	*fLog << warn << "Could not fill Lo Gain Abs. Time of pixel: "
615	<< pixid << " time = " << abstime << endl;
616	if (!pix.FillRelTimeLoGain(reltime))
617	*fLog << warn << "Could not fill Lo Gain Rel. Time of pixel: "
618	<< pixid << " time = " << reltime << endl;
619	}
620	} /* if (sig.IsLoGainUsed()) */
621	else
622	{
623	if (!pix.FillChargeHiGain(sumhi))
624	*fLog << warn << "Could not fill Hi Gain Charge of pixel: " << pixid
625	<< " signal = " << sumhi << endl;
626
627	if (TESTBIT(fFlags,kUseTimes))
628	{
629	if (!pix.FillAbsTimeHiGain(abstime))
630	*fLog << warn << "Could not fill Hi Gain Abs. Time of pixel: "
631	<< pixid << " time = " << abstime << endl;
632	if (!pix.FillRelTimeHiGain(reltime))
633	*fLog << warn << "Could not fill Hi Gain Rel. Time of pixel: "
634	<< pixid << " time = " << reltime << endl;
635	}
636	} /* else (sig.IsLoGainUsed()) */
637	break;
638
639	} /* switch(pixid) */
640
641	} /* while (pixel.Next()) */
642
643	return kTRUE;
644	}
645
646	Int_t MCalibrationCalc::PostProcess()
647	{
648
649	*fLog << inf << endl;
650
651	if (fEvents == 0)
652	{
653	*fLog << err << GetDescriptor()
654	<< ": This run contains only cosmics or pedestals, "
655	<< "cannot find events with more than 350 illuminated pixels. " << endl;
656	return kFALSE;
657	}
658
659	if (fEvents < fCosmics)
660	*fLog << warn << GetDescriptor()
661	<< ": WARNING: Run contains more cosmics or pedestals than calibration events " << endl;
662
663
664	*fLog << inf << GetDescriptor() << ": Cut Histogram Edges" << endl;
665
666	//
667	// Cut edges to make fits and viewing of the hists easier
668	//
669	fCalibrations->CutEdges();
670
671	//
672	// Fit the blind pixel
673	//
674	if (TESTBIT(fFlags,kUseBlindPixelFit))
675	{
676	//
677	// Get pointer to blind pixel
678	//
679	MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
680
681	*fLog << inf << GetDescriptor() << ": Fitting the Blind Pixel" << endl;
682
683	//
684	// retrieve mean and sigma of the blind pixel pedestal,
685	// so that we can use it for the fit
686	//
687	if (fPedestals->GetHistSize() > fBlindPixelId)
688	{
689	//
690	// retrieve the pedestal pix of the blind pixel
691	//
692	MHPedestalPixel &pedhist = (*fPedestals)(fBlindPixelId);
693	MPedestalPix &pedpix = (*fPedestals)[fBlindPixelId];
694	//
695	// retrieve the histogram containers
696	//
697	MHCalibrationBlindPixel *hist = blindpixel.GetHist();
698	//
699	// Set the corresponding values
700	//
701	const Float_t nslices = (Float_t)(fgBlindPixelLast-fgBlindPixelFirst+1);
702	const Float_t sqrslice = TMath::Sqrt(nslices);
703	const ULong_t nentries = fPedestals->GetTotalEntries();
704
705	const Float_t peddiff = (pedhist.GetChargeMean()-pedpix.GetPedestal())*nslices;
706
707	Float_t pederr = pedhist.GetChargeMeanErr()*pedhist.GetChargeMeanErr();
708	pederr += pedpix.GetPedestalRms()*pedpix.GetPedestalRms()/nentries;
709	pederr = TMath::Sqrt(pederr)*sqrslice;
710
711	//
712	// Fitted sigma: 1. one sqrt(Nr. slices) for the division which is not
713	// not appropriate: sigma(real)/slice = GetSigma*sqrt(nslices)
714	// 2. another sqrt(Nr. slices) to calculate back to number
715	// of slices
716	//
717	const Float_t pedsigma = pedhist.GetChargeSigma()*nslices;
718	const Float_t pedsigmaerr = pedhist.GetChargeSigmaErr()*nslices;
719
720	hist->SetMeanPedestal(peddiff);
721	hist->SetMeanPedestalErr(pederr);
722	hist->SetSigmaPedestal(pedsigma);
723	hist->SetSigmaPedestalErr(pedsigmaerr);
724	}
725
726	if (!blindpixel.FitCharge())
727	{
728	*fLog << warn << "Could not fit the blind pixel! " << endl;
729	*fLog << warn << "Setting bit kBlindPixelMethodValid to FALSE in MCalibrationCam" << endl;
730	fCalibrations->SetBlindPixelMethodValid(kFALSE);
731	}
732
733	fCalibrations->SetBlindPixelMethodValid(kTRUE);
734	blindpixel.DrawClone();
735	}
736	else
737	*fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Fit " << endl;
738
739
740	*fLog << inf << GetDescriptor() << ": Fitting the Normal Pixels" << endl;
741
742	//
743	// loop over the pedestal events and check if we have calibration
744	//
745	for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
746	{
747
748	MCalibrationPix &pix = (*fCalibrations)[pixid];
749
750	//
751	// Check if the pixel has been excluded from the fits
752	//
753	if (pix.IsExcluded())
754	continue;
755
756	//
757	// get the pedestals
758	//
759	const Float_t ped = (*fPedestals)[pixid].GetPedestal();
760	const Float_t prms = (*fPedestals)[pixid].GetPedestalRms();
761
762	//
763	// set them in the calibration camera
764	//
765	pix.SetPedestal(ped,prms,(Float_t)fNumHiGainSamples,(Float_t)fNumLoGainSamples);
766
767	//
768	// perform the Gauss fits to the charges
769	//
770	pix.FitCharge();
771
772	//
773	// Perform the Gauss fits to the arrival times
774	//
775	if (TESTBIT(fFlags,kUseTimes))
776	pix.FitTime();
777
778	}
779
780	if (TESTBIT(fFlags,kUseBlindPixelFit) && fCalibrations->IsBlindPixelMethodValid())
781	{
782
783	if (!fCalibrations->CalcNumPhotInsidePlexiglass())
784	{
785	*fLog << err
786	<< "Could not calculate the number of photons from the blind pixel " << endl;
787	*fLog << err
788	<< "You can try to calibrate using the MCalibrationCalc::SkipBlindPixelFit()" << endl;
789	fCalibrations->SetBlindPixelMethodValid(kFALSE);
790	}
791	}
792	else
793	*fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Calibration! " << endl;
794
795
796	if (TESTBIT(fFlags,kUsePinDiodeFit) && fCalibrations->IsPINDiodeMethodValid())
797	{
798
799	if (!fCalibrations->CalcNumPhotInsidePlexiglass())
800	{
801	*fLog << err
802	<< "Could not calculate the number of photons from the blind pixel " << endl;
803	*fLog << err
804	<< "You can try to calibrate using the MCalibrationCalc::SkipPINDiodeFit()" << endl;
805	fCalibrations->SetPINDiodeMethodValid(kFALSE);
806	}
807	}
808	else
809	*fLog << inf << GetDescriptor() << ": Skipping PIN Diode Calibration! " << endl;
810
811	fCalibrations->SetReadyToSave();
812
813	if (GetNumExecutions()==0)
814	return kTRUE;
815
816	*fLog << inf << endl;
817	*fLog << dec << setfill(' ') << fCosmics << " Events presumably cosmics" << endl;
818
819	return kTRUE;
820	}
821
822	Bool_t MCalibrationCalc::CalcSignalBlindPixel(Byte_t *ptr, Float_t &signal, const Float_t ped) const
823	{
824
825	Byte_t *start = ptr + fgBlindPixelFirst;
826	Byte_t *end = ptr + fgBlindPixelLast;
827
828	Byte_t sum = 0;
829	Int_t sat = 0;
830
831	ptr = start;
832
833	while (ptr<end)
834	{
835	sum += *ptr;
836
837	if (*ptr++ >= fgSaturationLimit)
838	sat++;
839	}
840
841	if (sat)
842	{
843	*fLog << err << "HI Gain Saturation occurred in the blind pixel! "
844	<< " Do not know yet how to treat this ... aborting " << endl;
845	return kFALSE;
846	}
847
848	signal = (Float_t)sum - ped*(Float_t)(fgBlindPixelLast-fgBlindPixelFirst+1);
849
850	return kTRUE;
851	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: