/* ======================================================================== *\ ! ! * ! * 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 // // Task to calculate the calibration conversion factors from the FADC // time slices. The integrated time slices have to be delivered by an // MExtractedSignalCam. The pedestals by an MPedestalCam. // // The output container MCalibrationCam holds one entry of type MCalibrationPix // for every pixel. It is filled in the following way: // // ProProcess: Search for MPedestalCam, MExtractedSignalCam // Initialize MCalibrationCam // Initialize pulser light wavelength // // ReInit: MCalibrationCam::InitSize(NumPixels) is called which allocates // memory in a TClonesArray of type MCalibrationPix // Initialize number of used FADC slices // Optionally exclude pixels from calibration // // Process: 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 histograms of the PIN Diode are fitted to Gaussians // // Fits can be excluded via the commands: // MalibrationCam::SkipBlindPixelFits() (skip all blind // pixel fits) // MalibrationCam::SkipPinDiodeFits() (skip all PIN Diode // fits) // // Hi-Gain vs. Lo-Gain Calibration (very memory-intensive) // can be skipped with the command: // MalibrationCam::SkipHiLoGainCalibration() // // Input Containers: // MRawEvtData // MPedestalCam // MExtractedSignalCam // // Output Containers: // MCalibrationCam // ////////////////////////////////////////////////////////////////////////////// #include "MCalibrationCalc.h" // FXIME: Usage of fstream is a preliminary workaround! #include // FXIME: This has to be removed!!!! (YES, WHEN WE HAVE ACCESS TO THE DATABASE!!!!!) #include "MCalibrationConfig.h" #include #include #include "MLog.h" #include "MLogManip.h" #include "MParList.h" #include "MGeomCam.h" #include "MRawRunHeader.h" #include "MRawEvtPixelIter.h" #include "MPedestalCam.h" #include "MPedestalPix.h" #include "MCalibrationCam.h" #include "MCalibrationPix.h" #include "MExtractedSignalCam.h" #include "MExtractedSignalPix.h" #include "MCalibrationBlindPix.h" #include "MCalibrationPINDiode.h" ClassImp(MCalibrationCalc); using namespace std; const UInt_t MCalibrationCalc::fBlindPixelId = 559; const UInt_t MCalibrationCalc::fPINDiodeId = 9999; const Byte_t MCalibrationCalc::fgSaturationLimit = 254; const Byte_t MCalibrationCalc::fgBlindPixelFirst = 3; const Byte_t MCalibrationCalc::fgBlindPixelLast = 12; // -------------------------------------------------------------------------- // // Default constructor. // MCalibrationCalc::MCalibrationCalc(const char *name, const char *title) : fPedestals(NULL), fCalibrations(NULL), fSignals(NULL), fRawEvt(NULL), fRunHeader(NULL), fEvtTime(NULL) { 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"); Clear(); SetBlindPixelRange(); } void MCalibrationCalc::Clear(const Option_t *o) { SETBIT(fFlags, kUseBlindPixelFit); SETBIT(fFlags, kUseQualityChecks); SETBIT(fFlags, kHiLoGainCalibration); CLRBIT(fFlags, kHiGainOverFlow); CLRBIT(fFlags, kLoGainOverFlow); // As long as we don't have the PIN Diode: CLRBIT(fFlags, kUsePinDiodeFit); fBlindPixelFirst = 0; fBlindPixelLast = 0; fNumHiGainSamples = 0; fNumLoGainSamples = 0; fConversionHiLo = 0; fNumExcludedPixels = 0; fColor = kECT1; } void MCalibrationCalc::SetBlindPixelRange(Byte_t first, Byte_t last) { fBlindPixelFirst = first; fBlindPixelLast = last; } MCalibrationBlindPix *MCalibrationCalc::GetBlindPixel() const { return fCalibrations->GetBlindPixel(); } MCalibrationPINDiode *MCalibrationCalc::GetPINDiode() const { return fCalibrations->GetPINDiode(); } // -------------------------------------------------------------------------- // // 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 // // The following output containers are only searched, but not created // // - MTime // Int_t MCalibrationCalc::PreProcess(MParList *pList) { fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData"); if (!fRawEvt) { *fLog << err << 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; } fEvtTime = (MTime*)pList->FindObject("MTime"); 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 << err << 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; } fCalibrations->SetGeomCam(cam); fNumHiGainSamples = fSignals->GetNumUsedHiGainFADCSlices(); fNumLoGainSamples = fSignals->GetNumUsedLoGainFADCSlices(); fSqrtHiGainSamples = TMath::Sqrt((Float_t)fNumHiGainSamples); UInt_t npixels = cam->GetNumPixels(); for (UInt_t i=0; iGetFirstUsedSliceHiGain(), fSignals->GetLastUsedSliceHiGain()); pix.SetAbsTimeBordersLoGain(fSignals->GetFirstUsedSliceLoGain(), fSignals->GetLastUsedSliceLoGain()); if (!TESTBIT(fFlags,kUseQualityChecks)) pix.SetExcludeQualityCheck(); // Exclude the blind pixel and the PIN Diode from normal pixel calibration: if (i == fBlindPixelId) pix.SetExcluded(); if (i == fPINDiodeId) pix.SetExcluded(); } // // Look for file to exclude pixels from analysis // if (!fExcludedPixelsFile.IsNull()) { fExcludedPixelsFile = gSystem->ExpandPathName(fExcludedPixelsFile.Data()); // // Initialize reading the file // ifstream in(fExcludedPixelsFile.Data(),ios::in); if (in) { *fLog << inf << "Use excluded pixels from file: '" << fExcludedPixelsFile.Data() << "'" << endl; // // Read the file and count the number of entries // UInt_t pixel = 0; while (++fNumExcludedPixels) { in >> pixel; if (!in.good()) break; // // Check for out of range // if (pixel > npixels) { *fLog << warn << "WARNING: To be excluded pixel: " << pixel << " is out of range " << endl; continue; } // // Exclude pixel // MCalibrationPix &pix = (*fCalibrations)[pixel]; pix.SetExcluded(); *fLog << GetDescriptor() << inf << ": Exclude Pixel: " << pixel << endl; } if (--fNumExcludedPixels == 0) *fLog << warn << "WARNING: File '" << fExcludedPixelsFile.Data() << "'" << " is empty " << endl; else fCalibrations->SetNumPixelsExcluded(fNumExcludedPixels); } else *fLog << warn << dbginf << "Cannot open file '" << fExcludedPixelsFile.Data() << "'" << endl; } 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() { // // Initialize pointers to blind pixel, PIN Diode and individual pixels // MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel()); MCalibrationPINDiode &pindiode = *(fCalibrations->GetPINDiode()); MRawEvtPixelIter pixel(fRawEvt); // // Create a (second) loop to do fill the calibration histograms // Search for: a signal in MExtractedSignalCam // Fill histograms with: // charge // charge vs. event nr. // relative arrival time w.r.t. pixel 1 // while (pixel.Next()) { const UInt_t pixid = pixel.GetPixelId(); MCalibrationPix &pix = (*fCalibrations)[pixid]; MExtractedSignalPix &sig = (*fSignals)[pixid]; const Float_t sumhi = sig.GetExtractedSignalHiGain(); const Float_t sumlo = sig.GetExtractedSignalLoGain(); Float_t abstime = 0.; #if 0 Float_t reltime = 0.; if (TESTBIT(fFlags,kUseTimes)) { // // Have a look in MArrivalTime, // otherwise search the position of maximum bin // in MRawEvtData // if (fArrivalTime) { abstime = (*fArrivalTime)[pixid]; reltime = abstime - (*fArrivalTime)[1]; } else { if (pixid == 1) referencetime = (Float_t)pixel.GetIdxMaxHiGainSample(); if (sig.IsLoGainUsed()) { abstime = (Float_t)pixel.GetIdxMaxLoGainSample(); // reltime = abstime - referencetime; } else { abstime = (Float_t)pixel.GetIdxMaxHiGainSample(); // reltime = abstime - referencetime; } // } // } /* if Use Times */ #endif switch(pixid) { case fBlindPixelId: if (TESTBIT(fFlags,kUseBlindPixelFit)) { Byte_t *ptr = pixel.GetHiGainSamples(); Float_t blindpixelsumhi = 0.; Float_t blindpixelsumlo = 0.; // // We need a dedicated signal extractor for the blind pixel // Int_t diff = 0; if (fBlindPixelLast > 15) { diff = fBlindPixelLast - 15; fBlindPixelLast = 15; } Byte_t *start = ptr + fBlindPixelFirst - 1; Byte_t *end = start + fBlindPixelLast - fBlindPixelFirst + 1; ptr = start; Int_t sum = 0; while (ptr 0) { ptr = pixel.GetLoGainSamples(); end = ptr + diff + 1; while (ptrCutEdges(); // // Fit the blind pixel // if (TESTBIT(fFlags,kUseBlindPixelFit)) { // // Get pointer to blind pixel // MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel()); *fLog << inf << GetDescriptor() << ": Fitting the Blind Pixel" << endl; // // retrieve mean and sigma of the blind pixel pedestal, // so that we can use it for the fit // if (fPedestals->GetHistSize() > fBlindPixelId) { Float_t pedestal; Float_t pederr; Float_t pedsigma; Float_t pedsigmaerr; const ULong_t nentries = fPedestals->GetTotalEntries(); const Float_t nslices = (Float_t)(fgBlindPixelLast-fgBlindPixelFirst+1); const Float_t sqrslice = TMath::Sqrt(nslices); // // retrieve the pedestal pix of the blind pixel // if (fPedestals->GetHistSize() != 0) { MHPedestalPixel &pedhist = (*fPedestals)(fBlindPixelId); pedestal = pedhist.GetChargeMean()*nslices; pederr = pedhist.GetChargeMeanErr()*nslices; // // Fitted sigma: 1. one sqrt(Nr. slices) for the division which is not // not appropriate: sigma(real)/slice = GetSigma*sqrt(nslices) // 2. another sqrt(Nr. slices) to calculate back to number // of slices // pedsigma = pedhist.GetChargeSigma()*nslices; pedsigmaerr = pedhist.GetChargeSigmaErr()*nslices; } else { MPedestalPix &pedpix = (*fPedestals)[fBlindPixelId]; pedestal = pedpix.GetPedestal()*nslices; pederr = pedpix.GetPedestalRms()*nslices/nentries; pedsigma = pedpix.GetPedestalRms()*sqrslice; pedsigmaerr = pederr/2.; } // // retrieve the histogram containers // MHCalibrationBlindPixel *hist = blindpixel.GetHist(); hist->SetMeanPedestal(pedestal); hist->SetMeanPedestalErr(pederr); hist->SetSigmaPedestal(pedsigma); hist->SetSigmaPedestalErr(pedsigmaerr); } if (!blindpixel.FitCharge()) { *fLog << warn << "Could not fit the blind pixel! " << endl; *fLog << warn << "Setting bit kBlindPixelMethodValid to FALSE in MCalibrationCam" << endl; fCalibrations->SetBlindPixelMethodValid(kFALSE); } else fCalibrations->SetBlindPixelMethodValid(kTRUE); if (blindpixel.CheckOscillations()) fCalibrations->SetBlindPixelMethodValid(kFALSE); blindpixel.DrawClone(); } else *fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Fit " << endl; *fLog << inf << 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]; // // Check if the pixel has been excluded from the fits // if (pix.IsExcluded()) continue; // // get the pedestals // const Float_t ped = (*fPedestals)[pixid].GetPedestal(); const Float_t prms = (*fPedestals)[pixid].GetPedestalRms(); // // set them in the calibration camera // pix.SetPedestal(ped,prms,(Float_t)fNumHiGainSamples,(Float_t)fNumLoGainSamples); // // perform the Gauss fits to the charges // pix.FitCharge(); // // check also for oscillations // pix.CheckOscillations(); } if (TESTBIT(fFlags,kUseBlindPixelFit) && fCalibrations->IsBlindPixelMethodValid()) { if (!fCalibrations->CalcFluxInsidePlexiglass()) { *fLog << err << "Could not calculate the number of photons from the blind pixel " << endl; *fLog << err << "You can try to calibrate using the MCalibrationCalc::SkipBlindPixelFit()" << endl; fCalibrations->SetBlindPixelMethodValid(kFALSE); } } else *fLog << inf << GetDescriptor() << ": Skipping Blind Pixel Calibration! " << endl; if (TESTBIT(fFlags,kUsePinDiodeFit) && fCalibrations->IsPINDiodeMethodValid()) { if (!fCalibrations->CalcFluxOutsidePlexiglass()) { *fLog << err << "Could not calculate the number of photons from the blind pixel " << endl; *fLog << err << "You can try to calibrate using the MCalibrationCalc::SkipPINDiodeFit()" << endl; fCalibrations->SetPINDiodeMethodValid(kFALSE); } } else *fLog << inf << GetDescriptor() << ": Skipping PIN Diode Calibration! " << endl; fCalibrations->SetReadyToSave(); return kTRUE; } Bool_t MCalibrationCalc::CalcSignalBlindPixel(Byte_t *ptr, Float_t &signal) const { Byte_t *newptr = ptr; Byte_t *start = ptr + fBlindPixelFirst-1; Byte_t *end = ptr + fBlindPixelLast; Byte_t sum = 0; Int_t sat = 0; newptr = start; while (ptr= fgSaturationLimit) sat++; } if (sat) { *fLog << err << "HI Gain Saturation occurred in the blind pixel! " << " Do not know yet how to treat this ... aborting " << endl; *fLog << err << "If you need absolutely any other kind of calibration, " << " use SkipBlindPixelFit() " << endl; return kFALSE; } signal = (Float_t)sum; return kTRUE; }