diff --git a/Scintillator/ScintRecAlgs/python/ScintRecAlgsConfig.py b/Scintillator/ScintRecAlgs/python/ScintRecAlgsConfig.py index 6319e9a01771a042f162c11f51d4f5b8d3bba181..f96dbcaf2f7ba283437116bfad0cd244af8220c6 100644 --- a/Scintillator/ScintRecAlgs/python/ScintRecAlgsConfig.py +++ b/Scintillator/ScintRecAlgs/python/ScintRecAlgsConfig.py @@ -29,7 +29,7 @@ def WaveformClockRecCfg(flags, name="ClockRecAlg", **kwargs): acc = ComponentAccumulator() tool = ClockReconstructionTool(name="ClockReconstructionTool") - + # tool.CheckResult = True kwargs.setdefault("ClockReconstructionTool", tool) recoAlg = CompFactory.ScintClockRecAlg(name, **kwargs) diff --git a/Scintillator/ScintRecTools/src/ClockReconstructionTool.cxx b/Scintillator/ScintRecTools/src/ClockReconstructionTool.cxx index 15a3e5f3b9452552163bec492cc6655625f04062..c4d1bf68f99a32d9e940c38de62bae9986e7464e 100644 --- a/Scintillator/ScintRecTools/src/ClockReconstructionTool.cxx +++ b/Scintillator/ScintRecTools/src/ClockReconstructionTool.cxx @@ -64,7 +64,8 @@ ClockReconstructionTool::reconstruct(const ScintWaveform& raw_wave, ATH_MSG_DEBUG("Created double array with length " << wave.size() ); ATH_MSG_DEBUG("First 10 elements:"); - for (unsigned int i=0; i < std::min(10, N); i++) + + for (int i=0; i < std::min(10, N); i++) ATH_MSG_DEBUG(" " << i << " " << wave[i]); // delta_nu = 1/T where T is the total waveform length @@ -138,5 +139,34 @@ ClockReconstructionTool::reconstruct(const ScintWaveform& raw_wave, delete fftr2c; + // Check the result for sanity + if (clockdata->amplitude() < m_amplitude_min || clockdata->frequency() < m_frequency_min || clockdata->frequency() > m_frequency_max) { + ATH_MSG_WARNING("Clock reconstructed with frequency = " << clockdata->frequency() << " amplitude = " << clockdata->amplitude() << "!"); + } + + if (m_checkResult) checkResult(raw_wave, clockdata); + return StatusCode::SUCCESS; } + +void +ClockReconstructionTool::checkResult(const ScintWaveform& raw_wave, + xAOD::WaveformClock* clockdata) const { + + // Go through each element in raw_wave and make sure time in clockdata matches + + float time; + for (unsigned int i=0; i<raw_wave.adc_counts().size(); i++) { + time = 2.*i; // Time in ns + + float dt = clockdata->time_from_clock(time); + + // Is raw_wave HI or LO? + bool hi = raw_wave.adc_counts()[i] > clockdata->dc_offset(); + + // Check for mismatch + if (((dt > 0.05) && (dt < 12.45) && !hi) || ((dt > 12.55) && (dt < 24.95) && hi) ) + ATH_MSG_WARNING("Clock Time: " << time << " Time to edge: " << dt << " found raw data: " << hi); + } + +} diff --git a/Scintillator/ScintRecTools/src/ClockReconstructionTool.h b/Scintillator/ScintRecTools/src/ClockReconstructionTool.h index 0b52229e22811b04849664565d167f0706a611f1..056d1a4e9e52b1a7ee6963a45b45c5872a46db99 100644 --- a/Scintillator/ScintRecTools/src/ClockReconstructionTool.h +++ b/Scintillator/ScintRecTools/src/ClockReconstructionTool.h @@ -45,6 +45,17 @@ class ClockReconstructionTool: public extends<AthAlgTool, IClockReconstructionTo /// Minimum samples in the input waveform array to try FFT IntegerProperty m_minimumSamples{this, "MinimumSamples", 40}; + /// Check reconstructed clock against waveform + BooleanProperty m_checkResult{this, "CheckResult", false}; + + void checkResult(const ScintWaveform& raw_wave, + xAOD::WaveformClock* clockdata) const; + + // Limits to print warnings + FloatProperty m_amplitude_min{this, "AmplitudeMin", 1000.}; + FloatProperty m_frequency_min{this, "FrequencyMin", 40.0}; + FloatProperty m_frequency_max{this, "FrequencyMax", 40.1}; + }; #endif // SCINTRECTOOLS_CLOCKRECONSTRUCTIONTOOL_H diff --git a/Scintillator/ScintRecTools/src/WaveformReconstructionTool.cxx b/Scintillator/ScintRecTools/src/WaveformReconstructionTool.cxx index 7b442b6950307f7da9d39de98cf037c5ce82bf32..e774f918de2e836faee2fa583258fcfbda1f0da2 100644 --- a/Scintillator/ScintRecTools/src/WaveformReconstructionTool.cxx +++ b/Scintillator/ScintRecTools/src/WaveformReconstructionTool.cxx @@ -20,6 +20,7 @@ #include <vector> #include <tuple> +#include <math.h> // Constructor WaveformReconstructionTool::WaveformReconstructionTool(const std::string& type, const std::string& name, const IInterface* parent) : @@ -168,8 +169,13 @@ WaveformReconstructionTool::reconstruct(const ScintWaveform& raw_wave, // Perform Gaussian fit to waveform WaveformFitResult gfit = fitGaussian(raw, wtime, wwave); if (! gfit.valid) { - ATH_MSG_WARNING( " Gaussian waveform fit failed! "); - hit->set_status_bit(xAOD::WaveformStatus::GFIT_FAILED); + // Lets try again with a more restricted width + ATH_MSG_WARNING( " Gaussian waveform fit failed with width " << raw.sigma << " try reducing width to 1 " ); + raw.sigma = 1.; + gfit = fitGaussian(raw, wtime, wwave); + if (!gfit.valid) { + hit->set_status_bit(xAOD::WaveformStatus::GFIT_FAILED); + } } // Fit results (or raw if it failed) @@ -455,6 +461,12 @@ WaveformReconstructionTool::findRawHitValues(const std::vector<float> time, cons ATH_MSG_DEBUG( "Initial Mean: " << rfit.mean << " RMS: " << rfit.sigma << " Peak: " << rfit.peak); + // Fix bad values + if (isnan(rfit.sigma)) { + ATH_MSG_DEBUG("Found sigma: " << rfit.sigma << " Replace with 2."); + rfit.sigma = 2.; + } + return rfit; } @@ -479,7 +491,7 @@ WaveformReconstructionTool::fitGaussian(const WaveformFitResult& raw, const std: gfunc.SetParameters(raw.peak, raw.mean, raw.sigma); gfunc.SetParError(0, std::sqrt(raw.peak)); gfunc.SetParError(1, raw.sigma); - gfunc.SetParError(2, raw.sigma / 10.); + gfunc.SetParError(2, raw.sigma / 5.); gfunc.SetParLimits(2, 0., 20.); // Constrain width TFitResultPtr gfitptr = tg.Fit(&gfunc, "QNS", ""); @@ -523,9 +535,9 @@ WaveformReconstructionTool::fitCBall(const WaveformFitResult& gfit, cbfit.peak = abs(gfit.peak); cbfit.mean = gfit.mean; cbfit.sigma = abs(gfit.sigma); - if (cbfit.sigma > 20.) cbfit.sigma = 5.; - cbfit.alpha = -0.25; // Negative to get tail on high side - cbfit.nval = 3.; + if (cbfit.sigma > 20.) cbfit.sigma = 2.; + cbfit.alpha = -0.5; // Negative to get tail on high side + cbfit.nval = 2.; // Define fit function and preset values TF1 cbfunc("cbfunc", "crystalball"); @@ -534,16 +546,27 @@ WaveformReconstructionTool::fitCBall(const WaveformFitResult& gfit, cbfunc.SetParameter(1, cbfit.mean); // Mean cbfunc.SetParError(1, cbfit.sigma); cbfunc.SetParameter(2, cbfit.sigma); // Width - cbfunc.SetParError(2, cbfit.sigma/10.); + cbfunc.SetParError(2, cbfit.sigma/5.); cbfunc.SetParLimits(2, 0., 20.); cbfunc.SetParameter(3, cbfit.alpha); // Tail parameter (negative for high-side tail) - cbfunc.SetParError(3, -cbfit.alpha/10.); + cbfunc.SetParError(3, -cbfit.alpha/5.); cbfunc.SetParLimits(3, -10., 0.); + // Try fixing the tail first + cbfunc.FixParameter(4, cbfit.nval); + + TFitResultPtr cbfitptr = tg.Fit(&cbfunc, "QNS", ""); + + if (!cbfitptr->IsValid()) { + ATH_MSG_WARNING( " First Crystal Ball waveform fit failed! "); + } + + // Now try releasing the tail parameter + cbfunc.ReleaseParameter(4); cbfunc.SetParameter(4, cbfit.nval); // Tail power - cbfunc.SetParError(4, cbfit.nval/10.); + cbfunc.SetParError(4, cbfit.nval/5.); cbfunc.SetParLimits(4, 0., 1.E3); - TFitResultPtr cbfitptr = tg.Fit(&cbfunc, "QNS", ""); + cbfitptr = tg.Fit(&cbfunc, "QNS", ""); cbfit.fit_status = cbfitptr; cbfit.valid = (cbfit.fit_status == 0); diff --git a/Scintillator/ScintRecTools/src/WaveformReconstructionTool.h b/Scintillator/ScintRecTools/src/WaveformReconstructionTool.h index 45a1e9c79d1229e337e6223090959d6b955f7b03..e8b2abf937c0b2934b6521b0786d0832240ffad3 100644 --- a/Scintillator/ScintRecTools/src/WaveformReconstructionTool.h +++ b/Scintillator/ScintRecTools/src/WaveformReconstructionTool.h @@ -80,7 +80,7 @@ class WaveformReconstructionTool: public extends<AthAlgTool, IWaveformReconstruc // // Fraction of peak to set local hit time - FloatProperty m_timingPeakFraction{this, "TimingPeakFraction", 0.2}; + FloatProperty m_timingPeakFraction{this, "TimingPeakFraction", 0.45}; // Baseline algorithms WaveformBaselineData& findSimpleBaseline(const ScintWaveform& wave) const; diff --git a/Tracker/TrackerRecAlgs/TrackerClusterFit/src/ClusterFitAlg.cxx b/Tracker/TrackerRecAlgs/TrackerClusterFit/src/ClusterFitAlg.cxx index 12fa1c7b945a077e99e6e89f1ba307a2186a0e98..8fe463790db2b8e474b911f3fa520f2e257234d0 100644 --- a/Tracker/TrackerRecAlgs/TrackerClusterFit/src/ClusterFitAlg.cxx +++ b/Tracker/TrackerRecAlgs/TrackerClusterFit/src/ClusterFitAlg.cxx @@ -377,9 +377,9 @@ ClusterFitAlg::AddTrack(std::unique_ptr<TrackCollection>& tracks, s->push_back( GetState(fitResult, fitCovariance, nullptr, station) ); for (const clusterInfo* cInfo : fitClusters) - { - s->push_back( GetState(fitResult, fitCovariance, cInfo->cluster, station) ); - } + { + s->push_back( GetState(fitResult, fitCovariance, cInfo->cluster, station) ); + } // Create and store track tracks->push_back(new Trk::Track(i, s , q)); @@ -495,4 +495,4 @@ ClusterFitAlg::Residuals(std::vector<const clusterInfo*>& fitClusters) const } } -} \ No newline at end of file +} diff --git a/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHit_v1.h b/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHit_v1.h index 976c51512ae734902db9a5b3e30a309c5610723e..8b69c776da83546cde3a787a25bcf2b3a830cf7f 100644 --- a/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHit_v1.h +++ b/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHit_v1.h @@ -105,7 +105,7 @@ namespace xAOD { this->set_status(this->status() | (1<<bit)); } bool status_bit(WaveformStatus bit) const { - return (this->status() | (1<<bit)); + return (this->status() & (1<<bit)); } bool threshold() const {