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

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