Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc (revision 5154) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc (revision 5155) @@ -60,15 +60,18 @@ #include "TString.h" +#include + ClassImp(MExtractTimeAndChargeDigitalFilter); using namespace std; -const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainFirst = 0; -const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainLast = 14; -const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainFirst = 3; -const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainLast = 14; -const Int_t MExtractTimeAndChargeDigitalFilter::fgHiGainWindowSize = 6; -const Int_t MExtractTimeAndChargeDigitalFilter::fgLoGainWindowSize = 6; -const Int_t MExtractTimeAndChargeDigitalFilter::fgBinningResolution = 10; +const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainFirst = 0; +const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainLast = 14; +const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainFirst = 3; +const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainLast = 14; +const Int_t MExtractTimeAndChargeDigitalFilter::fgWindowSizeHiGain = 6; +const Int_t MExtractTimeAndChargeDigitalFilter::fgWindowSizeLoGain = 6; +const Int_t MExtractTimeAndChargeDigitalFilter::fgBinningResolutionHiGain = 10; +const Int_t MExtractTimeAndChargeDigitalFilter::fgBinningResolutionLoGain = 10; // -------------------------------------------------------------------------- // @@ -92,5 +95,5 @@ SetBinningResolution(); - SetWeightsFile(""); + ReadWeightsFile(""); } @@ -108,4 +111,11 @@ { + if (windowh != fgWindowSizeHiGain) + *fLog << warn << GetDescriptor() + << ": ATTENTION!!! If you are not Hendrik Bartko, do NOT use a different window size than the default." << endl; + if (windowl != fgWindowSizeLoGain) + *fLog << warn << GetDescriptor() + << ": ATTENTION!!! If you are not Hendrik Bartko, do NOT use a different window size than the default" << endl; + fWindowSizeHiGain = windowh; fWindowSizeLoGain = windowl; @@ -162,11 +172,11 @@ // 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); +Bool_t MExtractTimeAndChargeDigitalFilter::ReadWeightsFile(TString filename) +{ + + fAmpWeightsHiGain .Set(fBinningResolutionHiGain*fWindowSizeHiGain); + fAmpWeightsLoGain .Set(fBinningResolutionLoGain*fWindowSizeLoGain); + fTimeWeightsHiGain.Set(fBinningResolutionHiGain*fWindowSizeHiGain); + fTimeWeightsLoGain.Set(fBinningResolutionLoGain*fWindowSizeLoGain); if (filename.IsNull()) @@ -182,37 +192,120 @@ fTimeWeightsLoGain[i] = 1.; } - return; - } - - TFile *f = new TFile(filename); - TH1F *hw_amp = (TH1F*)f->Get("hw_amp"); - TH1F *hw_time = (TH1F*)f->Get("hw_time"); - - if (!hw_amp) - { - *fLog << warn << GetDescriptor() - << ": Could not find hw_amp in " << filename << endl; - return; - } - - if (!hw_time) - { - *fLog << warn << GetDescriptor() - << ": Could not find hw_time in " << filename << endl; - return; - } - - 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; + return kTRUE; + } + + ifstream fin(filename.Data()); + + if (!fin) + { + *fLog << err << GetDescriptor() + << ": No weights file found: " << filename << endl; + return kFALSE; + } + + Int_t len = 0; + Int_t cnt = 0; + Bool_t hi = kFALSE; + Bool_t lo = kFALSE; + + TString str; + + while (1) + { + + str.ReadLine(fin); + if (!fin) + break; + + + if (str.Contains("# High Gain Weights:")) + { + str.ReplaceAll("# High Gain Weights:",""); + sscanf(str.Data(),"%2i%2i",&fWindowSizeHiGain,&fBinningResolutionHiGain); + *fLog << inf << "Found number of High Gain slices: " << fWindowSizeHiGain + << " and High Gain resolution: " << fBinningResolutionHiGain << endl; + len = fBinningResolutionHiGain*fWindowSizeHiGain; + fAmpWeightsHiGain .Set(len); + fTimeWeightsHiGain.Set(len); + hi = kTRUE; + continue; + } + + if (str.Contains("# Low Gain Weights:")) + { + str.ReplaceAll("# Low Gain Weights:",""); + sscanf(str.Data(),"%2i%2i",&fWindowSizeLoGain,&fBinningResolutionLoGain); + *fLog << inf << "Found number of Low Gain slices: " << fWindowSizeLoGain + << " and Low Gain resolution: " << fBinningResolutionLoGain << endl; + len = fBinningResolutionLoGain*fWindowSizeHiGain; + fAmpWeightsLoGain .Set(len); + fTimeWeightsLoGain.Set(len); + lo = kTRUE; + continue; + } + + if (str.Contains("#")) + continue; + + if (len == 0) + continue; + + sscanf(str.Data(),"\t%f\t%f",lo ? &fAmpWeightsLoGain [cnt] : &fAmpWeightsHiGain [cnt], + lo ? &fTimeWeightsLoGain[cnt] : &fTimeWeightsHiGain[cnt]); + + if (++cnt == len) + { + len = 0; + cnt = 0; + } + } + + if (cnt != len) + { + *fLog << err << GetDescriptor() + << ": Size mismatch in weights file " << filename << endl; + return kFALSE; + } + + if (!hi) + { + *fLog << err << GetDescriptor() + << ": No correct header found in weights file " << filename << endl; + return kFALSE; + } + + return kTRUE; +} + +//---------------------------------------------------------------------------- +// +// Read a pre-defined weights file into the class. +// This is mandatory for the extraction +// +Bool_t MExtractTimeAndChargeDigitalFilter::WriteWeightsFile(TString filename) +{ + + + Float_t amp[60]; + Float_t tim[60]; + + ofstream fn(filename.Data()); + + fn << "# High Gain Weights: 6 10" << endl; + fn << "# (Amplitude) (Time) " << endl; + + for (UInt_t i=0; i<60; i++) + { + fn << "\t" << amp[i] << "\t" << tim[i] << endl; + } + + fn << "# Low Gain Weights: 6 10" << endl; + fn << "# (Amplitude) (Time) " << endl; + + for (UInt_t i=0; i<60; i++) + { + fn << "\t" << amp[i] << "\t" << tim[i] << endl; + } + return kTRUE; } @@ -238,4 +331,7 @@ fLoGainSignal.Set(range); + + fTimeShiftHiGain = (Float_t)fHiGainFirst + 0.5 + 1./fBinningResolutionHiGain; + fTimeShiftLoGain = (Float_t)fLoGainFirst + 0.5 + 1./fBinningResolutionLoGain; return kTRUE; @@ -300,12 +396,4 @@ return; -#if 0 - if (maxpos < 2) - { - time = (Float_t)fHiGainFirst; - return; - } -#endif - Float_t time_sum = 0.; Float_t fmax = 0.; @@ -334,5 +422,5 @@ for (Int_t sample=0; sample < fWindowSizeHiGain; sample++) { - const Int_t idx = fBinningResolution*sample+fBinningResolutionHalf; + const Int_t idx = fBinningResolutionHiGain*sample+fBinningResolutionHalfHiGain; const Float_t pex = fHiGainSignal[sample+i]-PedMean[(sample+i+abflag) & 0x1]; sum += fAmpWeightsHiGain [idx]*pex; @@ -351,18 +439,18 @@ { ftime_max /= fmax; - Int_t t_iter = Int_t(ftime_max*fBinningResolution); + Int_t t_iter = Int_t(ftime_max*fBinningResolutionHiGain); Int_t sample_iter = 0; - while ( t_iter > fBinningResolutionHalf-1 || t_iter < -fBinningResolutionHalf ) - { - if (t_iter > fBinningResolutionHalf-1) + while ( t_iter > fBinningResolutionHalfHiGain-1 || t_iter < -fBinningResolutionHalfHiGain ) + { + if (t_iter > fBinningResolutionHalfHiGain-1) { - t_iter -= fBinningResolution; + t_iter -= fBinningResolutionHiGain; max_p--; sample_iter--; } - if (t_iter < -fBinningResolutionHalf) + if (t_iter < -fBinningResolutionHalfHiGain) { - t_iter += fBinningResolution; + t_iter += fBinningResolutionHiGain; max_p++; sample_iter++; @@ -377,5 +465,5 @@ for (Int_t sample=0; sample < fWindowSizeHiGain; sample++) { - const Int_t idx = fBinningResolution*sample + fBinningResolutionHalf + t_iter; + const Int_t idx = fBinningResolutionHiGain*sample + fBinningResolutionHalfHiGain + t_iter; const Int_t ids = max_p + sample; const Float_t pex = ids < 0 ? 0. : @@ -383,11 +471,9 @@ 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; + time = max_p + fTimeShiftHiGain /* this shifts the time to the start of the rising edge */ + - ((Float_t)t_iter)/fBinningResolutionHiGain - time_sum/sum; else time = 0.; @@ -404,12 +490,35 @@ { - 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(); - + Int_t range = fLoGainLast - fLoGainFirst + 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++; + } + + Float_t time_sum = 0.; + Float_t fmax = 0.; + Float_t ftime_max = 0.; + Int_t max_p = 0; + + const Float_t pedes = ped.GetPedestal(); const Float_t ABoffs = ped.GetPedestalABoffset(); @@ -421,25 +530,19 @@ // 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; - + for (Int_t i=0;ifBinningResolutionHalf-1 || t_iter<-fBinningResolutionHalf) - { - if (t_iter>fBinningResolutionHalf-1) + + while ( t_iter > fBinningResolutionHalfLoGain-1 || t_iter < -fBinningResolutionHalfLoGain ) + { + if (t_iter > fBinningResolutionHalfLoGain-1) { - t_iter -= fBinningResolution; - p--; + t_iter -= fBinningResolutionLoGain; + max_p--; sample_iter--; } - if (t_iter < -fBinningResolutionHalf) + if (t_iter < -fBinningResolutionHalfLoGain) { - t_iter += fBinningResolution; - p++; + t_iter += fBinningResolutionLoGain; + max_p++; sample_iter++; } } - sum = 0; - Byte_t *p2 = p; - + sum = 0.; // // Slide with a window of size fWindowSizeLoGain over the sample @@ -484,15 +584,15 @@ 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 += fAmpWeightsLoGain [idx]*(*p2-PedMean[adx]); - time_sum += fTimeWeightsLoGain[idx]*(*p2-PedMean[adx]); - p2++; + const Int_t idx = fBinningResolutionLoGain*sample + fBinningResolutionHalfLoGain + t_iter; + const Int_t ids = max_p + sample; + const Float_t pex = ids < 0 ? 0. : + ( ids > range ? 0. : fLoGainSignal[ids]-PedMean[(ids+abflag) & 0x1]); + sum += fAmpWeightsLoGain [idx]*pex; + time_sum += fTimeWeightsLoGain[idx]*pex; } - if (sum != 0) - time = (Float_t)(max_p - ptr) + 1. + sample_iter - - (Float_t)t_iter/fBinningResolution - 0.4 - time_sum/sum - + (Float_t)fLoGainFirst; + if (sum != 0.) + time = max_p + fTimeShiftLoGain /* this shifts the time to the start of the rising edge */ + - ((Float_t)t_iter)/fBinningResolutionLoGain - time_sum/sum; else time = 0.; @@ -527,5 +627,6 @@ if (IsEnvDefined(env, prefix, "BinningResolution", print)) { - SetBinningResolution(GetEnvValue(env, prefix, "BinningResolution", fBinningResolution)); + SetBinningResolution(GetEnvValue(env, prefix, "BinningResolutionHiGain", fBinningResolutionHiGain), + GetEnvValue(env, prefix, "BinningResolutionLoGain", fBinningResolutionLoGain)); rc = kTRUE; } @@ -536,2 +637,18 @@ } +Int_t MExtractTimeAndChargeDigitalFilter::PostProcess() +{ + + *fLog << endl; + *fLog << inf << "Used High Gain weights in the extractor: " << endl; + + for (Int_t i=0;i