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

Last change on this file since 2590 was 2581, checked in by gaug, 22 years ago
*** empty log message ***
File size: 15.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// HQ(npix) (distribution of summed FADC time slice entries) //
56// HT(npix) (distribution of position of maximum) //
57// HQvsN(npix) (distribution of charges vs. event number. //
58// //
59// PostProcess: All histograms HQ(npix) are fitted to a Gaussian //
60// All histograms HT(npix) are fitted to a Gaussian //
61// The histogram HBPQ (blind pixel) is fitted to a single //
62// PhE fit //
63// The histogram HBPT (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::SetSkipTFits() (skip all time fits) //
68// MalibrationCam::SetSkipBPFits() (skip all blind pixel fits) //
69// MalibrationCam::SetSkipPDFits() (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
97ClassImp(MCalibrationCalc);
98
99using namespace std;
100// --------------------------------------------------------------------------
101//
102// Default constructor. b is the number of slices before the maximum slice,
103// a the number of slices behind the maximum slice which is taken as signal.
104//
105MCalibrationCalc::MCalibrationCalc(const char *name, const char *title)
106 : fColor(kEBlue)
107{
108
109 fName = name ? name : "MCalibrationCalc";
110 fTitle = title ? title : "Task to calculate the calibration constants and MCalibrationCam ";
111
112 AddToBranchList("MRawEvtData.fHiGainPixId");
113 AddToBranchList("MRawEvtData.fLoGainPixId");
114 AddToBranchList("MRawEvtData.fHiGainFadcSamples");
115 AddToBranchList("MRawEvtData.fLoGainFadcSamples");
116
117 SETBIT(fFlags, kUseTFits);
118 SETBIT(fFlags, kUseBPFit);
119 SETBIT(fFlags, kUsePDFit);
120}
121
122// --------------------------------------------------------------------------
123//
124// The PreProcess searches for the following input containers:
125// - MRawEvtData
126// - MPedestalCam
127//
128// The following output containers are also searched and created if
129// they were not found:
130//
131// - MHCalibrationBlindPixel
132// - MCalibrationCam
133// - MTime
134//
135Int_t MCalibrationCalc::PreProcess(MParList *pList)
136{
137
138 fHistOverFlow = 0;
139 fEvents = 0;
140 fCosmics = 0;
141
142 fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
143 if (!fRawEvt)
144 {
145 *fLog << dbginf << "MRawEvtData not found... aborting." << endl;
146 return kFALSE;
147 }
148
149 const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
150 if (!runheader)
151 *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl;
152 else
153 if (runheader->GetRunType() == kRTMonteCarlo)
154 {
155 return kTRUE;
156 }
157
158 fCalibrations = (MCalibrationCam*)pList->FindCreateObj("MCalibrationCam");
159 if (!fCalibrations)
160 {
161 *fLog << err << dbginf << "MCalibrationCam could not be created ... aborting." << endl;
162 return kFALSE;
163 }
164
165 //
166 // fRawEvt->GetNumPixels() does not work !!!!!!!!!!
167 //
168 fCalibrations->InitSize(577);
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 }
181
182
183
184 fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
185 if (!fPedestals)
186 {
187 *fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl;
188 return kFALSE;
189 }
190
191
192 // MTime does not work!!!!
193
194 // fEvtTime = (MTime*)pList->FindObject("MRawEvtTime");
195 // if (!fEvtTime)
196 // {
197 // *fLog << dbginf << "MTime could not be created ... ¡!¡!¡!¡!" << endl;
198 // }
199
200 // fEvtTime->SetTime(0.);
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
225 return kTRUE;
226
227}
228
229
230// --------------------------------------------------------------------------
231//
232// Calculate the integral of the FADC time slices and store them as a new
233// pixel in the MCerPhotEvt container.
234//
235Int_t MCalibrationCalc::Process()
236{
237
238 fEvents++;
239
240 Bool_t cosmic = kFALSE;
241
242 MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
243 MCalibrationPINDiode &pindiode = *(fCalibrations->GetPINDiode());
244
245 MRawEvtPixelIter pixel(fRawEvt);
246
247 while (pixel.Next())
248 {
249
250 UShort_t sat = 0;
251 const Int_t pixid = pixel.GetPixelId();
252
253 if (!fCalibrations->IsPixelUsed(pixid))
254 if (!fCalibrations->AddPixel(pixid))
255 *fLog << err << dbginf << "MCalibrationPix(" << pixid << ") could not be created !!" << endl;
256
257 Byte_t *ptr = pixel.GetHiGainSamples();
258 Byte_t mid = pixel.GetIdxMaxHiGainSample();
259 UInt_t max = pixel.GetMaxHiGainSample();
260
261 Int_t sum = (max > gkSaturationLimit // overflow ?
262 ? sat++, pixel.GetSumLoGainSamples() // take Low Gain
263 : pixel.GetSumHiGainSamples()); // no overflow
264
265 if (sat)
266 {
267
268 ptr = pixel.GetLoGainSamples();
269 max = pixel.GetMaxLoGainSample();
270 mid = pixel.GetIdxMaxLoGainSample();
271
272 sum = (max == gkSaturationLimit // overflow of LoGain ??? -> GimmeABreak!!!
273 ? fHistOverFlow++, gkLoGainOverFlow // OUCH (Florian was maybe right)
274 : sum*gkConversionHiLo ); // OUFF (Florian was wrong) !!
275
276 // *fLog << warn << "Warning: Saturation of HiGain reached in slice " << (int)mid << " !!! " << endl;
277 // *fLog << warn << "Max = " << max << endl;
278
279 if (fHistOverFlow)
280 *fLog << err << dbginf << "Warning: Saturation of LoGain reached! "
281 << err << dbginf << "sum = " << sum << endl;
282
283 }
284
285
286 MPedestalPix &ped = (*fPedestals)[pixid];
287 MCalibrationPix &pix = (*fCalibrations)[pixid];
288
289 Float_t pedes = ped.GetPedestal();
290 Float_t pedrms = ped.GetPedestalRms();
291
292 //
293 // FIXME: This is preliminary, we will change to pedestals per slice!!!
294 // Assume pedestals per time slice ==> multiply with number of slices
295 //
296 pedes *= (sat ? fNumLoGainSamples : fNumHiGainSamples );
297 pedrms *= (sat ? fNumLoGainSamples : fNumHiGainSamples );
298
299 //
300 // This is a very primitive check for the number of cosmics
301 // The cut will be applied in the fit, but for the blind pixel,
302 // we need to remove this event
303 //
304 // FIXME: In the future need a much more sophisticated one!!!
305 //
306
307 if ((float)sum < pedes+5.*pedrms)
308 cosmic = kTRUE;
309
310 Float_t rsum = (float)sum - pedes;
311
312 switch(pixid)
313 {
314
315 case gkCalibrationBlindPixelId:
316 //
317 // FIXME: This works only when the blind pixel ID is much larger than
318 // the rest of the pixels (which is the case right now)
319 //
320// if (!cosmic)
321 if (!blindpixel.FillQ(sum))
322 *fLog << warn <<
323 "Overflow or Underflow occurred filling Blind Pixel sum = " << sum << endl;
324
325 if (!blindpixel.FillT((int)mid))
326 *fLog << warn <<
327 "Overflow or Underflow occurred filling Blind Pixel time = " << (int)mid << endl;
328
329 if (!blindpixel.FillRQvsT(rsum,fEvents))
330 *fLog << warn <<
331 "Overflow or Underflow occurred filling Blind Pixel eventnr = " << fEvents << endl;
332
333 case gkCalibrationPINDiodeId:
334 if (!pindiode.FillQ(sum))
335 *fLog << warn <<
336 "Overflow or Underflow occurred filling HQ: means = " << sum << endl;
337 if (!pindiode.FillT((int)mid))
338 *fLog << warn <<
339 "Overflow or Underflow occurred filling HT: time = " << (int)mid << endl;
340 if (!pindiode.FillRQvsT(rsum,fEvents))
341 *fLog << warn <<
342 "Overflow or Underflow occurred filling HQvsN: eventnr = " << fEvents << endl;
343
344 default:
345
346 if (!pix.FillQ(sum))
347 *fLog << warn << "Could not fill Q of pixel: " << pixid
348 << " signal = " << sum << endl;
349
350 //
351 // Fill the reduced charge into the control histo for better visibility
352 //
353 if (!pix.FillRQvsT(rsum,fEvents))
354 *fLog << warn << "Could not fill red. Q vs. EvtNr of pixel: " << pixid
355 << " signal = " << rsum << " event Nr: " << fEvents << endl;
356
357
358 if (!pix.FillT((int)mid))
359 *fLog << warn << "Could not fill T of pixel: " << pixid << " time = " << (int)mid << endl;
360
361
362 } /* switch(pixid) */
363
364 } /* while (pixel.Next()) */
365
366 if (cosmic)
367 fCosmics++;
368
369 fCalibrations->SetReadyToSave();
370
371 return kTRUE;
372}
373
374Int_t MCalibrationCalc::PostProcess()
375{
376
377 *fLog << inf << endl;
378 *fLog << GetDescriptor() << " Cut Histogram Edges" << endl;
379 //
380 // Cut edges to make fits and viewing of the hists easier
381 //
382 fCalibrations->CutEdges();
383
384 //
385 // Get pointer to blind pixel
386 //
387 MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
388
389 *fLog << GetDescriptor() << " Fitting the Blind Pixel" << endl;
390
391 //
392 // Fit the blind pixel
393 //
394 if (TESTBIT(fFlags,kUseBPFit))
395 {
396 if (!blindpixel.FitQ())
397 *fLog << err << dbginf << "Could not fit the blind pixel " << endl;
398
399 if (!blindpixel.FitT())
400 *fLog << warn << "Could not the Times of the blind pixel " << endl;
401
402 blindpixel.Draw();
403
404 }
405
406 *fLog << GetDescriptor() << " Fitting the Normal Pixels" << endl;
407 //
408 // loop over the pedestal events and check if we have calibration
409 //
410 for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
411 {
412
413 if (fCalibrations->IsPixelUsed(pixid))
414 {
415
416 MCalibrationPix &pix = (*fCalibrations)[pixid];
417
418 const Float_t ped = (*fPedestals)[pixid].GetPedestal() * fNumHiGainSamples;
419 const Float_t prms = (*fPedestals)[pixid].GetPedestalRms() * fNumHiGainSamples;
420
421 pix.SetPedestal(ped,prms);
422
423 if (TESTBIT(fFlags,kUseTFits))
424 pix.FitT();
425
426 if (!pix.FitQ())
427 continue;
428
429 }
430 }
431
432 fCalibrations->SetReadyToSave();
433
434 if (GetNumExecutions()==0)
435 return kTRUE;
436
437 *fLog << endl;
438 *fLog << dec << setfill(' ') << fCosmics << " Events presumably cosmics" << endl;
439
440 return kTRUE;
441}
Note: See TracBrowser for help on using the repository browser.