Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 5853) +++ trunk/MagicSoft/Mars/Changelog (revision 5854) @@ -31,4 +31,6 @@ ( SetSignalType(MCalibrateData::kPhe) ). Default is old version in photons + - speed up Process by storing pre-calculated calibration constants + in arrays (needed 40% of CPU time of the eventloop before, now: ) Index: trunk/MagicSoft/Mars/mcalib/MCalibrateData.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrateData.cc (revision 5853) +++ trunk/MagicSoft/Mars/mcalib/MCalibrateData.cc (revision 5854) @@ -280,4 +280,7 @@ } + fCalibConsts.Reset(); + fCalibFFactors.Reset(); + return kTRUE; } @@ -378,4 +381,39 @@ } + const Int_t npixels = fGeomCam->GetNumPixels(); + + if (fCalibrationMode > kNone) + { + + if (fCalibrations->GetSize() != npixels) + { + *fLog << err << GetDescriptor() + << ": Size mismatch between MGeomCam and MCalibrationChargeCam ... abort!" << endl; + return kFALSE; + } + + if (fBadPixels->GetSize() != npixels) + { + *fLog << err << GetDescriptor() + << ": Size mismatch between MGeomCam and MBadPixelsCam ... abort!" << endl; + return kFALSE; + } + + if (fBadPixels->GetSize() != npixels) + { + *fLog << err << GetDescriptor() + << ": Size mismatch between MGeomCam and MBadPixelsCam ... abort!" << endl; + return kFALSE; + } + } + + fCalibConsts .Set(npixels); + fCalibFFactors.Set(npixels); + fHiLoConv .Set(npixels); + fHiLoConvErr .Set(npixels); + + if (!UpdateConversionFactors()) + return kFALSE; + if (TestPedestalFlag(kRun)) Calibrate(kFALSE, kTRUE); @@ -386,9 +424,12 @@ // -------------------------------------------------------------------------- // -// Get conversion factor and its error from MCalibrationCam -// -Bool_t MCalibrateData::GetConversionFactor(UInt_t pixidx, Float_t &hiloconv, Float_t &hiloconverr, - Float_t &calibConv, Float_t &calibConvVar, - Float_t &calibFFactor) const +// Update the conversion factors and F-Factors from MCalibrationCams into +// the arrays. Later, the here pre-calcualted conversion factors get simply +// copied from memory. +// +// This function can be called from outside in case that the MCalibrationCams +// have been updated... +// +Bool_t MCalibrateData::UpdateConversionFactors() { // @@ -397,93 +438,131 @@ const Float_t zenith = -1.; - hiloconv = 1.; - hiloconverr = 0.; - calibConv = 1.; - calibConvVar = 0.; - calibFFactor = 0.; - - Float_t calibQE = 1.; - Float_t calibQEVar = 0.; - - if(fCalibrationMode!=kNone) - { - MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCalibrations)[pixidx]; - - hiloconv = pix.GetConversionHiLo (); - hiloconverr= pix.GetConversionHiLoErr(); - - if ((*fBadPixels)[pixidx].IsUnsuitable()) - return kFALSE; - - calibConv = pix.GetMeanConvFADC2Phe(); - calibConvVar = pix.GetMeanConvFADC2PheVar(); - calibFFactor = pix.GetMeanFFactorFADC2Phot(); - - MCalibrationQEPix &qe = (MCalibrationQEPix&) (*fQEs)[pixidx]; - - switch(fCalibrationMode) + UInt_t skip = 0; + + for (UInt_t pixidx=0; pixidxGetNumPixels(); pixidx++) + { + + Float_t hiloconv = 1.; + Float_t hiloconverr = 0.; + Float_t calibConv = 1.; + Float_t calibConvVar = 0.; + Float_t calibFFactor = 0.; + + Float_t calibQE = 1.; + Float_t calibQEVar = 0.; + + if(fCalibrationMode!=kNone) { - case kFlatCharge: - { - MCalibrationChargePix &avpix = (MCalibrationChargePix&)fCalibrations->GetAverageArea(0); - calibConv = avpix.GetMean() / (pix.GetMean() * fGeomCam->GetPixRatio(pixidx)); - calibConvVar = (avpix.GetMeanRelVar() + pix.GetMeanRelVar()) * calibConv * calibConv; - if (pix.IsFFactorMethodValid()) + MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCalibrations)[pixidx]; + + hiloconv = pix.GetConversionHiLo (); + hiloconverr= pix.GetConversionHiLoErr(); + + if ((*fBadPixels)[pixidx].IsUnsuitable()) + { + skip++; + continue; + } + + calibConv = pix.GetMeanConvFADC2Phe(); + calibConvVar = pix.GetMeanConvFADC2PheVar(); + calibFFactor = pix.GetMeanFFactorFADC2Phot(); + + MCalibrationQEPix &qe = (MCalibrationQEPix&) (*fQEs)[pixidx]; + + switch(fCalibrationMode) + { + case kFlatCharge: { - const Float_t convmin1 = qe.GetQECascadesFFactor(zenith)/pix.GetMeanConvFADC2Phe(); - if (convmin1 > 0) - calibFFactor *= TMath::Sqrt(convmin1); - else - calibFFactor = -1.; + MCalibrationChargePix &avpix = (MCalibrationChargePix&)fCalibrations->GetAverageArea(0); + calibConv = avpix.GetMean() / (pix.GetMean() * fGeomCam->GetPixRatio(pixidx)); + calibConvVar = (avpix.GetMeanRelVar() + pix.GetMeanRelVar()) * calibConv * calibConv; + if (pix.IsFFactorMethodValid()) + { + const Float_t convmin1 = qe.GetQECascadesFFactor(zenith)/pix.GetMeanConvFADC2Phe(); + if (convmin1 > 0) + calibFFactor *= TMath::Sqrt(convmin1); + else + calibFFactor = -1.; + } + break; } - break; - } - case kBlindPixel: - if (!qe.IsBlindPixelMethodValid()) - return kFALSE; - calibQE = qe.GetQECascadesBlindPixel ( zenith ); - calibQEVar = qe.GetQECascadesBlindPixelVar( zenith ); - break; - - case kPinDiode: - if (!qe.IsPINDiodeMethodValid()) - return kFALSE; - calibQE = qe.GetQECascadesPINDiode ( zenith ); - calibQEVar = qe.GetQECascadesPINDiodeVar( zenith ); - break; - - case kFfactor: - if (!pix.IsFFactorMethodValid()) - return kFALSE; - calibQE = qe.GetQECascadesFFactor ( zenith ); - calibQEVar = qe.GetQECascadesFFactorVar( zenith ); - break; - - case kCombined: - if (!qe.IsCombinedMethodValid()) - return kFALSE; - calibQE = qe.GetQECascadesCombined ( zenith ); - calibQEVar = qe.GetQECascadesCombinedVar( zenith ); - break; - - case kDummy: - hiloconv = 1.; - hiloconverr = 0.; - break; - } /* switch calibration mode */ - } /* if(fCalibrationMode!=kNone) */ - else - { - calibConv = 1./fGeomCam->GetPixRatio(pixidx); - } - - calibConv /= calibQE; - - if (calibConv != 0. && calibQE != 0.) - { - // Now doing: - // calibConvVar = calibConvVar/(calibConv*calibConv) + calibQEVar/(calibQE*calibQE); - // calibConvVar *= (calibConv*calibConv); - calibConvVar += calibQEVar*(calibConv*calibConv)/(calibQE*calibQE); + case kBlindPixel: + if (!qe.IsBlindPixelMethodValid()) + { + skip++; + continue; + } + calibQE = qe.GetQECascadesBlindPixel ( zenith ); + calibQEVar = qe.GetQECascadesBlindPixelVar( zenith ); + break; + + case kPinDiode: + if (!qe.IsPINDiodeMethodValid()) + { + skip++; + continue; + } + calibQE = qe.GetQECascadesPINDiode ( zenith ); + calibQEVar = qe.GetQECascadesPINDiodeVar( zenith ); + break; + + case kFfactor: + if (!pix.IsFFactorMethodValid()) + { + skip++; + continue; + } + calibQE = qe.GetQECascadesFFactor ( zenith ); + calibQEVar = qe.GetQECascadesFFactorVar( zenith ); + break; + + case kCombined: + if (!qe.IsCombinedMethodValid()) + { + skip++; + continue; + } + calibQE = qe.GetQECascadesCombined ( zenith ); + calibQEVar = qe.GetQECascadesCombinedVar( zenith ); + break; + + case kDummy: + hiloconv = 1.; + hiloconverr = 0.; + break; + } /* switch calibration mode */ + } /* if(fCalibrationMode!=kNone) */ + else + { + calibConv = 1./fGeomCam->GetPixRatio(pixidx); + } + + calibConv /= calibQE; + + if (calibConv != 0. && calibQE != 0.) + { + // Now doing: + calibConvVar = calibConvVar/(calibConv*calibConv) + calibQEVar/(calibQE*calibQE); + calibConvVar *= (calibConv*calibConv); + // The above two lines had been commented by TB and replaced by the following line + // (without notice to me!) but it is simply wrong! + // calibConvVar += calibQEVar*(calibConv*calibConv)/(calibQE*calibQE); + } + + calibConv *= fRenormFactor; + + fHiLoConv [pixidx] = hiloconv; + fHiLoConvErr [pixidx] = hiloconverr; + fCalibConsts [pixidx] = calibConv; + fCalibFFactors[pixidx] = calibFFactor; + + } /* for (Int_t pixidx=0; pixidxGetNumPixels(); pixidx++) */ + + if (skip>fGeomCam->GetNumPixels()*0.9) + { + *fLog << warn << GetDescriptor() + << ": WARNING - GetConversionFactor has skipped more than 90% of the pixels... abort." << endl; + return kFALSE; } @@ -491,4 +570,12 @@ } + +// -------------------------------------------------------------------------- +// +// Apply the conversion factors and F-Factors from the arrays to the data. +// +// The flags 'data' and 'pedestal' decide whether the signal and/or the pedetals +// shall be calibrated, respectively. +// Int_t MCalibrateData::Calibrate(Bool_t data, Bool_t pedestal) const { @@ -500,23 +587,8 @@ const Float_t sqrtslices = TMath::Sqrt(slices); - Float_t hiloconv; - Float_t hiloconverr; - Float_t calibConv; - Float_t calibConvErr; - Float_t calibFFactor; - - UInt_t skip = 0; for (UInt_t pixidx=0; pixidx 0 ? signalErr*signalErr / (signal * signal) : 0. - + calibConv > 0 ? calibConvVar / (calibConv * calibConv ) : 0.; - nphotErr = TMath::Sqrt(nphotErr) * nphot; - */ + const Float_t nphot = signal * fCalibConsts [pixidx]; + const Float_t nphotErr = TMath::Sqrt(TMath::Abs(nphot)) * fCalibFFactors[pixidx]; MCerPhotPix *cpix = fCerPhotEvt->AddPixel(pixidx, nphot, nphotErr); @@ -555,22 +617,10 @@ if (sig.GetNumLoGainSaturated() > 0) - cpix->SetPixelSaturated(); - } + cpix->SetPixelSaturated(); + } /* if (data) */ + if (pedestal) { - /* - // pedestals/(used FADC slices) in [ADC] counts - const Float_t pedes = (*fPedestalMean)[pixidx].GetPedestal() * slices; - const Float_t pedrms = (*fPedestalRms)[pixidx].GetPedestalRms() * sqrtslices; - - // - // pedestals/(used FADC slices) in [number of photons] - // - const Float_t pedphot = pedes * calibConv; - const Float_t pedphotrms = pedrms * calibConv; - - (*fPedPhot)[pixidx].Set(pedphot, pedphotrms); - */ TIter NextPed(&fPedestalCams); TIter NextPhot(&fPedPhotCams); @@ -578,26 +628,19 @@ MPedestalCam *pedestal = 0; MPedPhotCam *pedphot = 0; + + const Float_t pedmeancalib = slices *fCalibConsts[pixidx]; + const Float_t pedrmscalib = sqrtslices*fCalibConsts[pixidx]; while ((pedestal=(MPedestalCam*)NextPed()) && (pedphot =(MPedPhotCam*)NextPhot())) { - // pedestals/(used FADC slices) in [ADC] counts - const Float_t pedes = (*pedestal)[pixidx].GetPedestal() * slices; - const Float_t pedrms = (*pedestal)[pixidx].GetPedestalRms() * sqrtslices; - - // pedestals/(used FADC slices) in [number of photons] - const Float_t mean = pedes * calibConv; - const Float_t rms = pedrms * calibConv; - - (*pedphot)[pixidx].Set(mean, rms); - pedphot->SetReadyToSave(); + // pedestals/(used FADC slices) in [number of photons] + const Float_t mean = (*pedestal)[pixidx].GetPedestal() *pedmeancalib; + const Float_t rms = (*pedestal)[pixidx].GetPedestalRms()*pedrmscalib; + + (*pedphot)[pixidx].Set(mean, rms); + pedphot->SetReadyToSave(); } - } - } - - if (skip>npix*0.9) - { - *fLog << warn << "WARNING - GetConversionFactor has skipped more than 90% of the pixels... skip." << endl; - return kCONTINUE; + } /* if (pedestal) */ } @@ -617,16 +660,4 @@ Int_t MCalibrateData::Process() { - /* - if (fCalibrations->GetNumPixels() != (UInt_t)fSignals->GetSize()) - { - // FIXME: MExtractedSignal must be of variable size - - // like MCerPhotEvt - because we must be able - // to reduce size by zero supression - // For the moment this check could be done in ReInit... - *fLog << err << "MExtractedSignal and MCalibrationCam have different sizes... abort." << endl; - return kFALSE; - } - */ - return Calibrate(fCalibrationMode!=kSkip, TestPedestalFlag(kEvent)); } Index: trunk/MagicSoft/Mars/mcalib/MCalibrateData.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrateData.h (revision 5853) +++ trunk/MagicSoft/Mars/mcalib/MCalibrateData.h (revision 5854) @@ -21,4 +21,8 @@ #endif +#ifndef MARS_MArrayF +#include "MArrayF.h" +#endif + class MGeomCam; class MBadPixelsCam; @@ -34,4 +38,5 @@ { private: + MGeomCam *fGeomCam; //! Camera geometry container MBadPixelsCam *fBadPixels; //! Bad Pixels information @@ -51,7 +56,10 @@ TList fPedPhotCams; //! List of pointers to corresponding output MPedPhotCam - Int_t Calibrate(Bool_t data, Bool_t pedestal) const; + MArrayF fCalibConsts; //! Array of calibration constants for each pixel, calculated only once! + MArrayF fCalibFFactors; //! Array of calibration F-Factors for each pixel, calculated only once! + MArrayF fHiLoConv; //! Array of calibration constants for each pixel, calculated only once! + MArrayF fHiLoConvErr; //! Array of calibration F-Factors for each pixel, calculated only once! - Bool_t GetConversionFactor(UInt_t, Float_t &, Float_t &, Float_t &, Float_t &, Float_t &) const; + Int_t Calibrate(Bool_t data, Bool_t pedestal) const; Int_t PreProcess(MParList *pList); @@ -93,4 +101,7 @@ const char *name=NULL, const char *title=NULL); + void AddPedestal(const char *name="Cam"); + void AddPedestal(const char *pedestal, const char *pedphot); + void EnablePedestalType(PedestalType_t i) { fPedestalFlag |= i; } void SetPedestalFlag(PedestalType_t i=kRun) { fPedestalFlag = i; } @@ -100,7 +111,6 @@ void SetSignalType ( SignalType_t sigtype=kPhot ) { fSignalType =sigtype; } - void AddPedestal(const char *name="Cam"); - void AddPedestal(const char *pedestal, const char *pedphot); - + Bool_t UpdateConversionFactors(); + ClassDef(MCalibrateData, 1) // Task to calibrate FADC counts into Cherenkov photons };