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

Last change on this file since 2648 was 2642, checked in by gaug, 21 years ago
*** empty log message ***
File size: 14.9 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 "MExtractedSignalCam.h"
96#include "MExtractedSignalPix.h"
97
98#include "MTime.h"
99#include "TMath.h"
100
101ClassImp(MCalibrationCalc);
102
103using namespace std;
104// --------------------------------------------------------------------------
105//
106// Default constructor. b is the number of slices before the maximum slice,
107// a the number of slices behind the maximum slice which is taken as signal.
108//
109MCalibrationCalc::MCalibrationCalc(const char *name, const char *title)
110 : fColor(kEBlue)
111{
112
113 fName = name ? name : "MCalibrationCalc";
114 fTitle = title ? title : "Task to calculate the calibration constants and MCalibrationCam ";
115
116 AddToBranchList("MRawEvtData.fHiGainPixId");
117 AddToBranchList("MRawEvtData.fLoGainPixId");
118 AddToBranchList("MRawEvtData.fHiGainFadcSamples");
119 AddToBranchList("MRawEvtData.fLoGainFadcSamples");
120
121 SETBIT(fFlags, kUseTimeFits);
122 SETBIT(fFlags, kUseBlindPixelFit);
123 SETBIT(fFlags, kUsePinDiodeFit);
124}
125
126// --------------------------------------------------------------------------
127//
128// The PreProcess searches for the following input containers:
129// - MRawEvtData
130// - MPedestalCam
131//
132// The following output containers are also searched and created if
133// they were not found:
134//
135// - MHCalibrationBlindPixel
136// - MCalibrationCam
137// - MTime
138//
139Int_t MCalibrationCalc::PreProcess(MParList *pList)
140{
141
142 fHistOverFlow = 0;
143 fEvents = 0;
144 fCosmics = 0;
145
146 fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
147 if (!fRawEvt)
148 {
149 *fLog << dbginf << "MRawEvtData not found... aborting." << endl;
150 return kFALSE;
151 }
152
153 const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
154 if (!runheader)
155 *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl;
156 else
157 if (runheader->GetRunType() == kRTMonteCarlo)
158 {
159 return kTRUE;
160 }
161
162 fCalibrations = (MCalibrationCam*)pList->FindCreateObj("MCalibrationCam");
163 if (!fCalibrations)
164 {
165 *fLog << err << dbginf << "MCalibrationCam could not be created ... aborting." << endl;
166 return kFALSE;
167 }
168
169
170 switch (fColor)
171 {
172 case kEBlue:
173 fCalibrations->SetColor(MCalibrationCam::kECBlue);
174 break;
175 case kEGreen:
176 fCalibrations->SetColor(MCalibrationCam::kECGreen);
177 break;
178 case kEUV:
179 fCalibrations->SetColor(MCalibrationCam::kECUV);
180 break;
181 case kECT1:
182 fCalibrations->SetColor(MCalibrationCam::kECCT1);
183 break;
184 default:
185 fCalibrations->SetColor(MCalibrationCam::kECCT1);
186 }
187
188 fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
189 if (!fPedestals)
190 {
191 *fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl;
192 return kFALSE;
193 }
194
195
196 fSignals = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
197 if (!fSignals)
198 {
199 *fLog << err << dbginf << "Cannot find MExtractedSignalCam ... aborting" << endl;
200 return kFALSE;
201 }
202
203 return kTRUE;
204}
205
206
207// --------------------------------------------------------------------------
208//
209// The ReInit searches for the following input containers:
210// - MRawRunHeader
211//
212Bool_t MCalibrationCalc::ReInit(MParList *pList )
213{
214
215 fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
216 if (!fRunHeader)
217 {
218 *fLog << dbginf << "MRawRunHeader not found... aborting." << endl;
219 return kFALSE;
220 }
221
222 // fNumHiGainSamples = fRunHeader->GetNumSamplesHiGain();
223 // fNumLoGainSamples = fRunHeader->GetNumSamplesLoGain();
224 fNumHiGainSamples = 6;
225 fNumLoGainSamples = 6;
226
227 //
228 // FIXME: The next statement simply does not work:
229 // fRawEvt->GetNumPixels() returns always 0
230 // fCalibrations->InitSize(fRawEvt->GetNumPixels());
231 //
232
233 fCalibrations->InitSize(577);
234
235 for (Int_t i=0;i<577;i++)
236 {
237 MCalibrationPix &pix = (*fCalibrations)[i];
238 pix.DefinePixId(i);
239 }
240
241 return kTRUE;
242
243}
244
245
246// --------------------------------------------------------------------------
247//
248// Calculate the integral of the FADC time slices and store them as a new
249// pixel in the MCerPhotEvt container.
250//
251Int_t MCalibrationCalc::Process()
252{
253
254 fEvents++;
255
256 Int_t cosmicpix = 0;
257
258 MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
259 MCalibrationPINDiode &pindiode = *(fCalibrations->GetPINDiode());
260
261 MRawEvtPixelIter pixel(fRawEvt);
262
263 //
264 // Create a first loop to sort out the cosmics ...
265 //
266 // This is a very primitive check for the number of cosmicpixs
267 // The cut will be applied in the fit, but for the blind pixel,
268 // we need to remove this event
269 //
270 // FIXME: In the future need a much more sophisticated one!!!
271 //
272
273 while (pixel.Next())
274 {
275
276 const Int_t pixid = pixel.GetPixelId();
277
278 Int_t sum = pixel.GetSumHiGainSamples();
279
280 MPedestalPix &ped = (*fPedestals)[pixid];
281
282 Float_t pedes = ped.GetPedestal();
283 Float_t pedrms = ped.GetPedestalRms();
284
285 if ((float)sum < ((pedes*fNumHiGainSamples)+(2.*fNumHiGainSamples*pedrms)) )
286 cosmicpix++;
287 }
288
289
290 if (cosmicpix > 50.)
291 {
292 fCosmics++;
293 return kCONTINUE;
294 }
295
296 pixel.Reset();
297
298 //
299 // Create a second loop to do fill the calibration histograms
300 //
301
302 while (pixel.Next())
303 {
304
305 const UInt_t pixid = pixel.GetPixelId();
306
307 MExtractedSignalPix &sig = (*fSignals)[pixid];
308
309 Float_t sumhi = sig.GetExtractedSignalHiGain();
310 Float_t sumlo = sig.GetExtractedSignalLoGain();
311 Bool_t logain = sig.IsLoGainUsed();
312
313 Byte_t mid;
314
315 if (logain)
316 mid = pixel.GetIdxMaxLoGainSample();
317 else
318 mid = pixel.GetIdxMaxHiGainSample();
319
320 MCalibrationPix &pix = (*fCalibrations)[pixid];
321
322 switch(pixid)
323 {
324
325 case gkCalibrationBlindPixelId:
326
327 if (!blindpixel.FillCharge(sumhi))
328 *fLog << warn <<
329 "Overflow or Underflow occurred filling Blind Pixel sum = " << sumhi << endl;
330
331 if (!blindpixel.FillTime((int)mid))
332 *fLog << warn <<
333 "Overflow or Underflow occurred filling Blind Pixel time = " << (int)mid << endl;
334
335 if (!blindpixel.FillRChargevsTime(sumhi,fEvents))
336 *fLog << warn <<
337 "Overflow or Underflow occurred filling Blind Pixel eventnr = " << fEvents << endl;
338
339 case gkCalibrationPINDiodeId:
340 if (!pindiode.FillCharge(sumhi))
341 *fLog << warn <<
342 "Overflow or Underflow occurred filling HCharge: means = " << sumhi << endl;
343 if (!pindiode.FillTime((int)mid))
344 *fLog << warn <<
345 "Overflow or Underflow occurred filling HTime: time = " << (int)mid << endl;
346 if (!pindiode.FillRChargevsTime(sumhi,fEvents))
347 *fLog << warn <<
348 "Overflow or Underflow occurred filling HChargevsN: eventnr = " << fEvents << endl;
349
350 default:
351
352 pix.SetChargesInGraph(sumhi,sumlo);
353
354 if (logain)
355 {
356
357 if (!pix.FillChargeLoGain(sumlo))
358 *fLog << warn << "Could not fill Lo Gain Charge of pixel: " << pixid
359 << " signal = " << sumlo << endl;
360
361 if (!pix.FillTimeLoGain((int)mid))
362 *fLog << warn << "Could not fill Lo Gain Time of pixel: "
363 << pixid << " time = " << (int)mid << endl;
364
365 //
366 // Fill the reduced charge into the control histo for better visibility
367 //
368 if (!pix.FillRChargevsTimeLoGain(sumlo,fEvents))
369 *fLog << warn << "Could not fill Lo Gain Charge vs. EvtNr of pixel: "
370 << pixid << " signal = " << sumlo << " event Nr: " << fEvents << endl;
371
372 }
373 else
374 {
375 if (!pix.FillChargeHiGain(sumhi))
376 *fLog << warn << "Could not fill Hi Gain Charge of pixel: " << pixid
377 << " signal = " << sumhi << endl;
378
379 if (!pix.FillTimeHiGain((int)mid))
380 *fLog << warn << "Could not fill Hi Gain Time of pixel: "
381 << pixid << " time = " << (int)mid << endl;
382
383 if (!pix.FillRChargevsTimeHiGain(sumhi,fEvents))
384 *fLog << warn << "Could not fill Hi Gain Charge vs. EvtNr of pixel: "
385 << pixid << " signal = " << sumhi << " event Nr: " << fEvents << endl;
386 }
387
388 } /* switch(pixid) */
389
390 } /* while (pixel.Next()) */
391
392 return kTRUE;
393}
394
395Int_t MCalibrationCalc::PostProcess()
396{
397
398
399 *fLog << inf << endl;
400 *fLog << GetDescriptor() << " Cut Histogram Edges" << endl;
401
402 //
403 // Cut edges to make fits and viewing of the hists easier
404 //
405 fCalibrations->CutEdges();
406
407 //
408 // Get pointer to blind pixel
409 //
410 MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
411
412 *fLog << GetDescriptor() << " Fitting the Blind Pixel" << endl;
413
414 //
415 // Fit the blind pixel
416 //
417 if (TESTBIT(fFlags,kUseBlindPixelFit))
418 {
419 if (blindpixel.FitCharge())
420 {
421 if (!fCalibrations->CalcNrPhotInnerPixel())
422 *fLog << err << dbginf << "Could not calculate Number of photons from the blind pixel " << endl;
423 }
424 else
425 {
426 *fLog << err << dbginf << "Could not fit the blind pixel " << endl;
427 }
428
429 if (!blindpixel.FitTime())
430 *fLog << warn << "Could not the Times of the blind pixel " << endl;
431
432 blindpixel.Draw();
433
434 }
435
436 *fLog << GetDescriptor() << " Fitting the Normal Pixels" << endl;
437
438 //
439 // loop over the pedestal events and check if we have calibration
440 //
441 for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
442 {
443
444 MCalibrationPix &pix = (*fCalibrations)[pixid];
445
446 const Float_t ped = (*fPedestals)[pixid].GetPedestal() * fNumHiGainSamples;
447 const Float_t prms = (*fPedestals)[pixid].GetPedestalRms() * TMath::Sqrt((float)fNumHiGainSamples);
448
449 pix.SetPedestal(ped,prms);
450
451 if (TESTBIT(fFlags,kUseTimeFits))
452 pix.FitTime();
453
454 pix.FitCharge();
455 }
456
457 fCalibrations->SetReadyToSave();
458
459 if (GetNumExecutions()==0)
460 return kTRUE;
461
462 *fLog << endl;
463 *fLog << dec << setfill(' ') << fCosmics << " Events presumably cosmics" << endl;
464
465 return kTRUE;
466}
Note: See TracBrowser for help on using the repository browser.