Index: trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc =================================================================== --- trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc (revision 5150) +++ trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc (revision 5152) @@ -80,4 +80,6 @@ // - SetBinningResolution(); // +// Sets all weights to 1. +// MExtractTimeAndChargeDigitalFilter::MExtractTimeAndChargeDigitalFilter(const char *name, const char *title) { @@ -90,8 +92,5 @@ SetBinningResolution(); - for (int i=0; i<60; i++){ - fw_amp[i] = 1; - fw_time[i] = 1; - } + SetWeightsFile(""); } @@ -161,6 +160,28 @@ // This is mandatory for the extraction // +// If filenname is empty, then all weights will be set to 1. +// void MExtractTimeAndChargeDigitalFilter::SetWeightsFile(TString filename) { + + fAmpWeightsHiGain .Set(fBinningResolution*fWindowSizeHiGain); + fAmpWeightsLoGain .Set(fBinningResolution*fWindowSizeLoGain); + fTimeWeightsHiGain.Set(fBinningResolution*fWindowSizeHiGain); + fTimeWeightsLoGain.Set(fBinningResolution*fWindowSizeLoGain); + + if (filename.IsNull()) + { + for (UInt_t i=0; iGetBinContent(i+1); - fw_time[i] = hw_time->GetBinContent(i+1); - // y[i]=hshape->GetBinContent(i+1); - // derivative[i]=hderivative->GetBinContent(i+1); - - // cout << "w_amp[" << i << "] = " << fw_amp[i] << endl; - } + for (UInt_t i=0; iGetBinContent(i+1); + fTimeWeightsHiGain[i] = hw_time->GetBinContent(i+1); + } + for (UInt_t i=0; iGetBinContent(i+1); + fTimeWeightsLoGain[i] = hw_time->GetBinContent(i+1); + } + // cout << "w_amp[" << i << "] = " << fw_amp[i] << endl; delete f; +} + +// -------------------------------------------------------------------------- +// +// ReInit +// +// Calls: +// - MExtractor::ReInit(pList); +// - Creates new arrays according to the extraction range +// +Bool_t MExtractTimeAndChargeDigitalFilter::ReInit(MParList *pList) +{ + + if (!MExtractTimeAndCharge::ReInit(pList)) + return kFALSE; + + Int_t range = fHiGainLast - fHiGainFirst + 1 + fHiLoLast; + + fHiGainSignal.Set(range); + + range = fLoGainLast - fLoGainFirst + 1; + + fLoGainSignal.Set(range); + + return kTRUE; } @@ -199,12 +247,71 @@ { - const Byte_t *end = ptr + fHiGainLast - fHiGainFirst + 1; + Int_t range = fHiGainLast - fHiGainFirst + 1; + const Byte_t *end = ptr + range; + Byte_t *p = ptr; + Byte_t maxpos = 0; + Byte_t max = 0; + Int_t count = 0; + + // + // Check for saturation in all other slices + // + while (p max) + { + max = *p; + maxpos = p-ptr; + } + + if (*p++ >= fSaturationLimit) + sat++; + } + + if (fHiLoLast != 0) + { + + end = logain + fHiLoLast; + + while (logain max) + { + max = *logain; + maxpos = range; + } + + if (*logain++ >= fSaturationLimit) + sat++; + } + } + + // + // allow one saturated slice + // + if (sat > 0) + return; + +#if 0 + if (maxpos < 2) + { + time = (Float_t)fHiGainFirst; + return; + } +#endif Float_t time_sum = 0.; Float_t fmax = 0.; Float_t ftime_max = 0.; - - Float_t pedes = ped.GetPedestal(); - + Int_t max_p = 0; + + const Float_t pedes = ped.GetPedestal(); const Float_t ABoffs = ped.GetPedestalABoffset(); @@ -216,23 +323,122 @@ // Calculate the sum of the first fWindowSize slices // - Byte_t *p = ptr; - Byte_t *max_p = ptr; - - while (p < end-fWindowSizeHiGain) + for (Int_t i=0;ifmax) + { + fmax = sum; + ftime_max = time_sum; + max_p = i; + } + } /* for (Int_t i=0;i fBinningResolutionHalf-1 || t_iter < -fBinningResolutionHalf ) + { + if (t_iter > fBinningResolutionHalf-1) + { + t_iter -= fBinningResolution; + max_p--; + sample_iter--; + } + if (t_iter < -fBinningResolutionHalf) + { + t_iter += fBinningResolution; + max_p++; + sample_iter++; + } + } + + sum = 0.; + // + // Slide with a window of size fWindowSizeHiGain over the sample + // and multiply the entries with the corresponding weights + // + for (Int_t sample=0; sample < fWindowSizeHiGain; sample++) + { + const Int_t idx = fBinningResolution*sample + fBinningResolutionHalf + t_iter; + const Int_t ids = max_p + sample; + const Float_t pex = ids < 0 ? 0. : + ( ids > range ? 0. : fHiGainSignal[ids]-PedMean[(ids+abflag) & 0x1]); + sum += fAmpWeightsHiGain [idx]*pex; + time_sum += fTimeWeightsHiGain[idx]*pex; + // cout << idx << " " << fw_amp[idx] << " " << pex << " " << t_iter << " " << sum << " " << ((Int_t(p2-ptr)+abflag) & 0x1) << " " << PedMean[(Int_t(p2-ptr)+abflag) & 0x1] << endl; + } + + if (sum != 0.) + time = max_p + 1. + sample_iter + - ((Float_t)t_iter)/fBinningResolution - 0.4 - time_sum/sum + + (Float_t)fHiGainFirst; + else + time = 0.; + } /* if (max!=0) */ + else + time=0.; + + return; +} + +void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeLoGain(Byte_t *ptr, Float_t &sum, Float_t &dsum, + Float_t &time, Float_t &dtime, + Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag) +{ + + const Byte_t *end = ptr + fLoGainLast - fLoGainFirst + 1; + + Float_t time_sum = 0; + Float_t fmax = 0; + Float_t ftime_max = 0; + + Float_t pedes = ped.GetPedestal(); + + const Float_t ABoffs = ped.GetPedestalABoffset(); + + Float_t PedMean[2]; + PedMean[0] = pedes + ABoffs; + PedMean[1] = pedes - ABoffs; + + // + // Calculate the sum of the first fWindowSize slices + // + sat = 0; + Byte_t *p = ptr; + Byte_t *max_p = ptr; + + while (p < end-fWindowSizeLoGain) + { + sum = 0; + time_sum = 0; + + // + // Slide with a window of size fWindowSizeLoGain over the sample + // and multiply the entries with the corresponding weights + // + Byte_t *p1 = p++; + for (Int_t sample=0; sample < fWindowSizeLoGain; sample++) + { + const Int_t idx = fBinningResolution*sample+fBinningResolutionHalf; + const Int_t adx = (Int_t(p1-ptr)+abflag) & 0x1; + sum += fAmpWeightsLoGain [idx]*(*p1-PedMean[adx]); + time_sum += fTimeWeightsLoGain[idx]*(*p1-PedMean[adx]); p1++; } @@ -244,41 +450,5 @@ max_p = p; } - p++; - } /* while (pfmax) - { - fmax = sum; - ftime_max = time_sum; - max_p = p; - } - - } -#endif + } /* while (pfmax) - { - fmax = sum; - ftime_max = time_sum; - max_p = p; - } - } /* while (pfBinningResolutionHalf-1 || t_iter<-fBinningResolutionHalf) - { - if (t_iter>fBinningResolutionHalf-1) - { - t_iter -= fBinningResolution; - p--; - sample_iter--; - } - if (t_iter < -fBinningResolutionHalf) - { - t_iter += fBinningResolution; - p++; - sample_iter++; - } - } - - sum = 0; - Byte_t *p2 = p; - - // - // Slide with a window of size fWindowSizeLoGain over the sample - // and multiply the entries with the corresponding weights - // - for (Int_t sample=0; sample < fWindowSizeLoGain; sample++) - { - const Int_t idx = fBinningResolution*sample+fBinningResolutionHalf+t_iter; - const Int_t adx = (Int_t(p2-ptr)+abflag) & 0x1; - sum += fw_amp [idx]*(*p2-PedMean[adx]); - time_sum += fw_time[idx]*(*p2-PedMean[adx]); - p2++; - } - - if (sum != 0) - time = (Float_t)(max_p - ptr) + 1. + sample_iter - - (Float_t)t_iter/fBinningResolution - 0.4 - time_sum/sum ; - else - time = 0.; - } /* if (max!=0) */ - else - time=0.; - - return; -} - Int_t MExtractTimeAndChargeDigitalFilter::ReadEnv(const TEnv &env, TString prefix, Bool_t print) { Index: trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.h =================================================================== --- trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.h (revision 5150) +++ trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.h (revision 5152) @@ -4,4 +4,8 @@ #ifndef MARS_MExtractTimeAndCharge #include "MExtractTimeAndCharge.h" +#endif + +#ifndef MARS_MArrayF +#include "MArrayF.h" #endif @@ -19,4 +23,7 @@ static const Int_t fgBinningResolution; + MArrayF fHiGainSignal; //! Need fast access to the signals in a float way + MArrayF fLoGainSignal; //! Store them in separate arrays + Int_t fWindowSizeHiGain; Int_t fWindowSizeLoGain; @@ -25,6 +32,10 @@ Int_t fBinningResolutionHalf; - Float_t fw_amp [60]; - Float_t fw_time[60]; + MArrayF fAmpWeightsHiGain; + MArrayF fTimeWeightsHiGain; + MArrayF fAmpWeightsLoGain; + MArrayF fTimeWeightsLoGain; + + Bool_t ReInit( MParList *pList ); Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);