source: trunk/MagicSoft/Mars/manalysis/MCalibrationCalc.cc@ 2627

Last change on this file since 2627 was 2627, checked in by gaug, 21 years ago
*** empty log message ***
File size: 15.2 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// This is a task which calculates the number of photons from the FADC //
30// time slices. At the moment it integrates simply the FADC values. //
31// //
32// Input Containers: //
33// MRawEvtData //
34// //
35// Output Containers: //
36// MCalibrationCam //
37// //
38// //
39// The class MCalibrationCam hold one entry of type MCalibrationPix for //
40// every pixel. It is filled in the following way: //
41// PreParocess: MalibrationCam::InitSize(577) is called which allocates //
42// memory in an TClonesArray of type MCalibrationPix and //
43// all pointers to NULL. //
44// //
45// Process: The NULL pointer is tested on every pixel via //
46// MalibrationCam::IsPixelUsed(npix). //
47// //
48// In case, IsPixelUsed returns NULL, //
49// MalibrationCam::AddPixel(npix) is invoked which creates a //
50// new MCalibrationPix(npix) in the npix's entry //
51// of the TClonesArray. //
52// //
53// Every MCalibrationPix holds a histogram class, //
54// MHCalibrationPixel which itself hold histograms of type: //
55// HCharge(npix) (distribution of summed FADC time slice entries)
56// HTime(npix) (distribution of position of maximum)
57// HChargevsN(npix) (distribution of charges vs. event number.
58//
59// PostProcess: All histograms HCharge(npix) are fitted to a Gaussian
60// All histograms HTime(npix) are fitted to a Gaussian
61// The histogram HBlindPixelCharge (blind pixel) is fitted to a single
62// PhE fit
63// The histogram HBlindPixelTime (blind pixel) is fitted to a Gaussian
64// The histograms of the PIN Diode are fitted to Gaussians
65//
66// Fits can be excluded via the commands:
67// MalibrationCam::SetSkipTimeFits() (skip all time fits)
68// MalibrationCam::SetSkipBlindPixelFits() (skip all blind pixel fits)
69// MalibrationCam::SetSkipPinDiodeFits() (skip all PIN Diode fits)
70//
71//////////////////////////////////////////////////////////////////////////////
72
73#include "MCalibrationCalc.h"
74#include "MCalibrationConfig.h"
75#include "MCalibrationFits.h"
76
77#include "MCalibrationCam.h"
78#include "MCalibrationPix.h"
79#include "MCalibrationBlindPix.h"
80#include "MCalibrationPINDiode.h"
81
82#include "MPedestalCam.h"
83#include "MPedestalPix.h"
84
85#include "MLog.h"
86#include "MLogManip.h"
87
88#include "MParList.h"
89#include "MH.h"
90
91#include "MRawRunHeader.h"
92#include "MRawEvtData.h" // MRawEvtData::GetNumPixels
93#include "MRawEvtPixelIter.h"
94
95#include "MTime.h"
96#include "TMath.h"
97
98ClassImp(MCalibrationCalc);
99
100using namespace std;
101// --------------------------------------------------------------------------
102//
103// Default constructor. b is the number of slices before the maximum slice,
104// a the number of slices behind the maximum slice which is taken as signal.
105//
106MCalibrationCalc::MCalibrationCalc(const char *name, const char *title)
107 : fColor(kEBlue)
108{
109
110 fName = name ? name : "MCalibrationCalc";
111 fTitle = title ? title : "Task to calculate the calibration constants and MCalibrationCam ";
112
113 AddToBranchList("MRawEvtData.fHiGainPixId");
114 AddToBranchList("MRawEvtData.fLoGainPixId");
115 AddToBranchList("MRawEvtData.fHiGainFadcSamples");
116 AddToBranchList("MRawEvtData.fLoGainFadcSamples");
117
118 SETBIT(fFlags, kUseTimeFits);
119 SETBIT(fFlags, kUseBlindPixelFit);
120 SETBIT(fFlags, kUsePinDiodeFit);
121}
122
123// --------------------------------------------------------------------------
124//
125// The PreProcess searches for the following input containers:
126// - MRawEvtData
127// - MPedestalCam
128//
129// The following output containers are also searched and created if
130// they were not found:
131//
132// - MHCalibrationBlindPixel
133// - MCalibrationCam
134// - MTime
135//
136Int_t MCalibrationCalc::PreProcess(MParList *pList)
137{
138
139 fHistOverFlow = 0;
140 fEvents = 0;
141 fCosmics = 0;
142
143 fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
144 if (!fRawEvt)
145 {
146 *fLog << dbginf << "MRawEvtData not found... aborting." << endl;
147 return kFALSE;
148 }
149
150 const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
151 if (!runheader)
152 *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl;
153 else
154 if (runheader->GetRunType() == kRTMonteCarlo)
155 {
156 return kTRUE;
157 }
158
159 fCalibrations = (MCalibrationCam*)pList->FindCreateObj("MCalibrationCam");
160 if (!fCalibrations)
161 {
162 *fLog << err << dbginf << "MCalibrationCam could not be created ... aborting." << endl;
163 return kFALSE;
164 }
165
166
167 switch (fColor)
168 {
169 case kEBlue:
170 fCalibrations->SetColor(MCalibrationCam::kECBlue);
171 break;
172 case kEGreen:
173 fCalibrations->SetColor(MCalibrationCam::kECGreen);
174 break;
175 case kEUV:
176 fCalibrations->SetColor(MCalibrationCam::kECUV);
177 break;
178 case kECT1:
179 fCalibrations->SetColor(MCalibrationCam::kECCT1);
180 break;
181 default:
182 fCalibrations->SetColor(MCalibrationCam::kECCT1);
183 }
184
185 fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
186 if (!fPedestals)
187 {
188 *fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl;
189 return kFALSE;
190 }
191
192 return kTRUE;
193}
194
195
196// --------------------------------------------------------------------------
197//
198// The ReInit searches for the following input containers:
199// - MRawRunHeader
200//
201Bool_t MCalibrationCalc::ReInit(MParList *pList )
202{
203
204 fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
205 if (!fRunHeader)
206 {
207 *fLog << dbginf << "MRawRunHeader not found... aborting." << endl;
208 return kFALSE;
209 }
210
211 fNumHiGainSamples = fRunHeader->GetNumSamplesHiGain();
212 fNumLoGainSamples = fRunHeader->GetNumSamplesLoGain();
213
214 //
215 // FIXME: The next statement simply does not work:
216 // fRawEvt->GetNumPixels() returns always 0
217 // fCalibrations->InitSize(fRawEvt->GetNumPixels());
218 //
219
220 fCalibrations->InitSize(577);
221
222 for (Int_t i=0;i<577;i++)
223 {
224 MCalibrationPix &pix = (*fCalibrations)[i];
225 pix.DefinePixId(i);
226 }
227
228 return kTRUE;
229
230}
231
232
233// --------------------------------------------------------------------------
234//
235// Calculate the integral of the FADC time slices and store them as a new
236// pixel in the MCerPhotEvt container.
237//
238Int_t MCalibrationCalc::Process()
239{
240
241 fEvents++;
242
243 Int_t cosmicpix = 0;
244
245 MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
246 MCalibrationPINDiode &pindiode = *(fCalibrations->GetPINDiode());
247
248 MRawEvtPixelIter pixel(fRawEvt);
249
250 //
251 // Create a first loop to sort out the cosmics ...
252 //
253 // This is a very primitive check for the number of cosmicpixs
254 // The cut will be applied in the fit, but for the blind pixel,
255 // we need to remove this event
256 //
257 // FIXME: In the future need a much more sophisticated one!!!
258 //
259
260 while (pixel.Next())
261 {
262
263 const Int_t pixid = pixel.GetPixelId();
264
265 Int_t sum = pixel.GetSumHiGainSamples();
266
267 MPedestalPix &ped = (*fPedestals)[pixid];
268
269 Float_t pedes = ped.GetPedestal();
270 Float_t pedrms = ped.GetPedestalRms();
271
272 if ((float)sum < ((pedes*fNumHiGainSamples)+(30.*pedrms)) )
273 cosmicpix++;
274 }
275
276
277 if (cosmicpix > 50.)
278 {
279 fCosmics++;
280 return kCONTINUE;
281 }
282
283 pixel.Reset();
284
285 //
286 // Create a second loop to do fill the calibration histograms
287 //
288
289 while (pixel.Next())
290 {
291
292 UShort_t sat = 0;
293
294 const Int_t pixid = pixel.GetPixelId();
295
296 Byte_t mid = pixel.GetIdxMaxHiGainSample();
297 UInt_t max = pixel.GetMaxHiGainSample();
298
299 MPedestalPix &ped = (*fPedestals)[pixid];
300 MCalibrationPix &pix = (*fCalibrations)[pixid];
301
302 Float_t pedes = ped.GetPedestal();
303 Float_t sum;
304
305 //
306 // FIXME: This is preliminary, we will change to pedestals per slice!!!
307 // Assume pedestals per time slice ==> multiply with number of slices
308 //
309
310 if (max < gkSaturationLimit) // take Hi Gain
311 {
312 sum = (float)pixel.GetSumHiGainSamples() - pedes*fNumHiGainSamples;
313 }
314 else // Lo Gain
315 {
316
317 sat++;
318 pix.SetHiGainSaturation();
319
320 sum = (float)pixel.GetSumLoGainSamples() - pedes*fNumLoGainSamples ;
321 sum *= gkConversionHiLo;
322
323 max = pixel.GetMaxLoGainSample();
324 mid = pixel.GetIdxMaxLoGainSample();
325
326 if (max > gkSaturationLimit)
327 {
328 *fLog << err << dbginf
329 << "Warning: Saturation of Lo Gain reached in pixel: "
330 << pixid << " " << " sum = " << sum << endl;
331 fHistOverFlow++;
332 }
333
334 }
335
336 switch(pixid)
337 {
338
339 case gkCalibrationBlindPixelId:
340
341 if (!blindpixel.FillCharge(sum))
342 *fLog << warn <<
343 "Overflow or Underflow occurred filling Blind Pixel sum = " << sum << endl;
344
345 if (!blindpixel.FillTime((int)mid))
346 *fLog << warn <<
347 "Overflow or Underflow occurred filling Blind Pixel time = " << (int)mid << endl;
348
349 if (!blindpixel.FillRChargevsTime(sum,fEvents))
350 *fLog << warn <<
351 "Overflow or Underflow occurred filling Blind Pixel eventnr = " << fEvents << endl;
352
353 case gkCalibrationPINDiodeId:
354 if (!pindiode.FillCharge(sum))
355 *fLog << warn <<
356 "Overflow or Underflow occurred filling HCharge: means = " << sum << endl;
357 if (!pindiode.FillTime((int)mid))
358 *fLog << warn <<
359 "Overflow or Underflow occurred filling HTime: time = " << (int)mid << endl;
360 if (!pindiode.FillRChargevsTime(sum,fEvents))
361 *fLog << warn <<
362 "Overflow or Underflow occurred filling HChargevsN: eventnr = " << fEvents << endl;
363
364 default:
365
366 if (sat)
367 {
368
369 if (!pix.FillChargeLoGain(sum))
370 *fLog << warn << "Could not fill Lo Gain Charge of pixel: " << pixid
371 << " signal = " << sum << endl;
372
373 if (!pix.FillTimeLoGain((int)mid))
374 *fLog << warn << "Could not fill Lo Gain Time of pixel: "
375 << pixid << " time = " << (int)mid << endl;
376
377 //
378 // Fill the reduced charge into the control histo for better visibility
379 //
380 if (!pix.FillRChargevsTimeLoGain(sum,fEvents))
381 *fLog << warn << "Could not fill Lo Gain Charge vs. EvtNr of pixel: "
382 << pixid << " signal = " << sum << " event Nr: " << fEvents << endl;
383
384 }
385 else
386 {
387 if (!pix.FillChargeHiGain(sum))
388 *fLog << warn << "Could not fill Hi Gain Charge of pixel: " << pixid
389 << " signal = " << sum << endl;
390
391 if (!pix.FillTimeHiGain((int)mid))
392 *fLog << warn << "Could not fill Hi Gain Time of pixel: "
393 << pixid << " time = " << (int)mid << endl;
394
395 if (!pix.FillRChargevsTimeHiGain(sum,fEvents))
396 *fLog << warn << "Could not fill Hi Gain Charge vs. EvtNr of pixel: "
397 << pixid << " signal = " << sum << " event Nr: " << fEvents << endl;
398 }
399
400 } /* switch(pixid) */
401
402 } /* while (pixel.Next()) */
403
404 return kTRUE;
405}
406
407Int_t MCalibrationCalc::PostProcess()
408{
409
410 *fLog << inf << endl;
411 *fLog << GetDescriptor() << " Cut Histogram Edges" << endl;
412
413 //
414 // Cut edges to make fits and viewing of the hists easier
415 //
416 fCalibrations->CutEdges();
417
418 //
419 // Get pointer to blind pixel
420 //
421 MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
422
423 *fLog << GetDescriptor() << " Fitting the Blind Pixel" << endl;
424
425 //
426 // Fit the blind pixel
427 //
428 if (TESTBIT(fFlags,kUseBlindPixelFit))
429 {
430 if (blindpixel.FitCharge())
431 {
432 if (!fCalibrations->CalcNrPhotInnerPixel())
433 *fLog << err << dbginf << "Could not calculate Number of photons from the blind pixel " << endl;
434 }
435 else
436 {
437 *fLog << err << dbginf << "Could not fit the blind pixel " << endl;
438 }
439
440 if (!blindpixel.FitTime())
441 *fLog << warn << "Could not the Times of the blind pixel " << endl;
442
443 blindpixel.Draw();
444
445 }
446
447 *fLog << GetDescriptor() << " Fitting the Normal Pixels" << endl;
448
449 //
450 // loop over the pedestal events and check if we have calibration
451 //
452 for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
453 {
454
455 MCalibrationPix &pix = (*fCalibrations)[pixid];
456
457 const Float_t ped = (*fPedestals)[pixid].GetPedestal() * fNumHiGainSamples;
458 const Float_t prms = (*fPedestals)[pixid].GetPedestalRms() * TMath::Sqrt((float)fNumHiGainSamples);
459
460 pix.SetPedestal(ped,prms);
461
462 if (TESTBIT(fFlags,kUseTimeFits))
463 pix.FitTime();
464
465 pix.FitCharge();
466 }
467
468 fCalibrations->SetReadyToSave();
469
470 if (GetNumExecutions()==0)
471 return kTRUE;
472
473 *fLog << endl;
474 *fLog << dec << setfill(' ') << fCosmics << " Events presumably cosmics" << endl;
475
476 return kTRUE;
477}
Note: See TracBrowser for help on using the repository browser.