/* ======================================================================== *\ ! ! * ! * This file is part of MARS, the MAGIC Analysis and Reconstruction ! * Software. It is distributed to you in the hope that it can be a useful ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. ! * It is distributed WITHOUT ANY WARRANTY. ! * ! * Permission to use, copy, modify and distribute this software and its ! * documentation for any purpose is hereby granted without fee, ! * provided that the above copyright notice appear in all copies and ! * that both that copyright notice and this permission notice appear ! * in supporting documentation. It is provided "as is" without express ! * or implied warranty. ! * ! ! ! Author(s): Markus Gaug 09/2003 ! ! Copyright: MAGIC Software Development, 2000-2001 ! ! \* ======================================================================== */ ////////////////////////////////////////////////////////////////////////////// // // // MCalibrationCalc // // // // This is a task which calculates the number of photons from the FADC // // time slices. At the moment it integrates simply the FADC values. // // // // Input Containers: // // MRawEvtData // // // // Output Containers: // // MCalibrationCam // // // // // // The class MCalibrationCam hold one entry of type MCalibrationPix for // // every pixel. It is filled in the following way: // // PreParocess: MalibrationCam::InitSize(577) is called which allocates // // memory in an TClonesArray of type MCalibrationPix and // // all pointers to NULL. // // // // Process: The NULL pointer is tested on every pixel via // // MalibrationCam::IsPixelUsed(npix). // // // // In case, IsPixelUsed returns NULL, // // MalibrationCam::AddPixel(npix) is invoked which creates a // // new MCalibrationPix(npix) in the npix's entry // // of the TClonesArray. // // // // Every MCalibrationPix holds a histogram class, // // MHCalibrationPixel which itself hold histograms of type: // // HQ(npix) (distribution of summed FADC time slice entries) // // HT(npix) (distribution of position of maximum) // // HQvsN(npix) (distribution of charges vs. event number. // // // // PostProcess: All histograms HQ(npix) are fitted to a Gaussian // // All histograms HT(npix) are fitted to a Gaussian // // The histogram HBPQ (blind pixel) is fitted to a single // // PhE fit // // The histogram HBPT (blind pixel) is fitted to a Gaussian // // The histograms of the PIN Diode are fitted to Gaussians // // // // Fits can be excluded via the commands: // // MalibrationCam::SetSkipTFits() (skip all time fits) // // MalibrationCam::SetSkipBPFits() (skip all blind pixel fits) // // MalibrationCam::SetSkipPDFits() (skip all PIN Diode fits) // // // ////////////////////////////////////////////////////////////////////////////// #include "MCalibrationCalc.h" #include "MCalibrationConfig.h" #include "MCalibrationFits.h" #include "MCalibrationCam.h" #include "MCalibrationPix.h" #include "MCalibrationBlindPix.h" #include "MCalibrationPINDiode.h" #include "MPedestalCam.h" #include "MPedestalPix.h" #include "MLog.h" #include "MLogManip.h" #include "MParList.h" #include "MH.h" #include "MRawRunHeader.h" #include "MRawEvtData.h" // MRawEvtData::GetNumPixels #include "MRawEvtPixelIter.h" #include "MTime.h" ClassImp(MCalibrationCalc); using namespace std; // -------------------------------------------------------------------------- // // Default constructor. b is the number of slices before the maximum slice, // a the number of slices behind the maximum slice which is taken as signal. // MCalibrationCalc::MCalibrationCalc(const char *name, const char *title) : fColor(kEBlue) { fName = name ? name : "MCalibrationCalc"; fTitle = title ? title : "Task to calculate the calibration constants and MCalibrationCam "; AddToBranchList("MRawEvtData.fHiGainPixId"); AddToBranchList("MRawEvtData.fLoGainPixId"); AddToBranchList("MRawEvtData.fHiGainFadcSamples"); AddToBranchList("MRawEvtData.fLoGainFadcSamples"); SETBIT(fFlags, kUseTFits); SETBIT(fFlags, kUseBPFit); SETBIT(fFlags, kUsePDFit); } // -------------------------------------------------------------------------- // // The PreProcess searches for the following input containers: // - MRawEvtData // - MPedestalCam // // The following output containers are also searched and created if // they were not found: // // - MHCalibrationBlindPixel // - MCalibrationCam // - MTime // Int_t MCalibrationCalc::PreProcess(MParList *pList) { fOverFlow = 0; fNrEvents = 0; fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData"); if (!fRawEvt) { *fLog << dbginf << "MRawEvtData not found... aborting." << endl; return kFALSE; } const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); if (!runheader) *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl; else if (runheader->GetRunType() == kRTMonteCarlo) { return kTRUE; } fCalibrations = (MCalibrationCam*)pList->FindCreateObj("MCalibrationCam"); if (!fCalibrations) { *fLog << err << dbginf << "MCalibrationCam could not be created ... aborting." << endl; return kFALSE; } // // fRawEvt->GetNumPixels() does not work !!!!!!!!!! // fCalibrations->InitSize(577); switch (fColor) { case kEBlue: fCalibrations->SetColor(MCalibrationCam::kECBlue); break; case kEGreen: fCalibrations->SetColor(MCalibrationCam::kECGreen); break; case kEUV: fCalibrations->SetColor(MCalibrationCam::kECUV); } fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam"); if (!fPedestals) { *fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl; return kFALSE; } // MTime does not work!!!! // fEvtTime = (MTime*)pList->FindObject("MRawEvtTime"); // if (!fEvtTime) // { // *fLog << dbginf << "MTime could not be created ... ¡!¡!¡!¡!" << endl; // } // fEvtTime->SetTime(0.); return kTRUE; } // -------------------------------------------------------------------------- // // The ReInit searches for the following input containers: // - MRawRunHeader // Bool_t MCalibrationCalc::ReInit(MParList *pList ) { fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); if (!fRunHeader) { *fLog << dbginf << "MRawRunHeader not found... aborting." << endl; return kFALSE; } fNumHiGainSamples = fRunHeader->GetNumSamplesHiGain(); fNumLoGainSamples = fRunHeader->GetNumSamplesLoGain(); return kTRUE; } // -------------------------------------------------------------------------- // // Calculate the integral of the FADC time slices and store them as a new // pixel in the MCerPhotEvt container. // Int_t MCalibrationCalc::Process() { fNrEvents++; MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel()); MCalibrationPINDiode &pindiode = *(fCalibrations->GetPINDiode()); MRawEvtPixelIter pixel(fRawEvt); while (pixel.Next()) { UShort_t sat = 0; const Int_t pixid = pixel.GetPixelId(); if (!fCalibrations->IsPixelUsed(pixid)) if (!fCalibrations->AddPixel(pixid)) *fLog << err << dbginf << "MCalibrationPix(" << pixid << ") could not be created !!" << endl; Byte_t *ptr = pixel.GetHiGainSamples(); Byte_t mid = pixel.GetIdxMaxHiGainSample(); UInt_t max = pixel.GetMaxHiGainSample(); Int_t sum = (max > gkSaturationLimit // overflow ? ? sat++, pixel.GetSumLoGainSamples() // take Low Gain : pixel.GetSumHiGainSamples()); // no overflow if (sat) { ptr = pixel.GetLoGainSamples(); max = pixel.GetMaxLoGainSample(); mid = pixel.GetIdxMaxLoGainSample(); sum = (max == gkSaturationLimit // overflow of LoGain ??? -> GimmeABreak!!! ? fOverFlow++, gkLoGainOverFlow // OUCH (Florian was maybe right) : sum*gkConversionHiLo ); // OUFF (Florian was wrong) !! // *fLog << warn << "Warning: Saturation of HiGain reached in slice " << (int)mid << " !!! " << endl; // *fLog << warn << "Max = " << max << endl; if (fOverFlow) *fLog << err << dbginf << "Warning: Saturation of LoGain reached! " << err << dbginf << "sum = " << sum << endl; } // // sanity check (old MC files sometimes have pixids>577) // if (fPedestals && !fPedestals->CheckBounds(pixid)) { *fLog << inf << "Pixel ID larger than camera... skipping event." << endl; return kCONTINUE; } MPedestalPix &ped = (*fPedestals)[pixid]; MCalibrationPix &pix = (*fCalibrations)[pixid]; Float_t pedes = ped.GetPedestal(); // // FIXME: This is preliminary, we will change to pedestals per slice!!! // Assume pedestals per time slice ==> multiply with number of slices // pedes *= (sat ? fNumLoGainSamples : fNumHiGainSamples ); Float_t rsum = (float)sum - pedes; switch(pixid) { case gkCalibrationBlindPixelId: if (!blindpixel.FillQ(sum)) *fLog << warn << "Overflow or Underflow occurred filling Blind Pixel means = " << sum << endl; if (!blindpixel.FillT((int)mid)) *fLog << warn << "Overflow or Underflow occurred filling Blind Pixel time = " << (int)mid << endl; if (!blindpixel.FillRQvsT(rsum,fNrEvents)) *fLog << warn << "Overflow or Underflow occurred filling Blind Pixel eventnr = " << fNrEvents << endl; case gkCalibrationPINDiodeId: if (!pindiode.FillQ(sum)) *fLog << warn << "Overflow or Underflow occurred filling HQ: means = " << sum << endl; if (!pindiode.FillT((int)mid)) *fLog << warn << "Overflow or Underflow occurred filling HT: time = " << (int)mid << endl; if (!pindiode.FillRQvsT(rsum,fNrEvents)) *fLog << warn << "Overflow or Underflow occurred filling HQvsN: eventnr = " << fNrEvents << endl; default: if (!pix.FillQ(sum)) *fLog << warn << "Could not fill Q of pixel: " << pixid << " signal = " << sum << endl; // // Fill the reduced charge into the control histo for better visibility // if (!pix.FillRQvsT(rsum,fNrEvents)) *fLog << warn << "Could not fill red. Q vs. EvtNr of pixel: " << pixid << " signal = " << rsum << " event Nr: " << fNrEvents << endl; if (!pix.FillT((int)mid)) *fLog << warn << "Could not fill T of pixel: " << pixid << " time = " << (int)mid << endl; } /* switch(pixid) */ } /* while (pixel.Next()) */ fCalibrations->SetReadyToSave(); return kTRUE; } Int_t MCalibrationCalc::PostProcess() { *fLog << inf << endl; *fLog << GetDescriptor() << " Cut Histogram Edges" << endl; // // Cut edges to make fits and viewing of the hists easier // fCalibrations->CutEdges(); // // Get pointer to blind pixel // MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel()); *fLog << GetDescriptor() << " Fitting the Blind Pixel" << endl; // // Fit the blind pixel // if (TESTBIT(fFlags,kUseBPFit)) { if (!blindpixel.FitQ()) *fLog << err << dbginf << "Could not fit the blind pixel " << endl; if (!blindpixel.FitT()) *fLog << warn << "Could not the Times of the blind pixel " << endl; blindpixel.Draw(); } *fLog << GetDescriptor() << " Fitting the Normal Pixels" << endl; // // loop over the pedestal events and check if we have calibration // for (Int_t pixid=0; pixidGetSize(); pixid++) { if (fCalibrations->IsPixelUsed(pixid)) { MCalibrationPix &pix = (*fCalibrations)[pixid]; if (TESTBIT(fFlags,kUseTFits)) pix.FitT(); if (!pix.FitQ()) continue; const Float_t ped = (*fPedestals)[pixid].GetPedestal(); const Float_t prms = (*fPedestals)[pixid].GetPedestalRms(); pix.SetPedestal(ped,prms); } } fCalibrations->SetReadyToSave(); if (GetNumExecutions()==0 || fOverFlow==0) return kTRUE; *fLog << inf << endl; *fLog << GetDescriptor() << " execution statistics:" << endl; *fLog << dec << setfill(' '); *fLog << " " << setw(7) << fOverFlow << " (" << setw(3) << (int)(fOverFlow*100/GetNumExecutions()) << "%) Evts skipped due to: lo gain saturated." << endl; return kTRUE; }