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

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