/* ======================================================================== *\ ! ! * ! * 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: // // HCharge(npix) (distribution of summed FADC time slice entries) // HTime(npix) (distribution of position of maximum) // HChargevsN(npix) (distribution of charges vs. event number. // // PostProcess: All histograms HCharge(npix) are fitted to a Gaussian // All histograms HTime(npix) are fitted to a Gaussian // The histogram HBlindPixelCharge (blind pixel) is fitted to a single // PhE fit // The histogram HBlindPixelTime (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::SetSkipTimeFits() (skip all time fits) // MalibrationCam::SetSkipBlindPixelFits() (skip all blind pixel fits) // MalibrationCam::SetSkipPinDiodeFits() (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 "MGeomCam.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 "MExtractedSignalCam.h" #include "MExtractedSignalPix.h" #include "MTime.h" #include "TMath.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, kUseTimeFits); SETBIT(fFlags, kUseBlindPixelFit); SETBIT(fFlags, kUsePinDiodeFit); } // -------------------------------------------------------------------------- // // 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) { fHistOverFlow = 0; fEvents = 0; fCosmics = 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; } switch (fColor) { case kEBlue: fCalibrations->SetColor(MCalibrationCam::kECBlue); break; case kEGreen: fCalibrations->SetColor(MCalibrationCam::kECGreen); break; case kEUV: fCalibrations->SetColor(MCalibrationCam::kECUV); break; case kECT1: fCalibrations->SetColor(MCalibrationCam::kECCT1); break; default: fCalibrations->SetColor(MCalibrationCam::kECCT1); } fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam"); if (!fPedestals) { *fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl; return kFALSE; } fSignals = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam"); if (!fSignals) { *fLog << err << dbginf << "Cannot find MExtractedSignalCam ... aborting" << endl; return kFALSE; } 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; } MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam"); if (!cam) { *fLog << err << GetDescriptor() << ": No MGeomCam found... aborting." << endl; return kFALSE; } fNumHiGainSamples = fSignals->GetNumUsedHiGainFADCSlices(); fNumLoGainSamples = fSignals->GetNumUsedLoGainFADCSlices(); fCalibrations->InitSize(cam->GetNumPixels()); for (UInt_t i=0;iGetNumPixels();i++) { MCalibrationPix &pix = (*fCalibrations)[i]; pix.DefinePixId(i); } 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() { Int_t cosmicpix = 0; MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel()); MCalibrationPINDiode &pindiode = *(fCalibrations->GetPINDiode()); MRawEvtPixelIter pixel(fRawEvt); // // Create a first loop to sort out the cosmics ... // // This is a very primitive check for the number of cosmicpixs // The cut will be applied in the fit, but for the blind pixel, // we need to remove this event // // FIXME: In the future need a much more sophisticated one!!! // while (pixel.Next()) { const UInt_t pixid = pixel.GetPixelId(); MExtractedSignalPix &sig = (*fSignals)[pixid]; MPedestalPix &ped = (*fPedestals)[pixid]; Float_t pedrms = ped.GetPedestalRms(); Float_t sumhi = sig.GetExtractedSignalHiGain(); if (sumhi < 15.*pedrms ) // cut at 3.5 sigma cosmicpix++; } if (cosmicpix > 100.) { fCosmics++; return kCONTINUE; } pixel.Reset(); fEvents++; // // Create a second loop to do fill the calibration histograms // while (pixel.Next()) { const UInt_t pixid = pixel.GetPixelId(); MExtractedSignalPix &sig = (*fSignals)[pixid]; Float_t sumhi = sig.GetExtractedSignalHiGain(); Float_t sumlo = sig.GetExtractedSignalLoGain(); Float_t mtime = sig.GetMeanArrivalTime(); MCalibrationPix &pix = (*fCalibrations)[pixid]; switch(pixid) { case gkCalibrationBlindPixelId: if (!blindpixel.FillCharge(sumhi)) *fLog << err << "Overflow or Underflow occurred filling Blind Pixel sum = " << sumhi << endl; if (!blindpixel.FillTime((int)mtime)) *fLog << err << "Overflow or Underflow occurred filling Blind Pixel time = " << mtime << endl; if (!blindpixel.FillRChargevsTime(sumhi,fEvents)) *fLog << warn << "Overflow or Underflow occurred filling Blind Pixel eventnr = " << fEvents << endl; case gkCalibrationPINDiodeId: if (!pindiode.FillCharge(sumhi)) *fLog << warn << "Overflow or Underflow occurred filling HCharge: means = " << sumhi << endl; if (!pindiode.FillTime((int)mtime)) *fLog << warn << "Overflow or Underflow occurred filling HTime: time = " << mtime << endl; if (!pindiode.FillRChargevsTime(sumhi,fEvents)) *fLog << warn << "Overflow or Underflow occurred filling HChargevsN: eventnr = " << fEvents << endl; break; default: pix.SetChargesInGraph(sumhi,sumlo); if (!pix.FillRChargevsTimeLoGain(sumlo,fEvents)) *fLog << warn << "Could not fill Lo Gain Charge vs. EvtNr of pixel: " << pixid << " signal = " << sumlo << " event Nr: " << fEvents << endl; if (!pix.FillRChargevsTimeHiGain(sumhi,fEvents)) *fLog << warn << "Could not fill Hi Gain Charge vs. EvtNr of pixel: " << pixid << " signal = " << sumhi << " event Nr: " << fEvents << endl; if (sig.IsLoGainUsed()) { if (!pix.FillChargeLoGain(sumlo)) *fLog << warn << "Could not fill Lo Gain Charge of pixel: " << pixid << " signal = " << sumlo << endl; if (!pix.FillTimeLoGain((int)mtime)) *fLog << warn << "Could not fill Lo Gain Time of pixel: " << pixid << " time = " << mtime << endl; } else { if (!pix.FillChargeHiGain(sumhi)) *fLog << warn << "Could not fill Hi Gain Charge of pixel: " << pixid << " signal = " << sumhi << endl; if (!pix.FillTimeHiGain((int)mtime)) *fLog << warn << "Could not fill Hi Gain Time of pixel: " << pixid << " time = " << mtime << endl; } break; } /* switch(pixid) */ } /* while (pixel.Next()) */ 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,kUseBlindPixelFit)) { if (!blindpixel.FitCharge()) *fLog << err << dbginf << "Could not fit 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++) { MCalibrationPix &pix = (*fCalibrations)[pixid]; const Float_t ped = (*fPedestals)[pixid].GetPedestal() * fNumHiGainSamples; const Float_t prms = (*fPedestals)[pixid].GetPedestalRms() * TMath::Sqrt((float)fNumHiGainSamples); pix.SetPedestal(ped,prms); pix.FitCharge(); if (TESTBIT(fFlags,kUseTimeFits)) pix.FitTime(); } if (!fCalibrations->CalcNumPhotInsidePlexiglass()) *fLog << err << dbginf << "Could not calculate the number of photons from the blind pixel " << endl; fCalibrations->SetReadyToSave(); if (GetNumExecutions()==0) return kTRUE; *fLog << endl; *fLog << dec << setfill(' ') << fCosmics << " Events presumably cosmics" << endl; return kTRUE; }