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

Last change on this file since 2932 was 2931, checked in by gaug, 21 years ago
*** empty log message ***
File size: 19.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// Task to calculate the calibration conversion factors from the FADC
30// time slices. The integrated time slices have to be delivered by an
31// MExtractedSignalCam. The pedestals by an MPedestalCam and possibly
32// arrival times from MArrivalTime
33//
34// The output container MCalibrationCam holds one entry of type MCalibrationPix
35// for every pixel. It is filled in the following way:
36//
37// ProProcess: Search for MPedestalCam, MExtractedSignalCam
38// Initialize MCalibrationCam
39// Initialize MArrivalTime if exists
40// Initialize pulser light wavelength
41//
42// ReInit: MCalibrationCam::InitSize(NumPixels) is called which allocates
43// memory in a TClonesArray of type MCalibrationPix
44// Initialize number of used FADC slices
45// Optionally exclude pixels from calibration
46//
47// Process: Optionally, a cut on cosmics can be performed
48//
49// Every MCalibrationPix holds a histogram class,
50// MHCalibrationPixel which itself hold histograms of type:
51// HCharge(npix) (distribution of summed FADC time slice
52// entries)
53// HTime(npix) (distribution of position of maximum)
54// HChargevsN(npix) (distribution of charges vs. event number.
55//
56// PostProcess: All histograms HCharge(npix) are fitted to a Gaussian
57// All histograms HTime(npix) are fitted to a Gaussian
58// The histogram HBlindPixelCharge (blind pixel) is fitted to
59// a single PhE fit
60//
61// The histograms of the PIN Diode are fitted to Gaussians
62//
63// Fits can be excluded via the commands:
64// MalibrationCam::SkipTimeFits() (skip all time fits)
65// MalibrationCam::SkipBlindPixelFits() (skip all blind
66// pixel fits)
67// MalibrationCam::SkipPinDiodeFits() (skip all PIN Diode
68// fits)
69//
70// Input Containers:
71// MRawEvtData
72// MPedestalCam
73// MExtractedSignalCam
74//
75// Output Containers:
76// MCalibrationCam
77//
78//////////////////////////////////////////////////////////////////////////////
79#include "MCalibrationCalc.h"
80
81// FXIME: Usage of fstream is a preliminary workaround!
82#include <fstream>
83
84// FXIME: This has to be removed!!!! (YES, WHEN WE HAVE ACCESS TO THE DATABASE!!!!!)
85#include "MCalibrationConfig.h"
86
87#include <TSystem.h>
88
89#include "MLog.h"
90#include "MLogManip.h"
91
92#include "MParList.h"
93
94#include "MGeomCam.h"
95#include "MRawRunHeader.h"
96#include "MRawEvtPixelIter.h"
97
98#include "MPedestalCam.h"
99#include "MPedestalPix.h"
100
101#include "MCalibrationCam.h"
102#include "MCalibrationPix.h"
103
104#include "MExtractedSignalCam.h"
105#include "MExtractedSignalPix.h"
106
107#include "MCalibrationBlindPix.h"
108#include "MCalibrationPINDiode.h"
109
110#include "MArrivalTime.h"
111
112ClassImp(MCalibrationCalc);
113
114using namespace std;
115
116// --------------------------------------------------------------------------
117//
118// Default constructor.
119//
120MCalibrationCalc::MCalibrationCalc(const char *name, const char *title)
121 : fPedestals(NULL), fCalibrations(NULL), fSignals(NULL),
122 fRawEvt(NULL), fRunHeader(NULL), fArrivalTime(NULL), fEvtTime(NULL),
123 fEvents(0), fCosmics(0),
124 fNumHiGainSamples(0), fNumLoGainSamples(0), fConversionHiLo(0.),
125 fNumExcludedPixels(0),
126 fColor(kEBlue)
127{
128
129 fName = name ? name : "MCalibrationCalc";
130 fTitle = title ? title : "Task to calculate the calibration constants and MCalibrationCam ";
131
132 AddToBranchList("MRawEvtData.fHiGainPixId");
133 AddToBranchList("MRawEvtData.fLoGainPixId");
134 AddToBranchList("MRawEvtData.fHiGainFadcSamples");
135 AddToBranchList("MRawEvtData.fLoGainFadcSamples");
136
137 SETBIT(fFlags, kUseTimes);
138 SETBIT(fFlags, kUseBlindPixelFit);
139 SETBIT(fFlags, kUsePinDiodeFit);
140 SETBIT(fFlags, kUseCosmicsRejection);
141 SETBIT(fFlags, kUseQualityChecks);
142
143}
144
145MCalibrationBlindPix *MCalibrationCalc::GetBlindPixel() const
146{
147 return fCalibrations->GetBlindPixel();
148}
149
150MCalibrationPINDiode *MCalibrationCalc::GetPINDiode() const
151{
152 return fCalibrations->GetPINDiode();
153}
154
155// --------------------------------------------------------------------------
156//
157// The PreProcess searches for the following input containers:
158// - MRawEvtData
159// - MPedestalCam
160//
161// The following output containers are also searched and created if
162// they were not found:
163//
164// - MHCalibrationBlindPixel
165// - MCalibrationCam
166//
167// The following output containers are only searched, but not created
168//
169// - MArrivaltime
170// - MTime
171//
172Int_t MCalibrationCalc::PreProcess(MParList *pList)
173{
174
175 fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
176 if (!fRawEvt)
177 {
178 *fLog << err << dbginf << "MRawEvtData not found... aborting." << endl;
179 return kFALSE;
180 }
181
182 const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
183 if (!runheader)
184 *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl;
185 else
186 if (runheader->GetRunType() == kRTMonteCarlo)
187 {
188 return kTRUE;
189 }
190
191 fCalibrations = (MCalibrationCam*)pList->FindCreateObj("MCalibrationCam");
192 if (!fCalibrations)
193 {
194 *fLog << err << dbginf << "MCalibrationCam could not be created ... aborting." << endl;
195 return kFALSE;
196 }
197
198 fArrivalTime = (MArrivalTime*)pList->FindObject("MArrivalTime");
199
200 if (!fArrivalTime)
201 CLRBIT(fFlags,kUseTimes);
202
203 fEvtTime = (MTime*)pList->FindObject("MTime");
204
205 switch (fColor)
206 {
207 case kEBlue:
208 fCalibrations->SetColor(MCalibrationCam::kECBlue);
209 break;
210 case kEGreen:
211 fCalibrations->SetColor(MCalibrationCam::kECGreen);
212 break;
213 case kEUV:
214 fCalibrations->SetColor(MCalibrationCam::kECUV);
215 break;
216 case kECT1:
217 fCalibrations->SetColor(MCalibrationCam::kECCT1);
218 break;
219 default:
220 fCalibrations->SetColor(MCalibrationCam::kECCT1);
221 }
222
223 fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
224 if (!fPedestals)
225 {
226 *fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl;
227 return kFALSE;
228 }
229
230
231 fSignals = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
232 if (!fSignals)
233 {
234 *fLog << err << dbginf << "Cannot find MExtractedSignalCam ... aborting" << endl;
235 return kFALSE;
236 }
237
238 return kTRUE;
239}
240
241
242// --------------------------------------------------------------------------
243//
244// The ReInit searches for the following input containers:
245// - MRawRunHeader
246//
247Bool_t MCalibrationCalc::ReInit(MParList *pList )
248{
249
250 fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
251 if (!fRunHeader)
252 {
253 *fLog << err << dbginf << ": MRawRunHeader not found... aborting." << endl;
254 return kFALSE;
255 }
256
257
258 MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
259 if (!cam)
260 {
261 *fLog << err << GetDescriptor() << ": No MGeomCam found... aborting." << endl;
262 return kFALSE;
263 }
264
265 fNumHiGainSamples = fSignals->GetNumUsedHiGainFADCSlices();
266 fNumLoGainSamples = fSignals->GetNumUsedLoGainFADCSlices();
267 fSqrtHiGainSamples = TMath::Sqrt((Float_t)fNumHiGainSamples);
268
269 UInt_t npixels = cam->GetNumPixels();
270
271 fCalibrations->InitSize(npixels);
272
273 for (UInt_t i=0; i<npixels; i++)
274 {
275
276 MCalibrationPix &pix = (*fCalibrations)[i];
277 pix.DefinePixId(i);
278 MHCalibrationPixel *hist = pix.GetHist();
279
280 hist->SetTimeFitRangesHiGain(fSignals->GetFirstUsedSliceHiGain(),
281 fSignals->GetLastUsedSliceHiGain());
282 hist->SetTimeFitRangesLoGain(fSignals->GetFirstUsedSliceLoGain(),
283 fSignals->GetLastUsedSliceLoGain());
284
285 if (!TESTBIT(fFlags,kUseQualityChecks))
286 pix.SetExcludeQualityCheck();
287
288 }
289
290 //
291 // Look for file to exclude pixels from analysis
292 //
293 if (!fExcludedPixelsFile.IsNull())
294 {
295
296 fExcludedPixelsFile = gSystem->ExpandPathName(fExcludedPixelsFile.Data());
297
298 //
299 // Initialize reading the file
300 //
301 ifstream in(fExcludedPixelsFile.Data(),ios::in);
302
303 if (in)
304 {
305 *fLog << inf << "Use excluded pixels from file: '" << fExcludedPixelsFile.Data() << "'" << endl;
306 //
307 // Read the file and count the number of entries
308 //
309 UInt_t pixel = 0;
310
311 while (++fNumExcludedPixels)
312 {
313
314 in >> pixel;
315
316 if (!in.good())
317 break;
318 //
319 // Check for out of range
320 //
321 if (pixel > npixels)
322 {
323 *fLog << warn << "WARNING: To be excluded pixel: " << pixel
324 << " is out of range " << endl;
325 continue;
326 }
327 //
328 // Exclude pixel
329 //
330 MCalibrationPix &pix = (*fCalibrations)[pixel];
331 pix.SetExcluded();
332
333 *fLog << GetDescriptor() << inf << ": Exclude Pixel: " << pixel << endl;
334
335 }
336
337 if (--fNumExcludedPixels == 0)
338 *fLog << warn << "WARNING: File '" << fExcludedPixelsFile.Data()
339 << "'" << " is empty " << endl;
340 else
341 fCalibrations->SetNumPixelsExcluded(fNumExcludedPixels);
342
343 }
344 else
345 *fLog << warn << dbginf << "Cannot open file '" << fExcludedPixelsFile.Data() << "'" << endl;
346 }
347
348 return kTRUE;
349}
350
351
352// --------------------------------------------------------------------------
353//
354// Calculate the integral of the FADC time slices and store them as a new
355// pixel in the MCerPhotEvt container.
356//
357Int_t MCalibrationCalc::Process()
358{
359
360 //
361 // Initialize pointers to blind pixel, PIN Diode and individual pixels
362 //
363 MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
364 MCalibrationPINDiode &pindiode = *(fCalibrations->GetPINDiode());
365
366 MRawEvtPixelIter pixel(fRawEvt);
367
368 //
369 // Perform cosmics cut
370 //
371 if (TESTBIT(fFlags,kUseCosmicsRejection))
372 {
373
374 Int_t cosmicpix = 0;
375
376 //
377 // Create a first loop to sort out the cosmics ...
378 //
379 // This is a very primitive check for the number of cosmicpixs
380 // The cut will be applied in the fit, but for the blind pixel,
381 // we need to remove this event
382 //
383 // FIXME: In the future need a much more sophisticated one!!!
384 //
385
386 while (pixel.Next())
387 {
388
389 const UInt_t pixid = pixel.GetPixelId();
390
391 MExtractedSignalPix &sig = (*fSignals)[pixid];
392 MPedestalPix &ped = (*fPedestals)[pixid];
393 Float_t pedrms = ped.GetPedestalRms()*fSqrtHiGainSamples;
394 Float_t sumhi = sig.GetExtractedSignalHiGain();
395
396 //
397 // We consider a pixel as presumably due to cosmics
398 // if its sum of FADC slices is lower than 3 pedestal RMS
399 //
400 if (sumhi < 3.*pedrms )
401 cosmicpix++;
402 }
403
404 //
405 // If the camera contains more than 230
406 // (this is the number of outer pixels plus about 50 inner ones)
407 // presumed pixels due to cosmics, then the event is discarted.
408 // This procedure is more or less equivalent to keeping only events
409 // with at least 350 pixels with high signals.
410 //
411 if (cosmicpix > 230.)
412 {
413 fCosmics++;
414 return kCONTINUE;
415 }
416
417 pixel.Reset();
418 }
419
420 fEvents++;
421
422 //
423 // Create a (second) loop to do fill the calibration histograms
424 //
425
426 while (pixel.Next())
427 {
428
429 const UInt_t pixid = pixel.GetPixelId();
430
431 MCalibrationPix &pix = (*fCalibrations)[pixid];
432
433 if (pix.IsExcluded())
434 continue;
435
436 MExtractedSignalPix &sig = (*fSignals)[pixid];
437
438 const Float_t sumhi = sig.GetExtractedSignalHiGain();
439 const Float_t sumlo = sig.GetExtractedSignalLoGain();
440
441 Double_t mtime = 0.;
442
443 if (TESTBIT(fFlags,kUseTimes))
444 mtime = (*fArrivalTime)[pixid];
445
446 switch(pixid)
447 {
448
449 case gkCalibrationBlindPixelId:
450
451 if (!blindpixel.FillCharge(sumhi))
452 *fLog << warn <<
453 "Overflow or Underflow occurred filling Blind Pixel sum = " << sumhi << endl;
454
455 if (TESTBIT(fFlags,kUseTimes))
456 {
457 if (!blindpixel.FillTime(mtime))
458 *fLog << warn <<
459 "Overflow or Underflow occurred filling Blind Pixel time = " << mtime << endl;
460 }
461
462 if (!blindpixel.FillRChargevsTime(sumhi,fEvents))
463 *fLog << warn <<
464 "Overflow or Underflow occurred filling Blind Pixel eventnr = " << fEvents << endl;
465 break;
466
467 case gkCalibrationPINDiodeId:
468
469 if (!pindiode.FillCharge(sumhi))
470 *fLog << warn <<
471 "Overflow or Underflow occurred filling PINDiode: sum = " << sumhi << endl;
472
473 if (TESTBIT(fFlags,kUseTimes))
474 {
475 if (!pindiode.FillTime(mtime))
476 *fLog << warn <<
477 "Overflow or Underflow occurred filling PINDiode: time = " << mtime << endl;
478 }
479
480 if (!pindiode.FillRChargevsTime(sumhi,fEvents))
481 *fLog << warn <<
482 "Overflow or Underflow occurred filling PINDiode: eventnr = " << fEvents << endl;
483 break;
484
485 default:
486
487 pix.FillChargesInGraph(sumhi,sumlo);
488
489 if (!pix.FillRChargevsTimeLoGain(sumlo,fEvents))
490 *fLog << warn << "Could not fill Lo Gain Charge vs. EvtNr of pixel: "
491 << pixid << " signal = " << sumlo << " event Nr: " << fEvents << endl;
492
493 if (!pix.FillRChargevsTimeHiGain(sumhi,fEvents))
494 *fLog << warn << "Could not fill Hi Gain Charge vs. EvtNr of pixel: "
495 << pixid << " signal = " << sumhi << " event Nr: " << fEvents << endl;
496
497 if (sig.IsLoGainUsed())
498 {
499
500 if (!pix.FillChargeLoGain(sumlo))
501 *fLog << warn << "Could not fill Lo Gain Charge of pixel: " << pixid
502 << " signal = " << sumlo << endl;
503
504 if (TESTBIT(fFlags,kUseTimes))
505 {
506 if (!pix.FillTimeLoGain(mtime))
507 *fLog << warn << "Could not fill Lo Gain Time of pixel: "
508 << pixid << " time = " << mtime << endl;
509 }
510 } /* if (sig.IsLoGainUsed()) */
511 else
512 {
513 if (!pix.FillChargeHiGain(sumhi))
514 *fLog << warn << "Could not fill Hi Gain Charge of pixel: " << pixid
515 << " signal = " << sumhi << endl;
516
517 if (TESTBIT(fFlags,kUseTimes))
518 {
519 if (!pix.FillTimeHiGain(mtime))
520 *fLog << warn << "Could not fill Hi Gain Time of pixel: "
521 << pixid << " time = " << mtime << endl;
522 }
523 } /* else (sig.IsLoGainUsed()) */
524 break;
525
526 } /* switch(pixid) */
527
528 } /* while (pixel.Next()) */
529
530 return kTRUE;
531}
532
533Int_t MCalibrationCalc::PostProcess()
534{
535
536 *fLog << inf << endl;
537
538 if (fEvents == 0)
539 {
540 *fLog << err << GetDescriptor()
541 << ": This run contains only cosmics or pedestals, "
542 << "cannot find events with more than 350 illuminated pixels. " << endl;
543 return kFALSE;
544 }
545
546 if (fEvents < fCosmics)
547 *fLog << warn << GetDescriptor()
548 << ": WARNING: Run contains more cosmics or pedestals than calibration events " << endl;
549
550
551 *fLog << inf << GetDescriptor() << ": Cut Histogram Edges" << endl;
552
553 //
554 // Cut edges to make fits and viewing of the hists easier
555 //
556 fCalibrations->CutEdges();
557
558 //
559 // Fit the blind pixel
560 //
561 if (TESTBIT(fFlags,kUseBlindPixelFit))
562 {
563 //
564 // Get pointer to blind pixel
565 //
566 MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
567
568 *fLog << inf << GetDescriptor() << ": Fitting the Blind Pixel" << endl;
569
570 if (!blindpixel.FitCharge())
571 {
572 *fLog << err << dbginf << "Could not fit the blind pixel! " << endl;
573 blindpixel.DrawClone();
574 return kFALSE;
575 }
576
577 blindpixel.DrawClone();
578 }
579 else
580 *fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Fit " << endl;
581
582
583 *fLog << inf << GetDescriptor() << ": Fitting the Normal Pixels" << endl;
584
585 //
586 // loop over the pedestal events and check if we have calibration
587 //
588 for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
589 {
590
591 MCalibrationPix &pix = (*fCalibrations)[pixid];
592
593 //
594 // get the pedestals
595 //
596 const Float_t ped = (*fPedestals)[pixid].GetPedestal() * fNumHiGainSamples;
597 const Float_t prms = (*fPedestals)[pixid].GetPedestalRms() * fSqrtHiGainSamples;
598
599 //
600 // set them in the calibration camera
601 //
602 pix.SetPedestal(ped,prms);
603
604
605 //
606 // Check if the pixel has been excluded from the fits
607 //
608 if (pix.IsExcluded())
609 continue;
610
611 //
612 // perform the Gauss fits to the charges
613 //
614 pix.FitCharge();
615
616 //
617 // Perform the Gauss fits to the arrival times
618 //
619 if (TESTBIT(fFlags,kUseTimes))
620 pix.FitTime();
621
622 }
623
624 if (TESTBIT(fFlags,kUseBlindPixelFit))
625 {
626
627 if (!fCalibrations->CalcNumPhotInsidePlexiglass())
628 {
629 *fLog << err
630 << "Could not calculate the number of photons from the blind pixel " << endl;
631 *fLog << err
632 << "You can try to calibrate using the MCalibrationCalc::SkipBlindPixelFit()" << endl;
633 return kFALSE;
634 }
635 }
636 else
637 *fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Fit " << endl;
638
639 fCalibrations->SetReadyToSave();
640
641 if (GetNumExecutions()==0)
642 return kTRUE;
643
644 *fLog << inf << endl;
645 *fLog << dec << setfill(' ') << fCosmics << " Events presumably cosmics" << endl;
646
647 return kTRUE;
648}
649
Note: See TracBrowser for help on using the repository browser.