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

Last change on this file since 2929 was 2904, 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
266 fNumHiGainSamples = fSignals->GetNumUsedHiGainFADCSlices();
267 fNumLoGainSamples = fSignals->GetNumUsedLoGainFADCSlices();
268 fSqrtHiGainSamples = TMath::Sqrt((Float_t)fNumHiGainSamples);
269
270 UInt_t npixels = cam->GetNumPixels();
271
272 fCalibrations->InitSize(npixels);
273
274 for (UInt_t i=0; i<npixels; i++)
275 {
276
277 MCalibrationPix &pix = (*fCalibrations)[i];
278 pix.DefinePixId(i);
279 MHCalibrationPixel *hist = pix.GetHist();
280
281 hist->SetTimeFitRangesHiGain(fSignals->GetFirstUsedSliceHiGain(),
282 fSignals->GetLastUsedSliceHiGain());
283 hist->SetTimeFitRangesLoGain(fSignals->GetFirstUsedSliceLoGain(),
284 fSignals->GetLastUsedSliceLoGain());
285
286 if (!TESTBIT(fFlags,kUseQualityChecks))
287 pix.SetExcludeQualityCheck();
288
289 }
290
291 //
292 // Look for file to exclude pixels from analysis
293 //
294 if (!fExcludedPixelsFile.IsNull())
295 {
296
297 fExcludedPixelsFile = gSystem->ExpandPathName(fExcludedPixelsFile.Data());
298
299 //
300 // Initialize reading the file
301 //
302 ifstream in(fExcludedPixelsFile.Data(),ios::in);
303
304 if (in)
305 {
306 *fLog << inf << "Use excluded pixels from file: '" << fExcludedPixelsFile.Data() << "'" << endl;
307 //
308 // Read the file and count the number of entries
309 //
310 UInt_t pixel = 0;
311
312 while (++fNumExcludedPixels)
313 {
314
315 in >> pixel;
316
317 if (!in.good())
318 break;
319 //
320 // Check for out of range
321 //
322 if (pixel > npixels)
323 {
324 *fLog << warn << "WARNING: To be excluded pixel: " << pixel
325 << " is out of range " << endl;
326 continue;
327 }
328 //
329 // Exclude pixel
330 //
331 MCalibrationPix &pix = (*fCalibrations)[pixel];
332 pix.SetExcluded();
333
334 *fLog << GetDescriptor() << inf << ": Exclude Pixel: " << pixel << endl;
335
336 }
337
338 if (--fNumExcludedPixels == 0)
339 *fLog << warn << "WARNING: File '" << fExcludedPixelsFile.Data()
340 << "'" << " is empty " << endl;
341 else
342 fCalibrations->SetNumPixelsExcluded(fNumExcludedPixels);
343
344 }
345 else
346 *fLog << warn << dbginf << "Cannot open file '" << fExcludedPixelsFile.Data() << "'" << endl;
347 }
348
349 return kTRUE;
350}
351
352
353// --------------------------------------------------------------------------
354//
355// Calculate the integral of the FADC time slices and store them as a new
356// pixel in the MCerPhotEvt container.
357//
358Int_t MCalibrationCalc::Process()
359{
360
361 //
362 // Initialize pointers to blind pixel, PIN Diode and individual pixels
363 //
364 MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
365 MCalibrationPINDiode &pindiode = *(fCalibrations->GetPINDiode());
366
367 MRawEvtPixelIter pixel(fRawEvt);
368
369 //
370 // Perform cosmics cut
371 //
372 if (TESTBIT(fFlags,kUseCosmicsRejection))
373 {
374
375 Int_t cosmicpix = 0;
376
377 //
378 // Create a first loop to sort out the cosmics ...
379 //
380 // This is a very primitive check for the number of cosmicpixs
381 // The cut will be applied in the fit, but for the blind pixel,
382 // we need to remove this event
383 //
384 // FIXME: In the future need a much more sophisticated one!!!
385 //
386
387 while (pixel.Next())
388 {
389
390 const UInt_t pixid = pixel.GetPixelId();
391
392 MExtractedSignalPix &sig = (*fSignals)[pixid];
393 MPedestalPix &ped = (*fPedestals)[pixid];
394 Float_t pedrms = ped.GetPedestalRms()*fSqrtHiGainSamples;
395 Float_t sumhi = sig.GetExtractedSignalHiGain();
396
397 //
398 // We consider a pixel as presumably due to cosmics
399 // if its sum of FADC slices is lower than 3 pedestal RMS
400 //
401 if (sumhi < 3.*pedrms )
402 cosmicpix++;
403 }
404
405 //
406 // If the camera contains more than 230
407 // (this is the number of outer pixels plus about 50 inner ones)
408 // presumed pixels due to cosmics, then the event is discarted.
409 // This procedure is more or less equivalent to keeping only events
410 // with at least 350 pixels with high signals.
411 //
412 if (cosmicpix > 230.)
413 {
414 fCosmics++;
415 return kCONTINUE;
416 }
417
418 pixel.Reset();
419 }
420
421 fEvents++;
422
423 //
424 // Create a (second) loop to do fill the calibration histograms
425 //
426
427 while (pixel.Next())
428 {
429
430 const UInt_t pixid = pixel.GetPixelId();
431
432 MCalibrationPix &pix = (*fCalibrations)[pixid];
433
434 if (pix.IsExcluded())
435 continue;
436
437 MExtractedSignalPix &sig = (*fSignals)[pixid];
438
439 const Float_t sumhi = sig.GetExtractedSignalHiGain();
440 const Float_t sumlo = sig.GetExtractedSignalLoGain();
441
442 Double_t mtime = 0.;
443
444 if (TESTBIT(fFlags,kUseTimes))
445 mtime = (*fArrivalTime)[pixid];
446
447 switch(pixid)
448 {
449
450 case gkCalibrationBlindPixelId:
451
452 if (!blindpixel.FillCharge(sumhi))
453 *fLog << warn <<
454 "Overflow or Underflow occurred filling Blind Pixel sum = " << sumhi << endl;
455
456 if (TESTBIT(fFlags,kUseTimes))
457 {
458 if (!blindpixel.FillTime(mtime))
459 *fLog << warn <<
460 "Overflow or Underflow occurred filling Blind Pixel time = " << mtime << endl;
461 }
462
463 if (!blindpixel.FillRChargevsTime(sumhi,fEvents))
464 *fLog << warn <<
465 "Overflow or Underflow occurred filling Blind Pixel eventnr = " << fEvents << endl;
466 break;
467
468 case gkCalibrationPINDiodeId:
469
470 if (!pindiode.FillCharge(sumhi))
471 *fLog << warn <<
472 "Overflow or Underflow occurred filling PINDiode: sum = " << sumhi << endl;
473
474 if (TESTBIT(fFlags,kUseTimes))
475 {
476 if (!pindiode.FillTime(mtime))
477 *fLog << warn <<
478 "Overflow or Underflow occurred filling PINDiode: time = " << mtime << endl;
479 }
480
481 if (!pindiode.FillRChargevsTime(sumhi,fEvents))
482 *fLog << warn <<
483 "Overflow or Underflow occurred filling PINDiode: eventnr = " << fEvents << endl;
484 break;
485
486 default:
487
488 pix.FillChargesInGraph(sumhi,sumlo);
489
490 if (!pix.FillRChargevsTimeLoGain(sumlo,fEvents))
491 *fLog << warn << "Could not fill Lo Gain Charge vs. EvtNr of pixel: "
492 << pixid << " signal = " << sumlo << " event Nr: " << fEvents << endl;
493
494 if (!pix.FillRChargevsTimeHiGain(sumhi,fEvents))
495 *fLog << warn << "Could not fill Hi Gain Charge vs. EvtNr of pixel: "
496 << pixid << " signal = " << sumhi << " event Nr: " << fEvents << endl;
497
498 if (sig.IsLoGainUsed())
499 {
500
501 if (!pix.FillChargeLoGain(sumlo))
502 *fLog << warn << "Could not fill Lo Gain Charge of pixel: " << pixid
503 << " signal = " << sumlo << endl;
504
505 if (TESTBIT(fFlags,kUseTimes))
506 {
507 if (!pix.FillTimeLoGain(mtime))
508 *fLog << warn << "Could not fill Lo Gain Time of pixel: "
509 << pixid << " time = " << mtime << endl;
510 }
511 } /* if (sig.IsLoGainUsed()) */
512 else
513 {
514 if (!pix.FillChargeHiGain(sumhi))
515 *fLog << warn << "Could not fill Hi Gain Charge of pixel: " << pixid
516 << " signal = " << sumhi << endl;
517
518 if (TESTBIT(fFlags,kUseTimes))
519 {
520 if (!pix.FillTimeHiGain(mtime))
521 *fLog << warn << "Could not fill Hi Gain Time of pixel: "
522 << pixid << " time = " << mtime << endl;
523 }
524 } /* else (sig.IsLoGainUsed()) */
525 break;
526
527 } /* switch(pixid) */
528
529 } /* while (pixel.Next()) */
530
531 return kTRUE;
532}
533
534Int_t MCalibrationCalc::PostProcess()
535{
536
537 *fLog << inf << endl;
538
539 if (fEvents == 0)
540 {
541 *fLog << err << GetDescriptor()
542 << ": This run contains only cosmics or pedestals, "
543 << "cannot find events with more than 350 illuminated pixels. " << endl;
544 return kFALSE;
545 }
546
547 if (fEvents < fCosmics)
548 *fLog << warn << GetDescriptor()
549 << ": WARNING: Run contains more cosmics or pedestals than calibration events " << endl;
550
551
552 *fLog << inf << GetDescriptor() << ": Cut Histogram Edges" << endl;
553
554 //
555 // Cut edges to make fits and viewing of the hists easier
556 //
557 fCalibrations->CutEdges();
558
559 //
560 // Fit the blind pixel
561 //
562 if (TESTBIT(fFlags,kUseBlindPixelFit))
563 {
564 //
565 // Get pointer to blind pixel
566 //
567 MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
568
569 *fLog << inf << GetDescriptor() << ": Fitting the Blind Pixel" << endl;
570
571 if (!blindpixel.FitCharge())
572 {
573 *fLog << err << dbginf << "Could not fit the blind pixel! " << endl;
574 blindpixel.DrawClone();
575 return kFALSE;
576 }
577
578 blindpixel.DrawClone();
579 }
580 else
581 *fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Fit " << endl;
582
583
584 *fLog << inf << GetDescriptor() << ": Fitting the Normal Pixels" << endl;
585
586 //
587 // loop over the pedestal events and check if we have calibration
588 //
589 for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
590 {
591
592 MCalibrationPix &pix = (*fCalibrations)[pixid];
593
594 //
595 // get the pedestals
596 //
597 const Float_t ped = (*fPedestals)[pixid].GetPedestal() * fNumHiGainSamples;
598 const Float_t prms = (*fPedestals)[pixid].GetPedestalRms() * fSqrtHiGainSamples;
599
600 //
601 // set them in the calibration camera
602 //
603 pix.SetPedestal(ped,prms);
604
605
606 //
607 // Check if the pixel has been excluded from the fits
608 //
609 if (pix.IsExcluded())
610 continue;
611
612 //
613 // perform the Gauss fits to the charges
614 //
615 pix.FitCharge();
616
617 //
618 // Perform the Gauss fits to the arrival times
619 //
620 if (TESTBIT(fFlags,kUseTimes))
621 pix.FitTime();
622
623 }
624
625 if (TESTBIT(fFlags,kUseBlindPixelFit))
626 {
627
628 if (!fCalibrations->CalcNumPhotInsidePlexiglass())
629 {
630 *fLog << err
631 << "Could not calculate the number of photons from the blind pixel " << endl;
632 *fLog << err
633 << "You can try to calibrate using the MCalibrationCalc::SkipBlindPixelFit()" << endl;
634 return kFALSE;
635 }
636 }
637 else
638 *fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Fit " << endl;
639
640 fCalibrations->SetReadyToSave();
641
642 if (GetNumExecutions()==0)
643 return kTRUE;
644
645 *fLog << inf << endl;
646 *fLog << dec << setfill(' ') << fCosmics << " Events presumably cosmics" << endl;
647
648 return kTRUE;
649}
650
Note: See TracBrowser for help on using the repository browser.