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

Last change on this file since 2844 was 2840, checked in by gaug, 21 years ago
*** empty log message ***
File size: 16.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// 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#include "MCalibrationCalc.h"
73#include "MCalibrationConfig.h"
74
75#include "MCalibrationCam.h"
76#include "MCalibrationPix.h"
77#include "MCalibrationBlindPix.h"
78#include "MCalibrationPINDiode.h"
79
80#include "MPedestalCam.h"
81#include "MPedestalPix.h"
82
83#include "MGeomCam.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"
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.
107//
108MCalibrationCalc::MCalibrationCalc(const char *name, const char *title)
109 : fPedestals(NULL), fCalibrations(NULL), fSignals(NULL),
110 fRawEvt(NULL), fRunHeader(NULL), fEvtTime(NULL),
111 fEvents(0), fHistOverFlow(0), fCosmics(0),
112 fNumHiGainSamples(0), fNumLoGainSamples(0), fConversionHiLo(0.),
113 fColor(kEBlue)
114{
115
116 fName = name ? name : "MCalibrationCalc";
117 fTitle = title ? title : "Task to calculate the calibration constants and MCalibrationCam ";
118
119 AddToBranchList("MRawEvtData.fHiGainPixId");
120 AddToBranchList("MRawEvtData.fLoGainPixId");
121 AddToBranchList("MRawEvtData.fHiGainFadcSamples");
122 AddToBranchList("MRawEvtData.fLoGainFadcSamples");
123
124 SETBIT(fFlags, kUseTimeFits);
125 SETBIT(fFlags, kUseBlindPixelFit);
126 SETBIT(fFlags, kUsePinDiodeFit);
127
128}
129
130// --------------------------------------------------------------------------
131//
132// The PreProcess searches for the following input containers:
133// - MRawEvtData
134// - MPedestalCam
135//
136// The following output containers are also searched and created if
137// they were not found:
138//
139// - MHCalibrationBlindPixel
140// - MCalibrationCam
141// - MTime
142//
143Int_t MCalibrationCalc::PreProcess(MParList *pList)
144{
145
146 fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
147 if (!fRawEvt)
148 {
149 *fLog << err << 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 << err << dbginf << ": MRawRunHeader not found... aborting." << endl;
219 return kFALSE;
220 }
221
222
223 MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
224 if (!cam)
225 {
226 *fLog << err << GetDescriptor() << ": No MGeomCam found... aborting." << endl;
227 return kFALSE;
228 }
229
230
231 fNumHiGainSamples = fSignals->GetNumUsedHiGainFADCSlices();
232 fNumLoGainSamples = fSignals->GetNumUsedLoGainFADCSlices();
233
234 fSqrtHiGainSamples = TMath::Sqrt((Float_t)fNumHiGainSamples);
235
236 fCalibrations->InitSize(cam->GetNumPixels());
237
238 for (UInt_t i=0;i<cam->GetNumPixels();i++)
239 {
240
241 MCalibrationPix &pix = (*fCalibrations)[i];
242 pix.DefinePixId(i);
243 MHCalibrationPixel *hist = pix.GetHist();
244
245 hist->SetTimeFitRangesHiGain(fSignals->GetFirstUsedSliceHiGain(),
246 fSignals->GetLastUsedSliceHiGain());
247 hist->SetTimeFitRangesLoGain(fSignals->GetFirstUsedSliceLoGain(),
248 fSignals->GetLastUsedSliceLoGain());
249
250 }
251
252 return kTRUE;
253
254}
255
256// --------------------------------------------------------------------------
257//
258// Calculate the integral of the FADC time slices and store them as a new
259// pixel in the MCerPhotEvt container.
260//
261Int_t MCalibrationCalc::Process()
262{
263
264 Int_t cosmicpix = 0;
265
266 MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
267 MCalibrationPINDiode &pindiode = *(fCalibrations->GetPINDiode());
268
269 MRawEvtPixelIter pixel(fRawEvt);
270
271 //
272 // Create a first loop to sort out the cosmics ...
273 //
274 // This is a very primitive check for the number of cosmicpixs
275 // The cut will be applied in the fit, but for the blind pixel,
276 // we need to remove this event
277 //
278 // FIXME: In the future need a much more sophisticated one!!!
279 //
280
281 while (pixel.Next())
282 {
283
284 const UInt_t pixid = pixel.GetPixelId();
285
286 MExtractedSignalPix &sig = (*fSignals)[pixid];
287 MPedestalPix &ped = (*fPedestals)[pixid];
288 Float_t pedrms = ped.GetPedestalRms()*fSqrtHiGainSamples;
289 Float_t sumhi = sig.GetExtractedSignalHiGain();
290
291 //
292 // We consider a pixel as presumably due to cosmics
293 // if its sum of FADC slices is lower than 3 pedestal RMS
294 //
295 if (sumhi < 3.*pedrms )
296 cosmicpix++;
297 }
298
299 //
300 // If the camera contains more than 230
301 // (this is the number of outer pixels plus about 50 inner ones)
302 // presumed pixels due to cosmics, then the event is discarted.
303 // This procedure is more or less equivalent to keeping only events
304 // with at least 350 pixels with high signals.
305 //
306 if (cosmicpix > 230.)
307 {
308 fCosmics++;
309 return kCONTINUE;
310 }
311
312 pixel.Reset();
313 fEvents++;
314
315
316 //
317 // Create a second loop to do fill the calibration histograms
318 //
319
320 while (pixel.Next())
321 {
322
323 const UInt_t pixid = pixel.GetPixelId();
324
325 MExtractedSignalPix &sig = (*fSignals)[pixid];
326
327 Float_t sumhi = sig.GetExtractedSignalHiGain();
328 Float_t sumlo = sig.GetExtractedSignalLoGain();
329 Float_t mtime = sig.GetMeanArrivalTime();
330
331 MCalibrationPix &pix = (*fCalibrations)[pixid];
332
333 switch(pixid)
334 {
335
336 case gkCalibrationBlindPixelId:
337
338 if (!blindpixel.FillCharge(sumhi))
339 *fLog << warn <<
340 "Overflow or Underflow occurred filling Blind Pixel sum = " << sumhi << endl;
341
342 if (!blindpixel.FillTime((int)mtime))
343 *fLog << warn <<
344 "Overflow or Underflow occurred filling Blind Pixel time = " << mtime << endl;
345
346 if (!blindpixel.FillRChargevsTime(sumhi,fEvents))
347 *fLog << warn <<
348 "Overflow or Underflow occurred filling Blind Pixel eventnr = " << fEvents << endl;
349 break;
350
351 case gkCalibrationPINDiodeId:
352 if (!pindiode.FillCharge(sumhi))
353 *fLog << warn <<
354 "Overflow or Underflow occurred filling PINDiode: sum = " << sumhi << endl;
355 if (!pindiode.FillTime((int)mtime))
356 *fLog << warn <<
357 "Overflow or Underflow occurred filling PINDiode: time = " << mtime << endl;
358 if (!pindiode.FillRChargevsTime(sumhi,fEvents))
359 *fLog << warn <<
360 "Overflow or Underflow occurred filling PINDiode: eventnr = " << fEvents << endl;
361 break;
362
363 default:
364
365 pix.SetChargesInGraph(sumhi,sumlo);
366
367 if (!pix.FillRChargevsTimeLoGain(sumlo,fEvents))
368 *fLog << warn << "Could not fill Lo Gain Charge vs. EvtNr of pixel: "
369 << pixid << " signal = " << sumlo << " event Nr: " << fEvents << endl;
370
371 if (!pix.FillRChargevsTimeHiGain(sumhi,fEvents))
372 *fLog << warn << "Could not fill Hi Gain Charge vs. EvtNr of pixel: "
373 << pixid << " signal = " << sumhi << " event Nr: " << fEvents << endl;
374
375 if (sig.IsLoGainUsed())
376 {
377
378 if (!pix.FillChargeLoGain(sumlo))
379 *fLog << warn << "Could not fill Lo Gain Charge of pixel: " << pixid
380 << " signal = " << sumlo << endl;
381
382 if (!pix.FillTimeLoGain((int)mtime))
383 *fLog << warn << "Could not fill Lo Gain Time of pixel: "
384 << pixid << " time = " << mtime << endl;
385
386 }
387 else
388 {
389 if (!pix.FillChargeHiGain(sumhi))
390 *fLog << warn << "Could not fill Hi Gain Charge of pixel: " << pixid
391 << " signal = " << sumhi << endl;
392
393 if (!pix.FillTimeHiGain((int)mtime))
394 *fLog << warn << "Could not fill Hi Gain Time of pixel: "
395 << pixid << " time = " << mtime << endl;
396
397 }
398 break;
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
412 if (fEvents == 0)
413 {
414 *fLog << err << GetDescriptor()
415 << ": This run contains only cosmics or pedestals, "
416 << "cannot find events with more than 350 illuminated pixels. " << endl;
417 return kFALSE;
418 }
419
420 if (fEvents < fCosmics)
421 *fLog << warn << GetDescriptor()
422 << ": WARNING: Run contains more cosmics or pedestals than calibration events " << endl;
423
424
425 *fLog << inf << GetDescriptor() << ": Cut Histogram Edges" << endl;
426
427 //
428 // Cut edges to make fits and viewing of the hists easier
429 //
430 fCalibrations->CutEdges();
431
432 //
433 // Get pointer to blind pixel
434 //
435 MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
436
437 *fLog << inf << GetDescriptor() << ": Fitting the Blind Pixel" << endl;
438
439 //
440 // Fit the blind pixel
441 //
442 if (TESTBIT(fFlags,kUseBlindPixelFit))
443 {
444
445 if (!blindpixel.FitCharge())
446 {
447 *fLog << err << dbginf << "Could not fit the blind pixel! " << endl;
448 blindpixel.DrawClone();
449 return kFALSE;
450 }
451
452 blindpixel.DrawClone();
453 }
454
455 *fLog << inf << GetDescriptor() << ": Fitting the Normal Pixels" << endl;
456
457 //
458 // loop over the pedestal events and check if we have calibration
459 //
460 for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
461 {
462
463 MCalibrationPix &pix = (*fCalibrations)[pixid];
464
465 //
466 // get the pedestals
467 //
468 const Float_t ped = (*fPedestals)[pixid].GetPedestal() * fNumHiGainSamples;
469 const Float_t prms = (*fPedestals)[pixid].GetPedestalRms() * fSqrtHiGainSamples;
470
471 //
472 // set them in the calibration camera
473 //
474 pix.SetPedestal(ped,prms);
475
476
477 //
478 // Check if the pixel has been excluded from the fits
479 //
480 if (pix.IsExcluded())
481 continue;
482
483 //
484 // perform the Gauss fits to the charges
485 //
486 pix.FitCharge();
487
488 //
489 // Perform the Gauss fits to the arrival times
490 //
491 if (TESTBIT(fFlags,kUseTimeFits))
492 pix.FitTime();
493
494 }
495
496 if (!fCalibrations->CalcNumPhotInsidePlexiglass())
497 {
498 *fLog << err << GetDescriptor()
499 << ": Could not calculate the number of photons from the blind pixel " << endl;
500 return kFALSE;
501 }
502
503 fCalibrations->SetReadyToSave();
504
505 if (GetNumExecutions()==0)
506 return kTRUE;
507
508 *fLog << inf << endl;
509 *fLog << dec << setfill(' ') << fCosmics << " Events presumably cosmics" << endl;
510
511 return kTRUE;
512}
513
Note: See TracBrowser for help on using the repository browser.