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

Last change on this file since 2830 was 2809, checked in by gaug, 22 years ago
*** empty log message ***
File size: 16.0 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"
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.
109//
110MCalibrationCalc::MCalibrationCalc(const char *name, const char *title)
111 : fColor(kEBlue)
112{
113
114 fName = name ? name : "MCalibrationCalc";
115 fTitle = title ? title : "Task to calculate the calibration constants and MCalibrationCam ";
116
117 AddToBranchList("MRawEvtData.fHiGainPixId");
118 AddToBranchList("MRawEvtData.fLoGainPixId");
119 AddToBranchList("MRawEvtData.fHiGainFadcSamples");
120 AddToBranchList("MRawEvtData.fLoGainFadcSamples");
121
122 SETBIT(fFlags, kUseTimeFits);
123 SETBIT(fFlags, kUseBlindPixelFit);
124 SETBIT(fFlags, kUsePinDiodeFit);
125
126}
127
128// --------------------------------------------------------------------------
129//
130// The PreProcess searches for the following input containers:
131// - MRawEvtData
132// - MPedestalCam
133//
134// The following output containers are also searched and created if
135// they were not found:
136//
137// - MHCalibrationBlindPixel
138// - MCalibrationCam
139// - MTime
140//
141Int_t MCalibrationCalc::PreProcess(MParList *pList)
142{
143
144 fHistOverFlow = 0;
145 fEvents = 0;
146 fCosmics = 0;
147
148 fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
149 if (!fRawEvt)
150 {
151 *fLog << dbginf << "MRawEvtData not found... aborting." << endl;
152 return kFALSE;
153 }
154
155 const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
156 if (!runheader)
157 *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl;
158 else
159 if (runheader->GetRunType() == kRTMonteCarlo)
160 {
161 return kTRUE;
162 }
163
164 fCalibrations = (MCalibrationCam*)pList->FindCreateObj("MCalibrationCam");
165 if (!fCalibrations)
166 {
167 *fLog << err << dbginf << "MCalibrationCam could not be created ... aborting." << endl;
168 return kFALSE;
169 }
170
171
172 switch (fColor)
173 {
174 case kEBlue:
175 fCalibrations->SetColor(MCalibrationCam::kECBlue);
176 break;
177 case kEGreen:
178 fCalibrations->SetColor(MCalibrationCam::kECGreen);
179 break;
180 case kEUV:
181 fCalibrations->SetColor(MCalibrationCam::kECUV);
182 break;
183 case kECT1:
184 fCalibrations->SetColor(MCalibrationCam::kECCT1);
185 break;
186 default:
187 fCalibrations->SetColor(MCalibrationCam::kECCT1);
188 }
189
190 fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
191 if (!fPedestals)
192 {
193 *fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl;
194 return kFALSE;
195 }
196
197
198 fSignals = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
199 if (!fSignals)
200 {
201 *fLog << err << dbginf << "Cannot find MExtractedSignalCam ... aborting" << endl;
202 return kFALSE;
203 }
204
205 return kTRUE;
206}
207
208
209// --------------------------------------------------------------------------
210//
211// The ReInit searches for the following input containers:
212// - MRawRunHeader
213//
214Bool_t MCalibrationCalc::ReInit(MParList *pList )
215{
216
217 fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
218 if (!fRunHeader)
219 {
220 *fLog << dbginf << "MRawRunHeader not found... aborting." << endl;
221 return kFALSE;
222 }
223
224
225 MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
226 if (!cam)
227 {
228 *fLog << err << GetDescriptor() << ": No MGeomCam found... aborting." << endl;
229 return kFALSE;
230 }
231
232
233 fNumHiGainSamples = fSignals->GetNumUsedHiGainFADCSlices();
234 fNumLoGainSamples = fSignals->GetNumUsedLoGainFADCSlices();
235
236 fSqrtHiGainSamples = TMath::Sqrt((Float_t)fNumHiGainSamples);
237
238 fCalibrations->InitSize(cam->GetNumPixels());
239
240 for (UInt_t i=0;i<cam->GetNumPixels();i++)
241 {
242
243 MCalibrationPix &pix = (*fCalibrations)[i];
244 pix.DefinePixId(i);
245 MHCalibrationPixel *hist = pix.GetHist();
246
247 hist->SetTimeFitRangesHiGain(fSignals->GetFirstUsedSliceHiGain(),
248 fSignals->GetLastUsedSliceHiGain());
249 hist->SetTimeFitRangesLoGain(fSignals->GetFirstUsedSliceLoGain(),
250 fSignals->GetLastUsedSliceLoGain());
251
252 }
253
254 return kTRUE;
255
256}
257
258
259// --------------------------------------------------------------------------
260//
261// Calculate the integral of the FADC time slices and store them as a new
262// pixel in the MCerPhotEvt container.
263//
264Int_t MCalibrationCalc::Process()
265{
266
267 Int_t cosmicpix = 0;
268
269 MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
270 MCalibrationPINDiode &pindiode = *(fCalibrations->GetPINDiode());
271
272 MRawEvtPixelIter pixel(fRawEvt);
273
274 //
275 // Create a first loop to sort out the cosmics ...
276 //
277 // This is a very primitive check for the number of cosmicpixs
278 // The cut will be applied in the fit, but for the blind pixel,
279 // we need to remove this event
280 //
281 // FIXME: In the future need a much more sophisticated one!!!
282 //
283
284 while (pixel.Next())
285 {
286
287 const UInt_t pixid = pixel.GetPixelId();
288
289 MExtractedSignalPix &sig = (*fSignals)[pixid];
290 MPedestalPix &ped = (*fPedestals)[pixid];
291 Float_t pedrms = ped.GetPedestalRms()*fSqrtHiGainSamples;
292 Float_t sumhi = sig.GetExtractedSignalHiGain();
293
294 //
295 // We consider a pixel as presumably due to cosmics
296 // if its sum of FADC slices is lower than 3 pedestal RMS
297 //
298 if (sumhi < 3.*pedrms )
299 cosmicpix++;
300 }
301
302 //
303 // If the camera contains more than 230
304 // (this is the number of outer pixels plus about 50 inner ones)
305 // presumed pixels due to cosmics, then the event is discarted.
306 // This procedure is more or less equivalent to keeping only events
307 // with at least 350 pixels with high signals.
308 //
309 if (cosmicpix > 230.)
310 {
311 fCosmics++;
312 return kCONTINUE;
313 }
314
315 pixel.Reset();
316 fEvents++;
317
318
319 //
320 // Create a second loop to do fill the calibration histograms
321 //
322
323 while (pixel.Next())
324 {
325
326 const UInt_t pixid = pixel.GetPixelId();
327
328 MExtractedSignalPix &sig = (*fSignals)[pixid];
329
330 Float_t sumhi = sig.GetExtractedSignalHiGain();
331 Float_t sumlo = sig.GetExtractedSignalLoGain();
332 Float_t mtime = sig.GetMeanArrivalTime();
333
334 MCalibrationPix &pix = (*fCalibrations)[pixid];
335
336 switch(pixid)
337 {
338
339 case gkCalibrationBlindPixelId:
340
341 if (!blindpixel.FillCharge(sumhi))
342 *fLog << err <<
343 "Overflow or Underflow occurred filling Blind Pixel sum = " << sumhi << endl;
344
345 if (!blindpixel.FillTime((int)mtime))
346 *fLog << err <<
347 "Overflow or Underflow occurred filling Blind Pixel time = " << mtime << endl;
348
349 if (!blindpixel.FillRChargevsTime(sumhi,fEvents))
350 *fLog << warn <<
351 "Overflow or Underflow occurred filling Blind Pixel eventnr = " << fEvents << endl;
352 break;
353
354 case gkCalibrationPINDiodeId:
355 if (!pindiode.FillCharge(sumhi))
356 *fLog << warn <<
357 "Overflow or Underflow occurred filling PINDiode: sum = " << sumhi << endl;
358 if (!pindiode.FillTime((int)mtime))
359 *fLog << warn <<
360 "Overflow or Underflow occurred filling PINDiode: time = " << mtime << endl;
361 if (!pindiode.FillRChargevsTime(sumhi,fEvents))
362 *fLog << warn <<
363 "Overflow or Underflow occurred filling PINDiode: eventnr = " << fEvents << endl;
364 break;
365
366 default:
367
368 pix.SetChargesInGraph(sumhi,sumlo);
369
370 if (!pix.FillRChargevsTimeLoGain(sumlo,fEvents))
371 *fLog << warn << "Could not fill Lo Gain Charge vs. EvtNr of pixel: "
372 << pixid << " signal = " << sumlo << " event Nr: " << fEvents << endl;
373
374 if (!pix.FillRChargevsTimeHiGain(sumhi,fEvents))
375 *fLog << warn << "Could not fill Hi Gain Charge vs. EvtNr of pixel: "
376 << pixid << " signal = " << sumhi << " event Nr: " << fEvents << endl;
377
378 if (sig.IsLoGainUsed())
379 {
380
381 if (!pix.FillChargeLoGain(sumlo))
382 *fLog << warn << "Could not fill Lo Gain Charge of pixel: " << pixid
383 << " signal = " << sumlo << endl;
384
385 if (!pix.FillTimeLoGain((int)mtime))
386 *fLog << warn << "Could not fill Lo Gain Time of pixel: "
387 << pixid << " time = " << mtime << endl;
388
389 }
390 else
391 {
392 if (!pix.FillChargeHiGain(sumhi))
393 *fLog << warn << "Could not fill Hi Gain Charge of pixel: " << pixid
394 << " signal = " << sumhi << endl;
395
396 if (!pix.FillTimeHiGain((int)mtime))
397 *fLog << warn << "Could not fill Hi Gain Time of pixel: "
398 << pixid << " time = " << mtime << endl;
399
400 }
401 break;
402
403 } /* switch(pixid) */
404
405 } /* while (pixel.Next()) */
406
407 return kTRUE;
408}
409
410Int_t MCalibrationCalc::PostProcess()
411{
412
413
414 *fLog << inf << endl;
415
416 if (fEvents == 0)
417 {
418 *fLog << err << GetDescriptor()
419 << ": This run contains only cosmics or pedestals, "
420 << "cannot find events with more than 350 illuminated pixels. " << endl;
421 return kFALSE;
422 }
423
424 if (fEvents < fCosmics)
425 *fLog << warn << GetDescriptor()
426 << ": WARNING: Run contains more cosmics or pedestals than calibration events " << endl;
427
428
429 *fLog << GetDescriptor() << " Cut Histogram Edges" << endl;
430
431 //
432 // Cut edges to make fits and viewing of the hists easier
433 //
434 fCalibrations->CutEdges();
435
436 //
437 // Get pointer to blind pixel
438 //
439 MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
440
441 *fLog << GetDescriptor() << " Fitting the Blind Pixel" << endl;
442
443 //
444 // Fit the blind pixel
445 //
446 if (TESTBIT(fFlags,kUseBlindPixelFit))
447 {
448
449 if (!blindpixel.FitCharge())
450 *fLog << err << dbginf << "Could not fit the blind pixel " << endl;
451
452 blindpixel.DrawClone();
453 }
454
455 *fLog << 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 // perform the Gauss fits to the charges
478 //
479 pix.FitCharge();
480
481 //
482 // Perform the Gauss fits to the arrival times
483 //
484 if (TESTBIT(fFlags,kUseTimeFits))
485 pix.FitTime();
486
487 }
488
489 if (!fCalibrations->CalcNumPhotInsidePlexiglass())
490 *fLog << err << dbginf << "Could not calculate the number of photons from the blind pixel " << endl;
491
492 fCalibrations->SetReadyToSave();
493
494 if (GetNumExecutions()==0)
495 return kTRUE;
496
497 *fLog << endl;
498 *fLog << dec << setfill(' ') << fCosmics << " Events presumably cosmics" << endl;
499
500 return kTRUE;
501}
Note: See TracBrowser for help on using the repository browser.