31 #include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h" 32 #include "larevt/CalibrationDBI/Interface/ChannelStatusService.h" 44 : FitCache{
"GausFitCache_CCHitFinderAlg"}
46 if (pset.has_key(
"MinSigInd"))
48 <<
"CCHitFinderAlg: Using no-longer-valid fcl input: MinSigInd, MinSigCol, etc";
50 fMinPeak = pset.get<std::vector<float>>(
"MinPeak");
51 fMinRMS = pset.get<std::vector<float>>(
"MinRMS");
52 fMaxBumps = pset.get<
unsigned short>(
"MaxBumps");
55 fChiNorms = pset.get<std::vector<float>>(
"ChiNorms");
58 fStudyHits = pset.get<
bool>(
"StudyHits",
false);
60 fUWireRange = pset.get<std::vector<short>>(
"UWireRange");
61 fUTickRange = pset.get<std::vector<short>>(
"UTickRange");
62 fVWireRange = pset.get<std::vector<short>>(
"VWireRange");
63 fVTickRange = pset.get<std::vector<short>>(
"VTickRange");
64 fWWireRange = pset.get<std::vector<short>>(
"WWireRange");
65 fWTickRange = pset.get<std::vector<short>>(
"WTickRange");
67 if (
fMinPeak.size() != fMinRMS.size()) {
77 <<
"CCHitFinder algorithm is currently hard-coded to support at most " <<
MaxGaussians 78 <<
" bumps per region of interest, but " << fMaxBumps <<
" have been requested.\n" 80 <<
". If this is not acceptable, increase CCHitFinderAlg::MaxGaussians" 81 <<
" value and recompile.";
90 if (fUWireRange.size() != 2 || fUTickRange.size() != 2 || fVWireRange.size() != 2 ||
91 fVTickRange.size() != 2 || fWWireRange.size() != 2 || fWTickRange.size() != 2) {
92 mf::LogError(
"CCHF") <<
"Invalid vector size for StudyHits. Must be 2";
102 constexpr
unsigned short maxticks = 1000;
103 float*
ticks =
new float[maxticks];
105 for (
unsigned short ii = 0; ii < maxticks; ++ii) {
108 float* signl =
new float[maxticks];
114 lariov::ChannelStatusProvider
const& channelStatus =
117 for (
size_t wireIter = 0; wireIter < Wires.size(); wireIter++) {
121 if (channelStatus.IsBad(
theChannel))
continue;
138 std::vector<float> signal(theWire.
Signal());
140 unsigned short nabove = 0;
141 unsigned short tstart = 0;
142 unsigned short maxtime = signal.size() - 2;
144 unsigned short mintime = 3;
145 for (
unsigned short time = 3; time < maxtime; ++time) {
151 for (
unsigned short time = mintime; time < maxtime; ++time) {
153 if (nabove == 0) tstart = time;
158 if (nabove > minSamples) {
160 if (nabove > maxticks)
162 <<
"Long RAT " << nabove <<
" " << maxticks <<
" No signal on wire " <<
theWireNum 163 <<
" after time " << time;
164 if (nabove > maxticks)
break;
165 unsigned short npt = 0;
169 for (
unsigned short ii = tstart; ii < time; ++ii) {
170 signl[npt] = signal[ii];
171 adcsum += signl[npt];
172 if (signal[ii] > signal[ii - 1] && signal[ii - 1] > signal[ii - 2] &&
173 signal[ii] > signal[ii + 1] && signal[ii + 1] > signal[ii + 2])
174 bumps.push_back(npt);
182 StoreHits(tstart, npt, WireInfo, adcsum);
187 unsigned short nHitsFit =
bumps.size();
188 unsigned short nfit = 0;
191 bool HitStored =
false;
195 while (nHitsFit <= nMaxFit) {
197 FitNG(nHitsFit, npt, ticks, signl);
204 StoreHits(tstart, npt, WireInfo, adcsum);
212 if (!HitStored && npt < maxticks) {
215 StoreHits(tstart, npt, WireInfo, adcsum);
217 else if (nHitsFit > 0)
237 std::array<double, 3>& params,
238 std::array<double, 3>& paramerrors,
246 const float time_shift = (ticks[npt - 1] + ticks[0]) / 2.F;
249 for (
size_t i = 0; i < npt; ++i) {
252 <<
"Non-positive charge encountered. Backing up to ROOT fit.";
256 fitter.
add(ticks[i] - time_shift, signl[i]);
262 MF_LOG_DEBUG(
"CCHitFinderAlg") <<
"Fast Gaussian fit failed.";
269 chidof = chi2 / fitter.
NDF();
272 params[1] += time_shift;
277 for (
double&
par : paramerrors)
278 par *= std::sqrt(chidof);
288 dof = npt - 3 * nGaus;
293 if (
bumps.size() == 0)
return;
296 std::vector<double> partmp;
297 std::vector<double> partmperr;
307 std::array<double, 3> params, paramerrors;
314 std::copy(params.begin(), params.end(), partmp.begin());
316 std::copy(paramerrors.begin(), paramerrors.end(), partmperr.begin());
331 std::stringstream numConv;
332 std::string eqn =
"gaus";
333 if (nGaus > 1) eqn =
"gaus(0)";
334 for (
unsigned short ii = 3; ii < nGaus * 3; ii += 3) {
335 eqn.append(
" + gaus(");
338 eqn.append(numConv.str());
342 auto Gn = std::make_unique<TF1>(
"gn", eqn.c_str());
343 TGraph* fitn =
new TGraph(npt, ticks, signl);
346 for (
unsigned short ii = 0; ii <
bumps.size(); ++ii) {
347 unsigned short index = ii * 3;
348 unsigned short bumptime =
bumps[ii];
349 double amp = signl[bumptime];
350 Gn->SetParameter(index, amp);
351 Gn->SetParLimits(index, 0., 9999.);
352 Gn->SetParameter(index + 1, (
double)bumptime);
353 Gn->SetParLimits(index + 1, 0, (
double)npt);
355 Gn->SetParLimits(index + 2, 1., 3 * (
double)
fMinRMS[thePlane]);
359 for (
unsigned short ii =
bumps.size(); ii < nGaus; ++ii) {
362 unsigned short imbig = 0;
363 for (
unsigned short jj = 0; jj < npt; ++jj) {
364 float diff = signl[jj] - Gn->Eval((Double_t)jj, 0, 0, 0);
372 unsigned short index = ii * 3;
373 Gn->SetParameter(index, (
double)big);
374 Gn->SetParLimits(index, 0., 9999.);
375 Gn->SetParameter(index + 1, (
double)imbig);
376 Gn->SetParLimits(index + 1, 0, (
double)npt);
378 Gn->SetParLimits(index + 2, 1., 5 * (
double)
fMinRMS[thePlane]);
384 fitn->Fit(&*Gn,
"WNQB");
386 for (
unsigned short ipar = 0; ipar < 3 * nGaus; ++ipar) {
387 partmp.push_back(Gn->GetParameter(ipar));
388 partmperr.push_back(Gn->GetParError(ipar));
399 std::vector<std::pair<unsigned short, unsigned short>> times;
401 for (
unsigned short ii = 0; ii < nGaus; ++ii) {
402 unsigned short index = ii * 3;
403 times.push_back(std::make_pair(partmp[index + 1], ii));
405 std::sort(times.begin(), times.end());
408 for (
unsigned short ii = 0; ii < nGaus; ++ii) {
409 if (times[ii].
second != ii) {
416 std::vector<double> partmpt;
417 std::vector<double> partmperrt;
418 for (
unsigned short ii = 0; ii < nGaus; ++ii) {
419 unsigned short index = times[ii].second * 3;
420 partmpt.push_back(partmp[index]);
421 partmpt.push_back(partmp[index + 1]);
422 partmpt.push_back(partmp[index + 2]);
423 partmperrt.push_back(partmperr[index]);
424 partmperrt.push_back(partmperr[index + 1]);
425 partmperrt.push_back(partmperr[index + 2]);
428 partmperr = partmperrt;
434 for (
unsigned short ii = 0; ii < nGaus; ++ii) {
435 unsigned short index = ii * 3;
437 short fittime = partmp[index + 1];
438 if (fittime < 0 || fittime > npt - 1) {
448 float rms = partmp[index + 2];
449 if (rms < 0.5 *
fMinRMS[thePlane] || rms > 5 *
fMinRMS[thePlane]) {
454 for (
unsigned short jj = 0; jj < nGaus; ++jj) {
455 if (jj == ii)
continue;
456 unsigned short jndex = jj * 3;
457 float timediff =
std::abs(partmp[jndex + 1] - partmp[index + 1]);
482 for (
unsigned short ii = 0; ii < npt; ++ii) {
484 sumST += signl[ii] * ticks[ii];
486 float mean = sumST / sumS;
488 for (
unsigned short ii = 0; ii < npt; ++ii) {
489 float arg = ticks[ii] -
mean;
490 rms += signl[ii] * arg * arg;
492 rms = std::sqrt(rms / sumS);
493 float amp = sumS / (
Sqrt2Pi * rms);
501 float meanerr = std::sqrt(1 / sumS);
502 float rmserr = 0.2 * rms;
504 parerr.push_back(meanerr);
517 size_t nhits =
par.size() / 3;
520 mf::LogError(
"CCHitFinder") <<
"Too many hits: existing " <<
allhits.size() <<
" plus new " 521 << nhits <<
" beyond the maximum " <<
allhits.max_size();
525 if (nhits == 0)
return;
530 const float loTime = TStart;
531 const float hiTime = TStart + npt;
535 for (
size_t hit = 0;
hit < nhits; ++
hit) {
536 const unsigned short index = 3 *
hit;
539 for (
size_t hit = 0;
hit < nhits; ++
hit) {
540 const size_t index = 3 *
hit;
542 const float charge_err =
548 par[index + 1] + TStart,
553 adcsum * charge / gsum,
554 adcsum * charge / gsum,
572 unsigned short tstart)
589 for (
unsigned short ipl = 0; ipl < 3; ++ipl) {
645 for (
unsigned short ii = 0; ii < npt; ++ii) {
646 if (signl[ii] > big) {
670 for (
unsigned short ii = 0; ii < npt; ++ii) {
672 sumt += signl[ii] * ii;
674 float aveb = sumt /
sum;
677 for (
unsigned short ii = 0; ii < npt; ++ii) {
678 float dbin = (float)ii - aveb;
679 sumt += signl[ii] * dbin * dbin;
689 if (
par.size() == 3) {
701 std::cout <<
"Check lo and hi W/T for each plane" << std::endl;
702 for (
unsigned short ipl = 0; ipl < 3; ++ipl) {
703 std::cout << ipl <<
" lo " <<
loWire[ipl] <<
" " <<
loTime[ipl] <<
" hi " <<
hiWire[ipl]
704 <<
" " <<
hiTime[ipl] << std::endl;
706 std::cout <<
" ipl nRAT bCnt bChi bRMS hCnt hRMS dT/dW New_ChiNorm" << std::endl;
707 for (
unsigned short ipl = 0; ipl < 3; ++ipl) {
714 std::cout << ipl <<
std::right << std::setw(5) <<
RATCnt[ipl] << std::setw(5)
715 <<
bumpCnt[ipl] << std::setw(7) << std::fixed << std::setprecision(2)
717 << std::setw(7) << std::setprecision(1) <<
hitRMS[ipl] << std::setw(7) << dTdW
722 std::cout <<
"nRAT is the number of Regions Above Threshold (RAT) used in the study.\n";
723 std::cout <<
"bCnt is the number of single bumps that were successfully fitted \n";
724 std::cout <<
"bChi is the average chisq/DOF of the first fit\n";
725 std::cout <<
"bRMS is the average calculated RMS of the bumps\n";
726 std::cout <<
"hCnt is the number of RATs that have a single hit\n";
727 std::cout <<
"hRMS is the average RMS from the Gaussian fit -> use this value for " 728 "fMinRMS[plane] in the fcl file\n";
729 std::cout <<
"dTdW is the slope of the track\n";
731 <<
"New_ChiNorm is the recommended values of ChiNorm that should be used in the fcl file\n";
748 if (nGaus == 0)
return;
749 MultiGausFits.resize(nGaus);
750 std::fill(MultiGausFits.begin(), MultiGausFits.end(), 0);
756 ++MultiGausFits[std::min(nGaus, (
unsigned int)MultiGausFits.size()) - 1];
void RunCCHitFinder(std::vector< recob::Wire > const &Wires)
std::vector< short > fWTickRange
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
bool add(Data_t x, Data_t y, Data_t sy=Data_t(1.0))
Clears all the input statistics.
void AddMultiGaus(unsigned int nGaus)
float fChiSplit
Estimated noise error on the Signal.
static constexpr float SqrtPi
void FitNG(unsigned short nGaus, unsigned short npt, float *ticks, float *signl)
raw::ChannelID_t theChannel
std::vector< float > loWire
void StoreHits(unsigned short TStart, unsigned short npt, HitChannelInfo_t info, float adcsum)
std::vector< float > hiWire
std::vector< unsigned short > bumps
std::vector< float > fChiNorms
std::vector< float > bumpRMS
bool fUseFastFit
whether to attempt using a fast fit on single gauss.
constexpr auto abs(T v)
Returns the absolute value of the argument.
CCHitFinderAlg(fhicl::ParameterSet const &pset)
SigType_t SignalType(PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
FitStats_t TriedFitStats
counts of the tried fits
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
tick ticks
Alias for common language habits.
Classes performing simple fits.
geo::View_t View() const
Returns the view the channel belongs to.
virtual bool FillResults(FitParameters_t ¶ms, FitMatrix_t &Xmat, Data_t &det, FitMatrix_t &Smat) const override
Fills the specified parameters.
std::vector< short > fUTickRange
std::vector< recob::Hit > allhits
raw::ChannelID_t Channel() const
Returns the ID of the channel (or InvalidChannelID)
geo::WireReadoutGeom const * wireReadoutGeom
void Reset(unsigned int nGaus)
void fill(const art::PtrVector< recob::Hit > &hits, int only_plane)
std::vector< float > hiTime
std::vector< short > fUWireRange
virtual int NDF() const override
Returns the degrees of freedom in the determination of the fit.
std::vector< int > bumpCnt
Detector simulation of raw signals on wires.
std::vector< float > Signal() const
Return a zero-padded full length vector filled with RoI signal.
std::vector< float > loTime
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
static bool FastGaussianFit(unsigned short npt, float const *ticks, float const *signl, std::array< double, 3 > ¶ms, std::array< double, 3 > ¶merrors, float &chidof)
Performs a "fast" fit.
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
void StudyHits(unsigned short flag, unsigned short npt=0, float *signl=nullptr, unsigned short tstart=0)
virtual std::vector< WireID > ChannelToWire(raw::ChannelID_t channel) const =0
std::vector< short > fWWireRange
unsigned short theWireNum
std::vector< float > fMinPeak
Class holding the regions of interest of signal from a channel.
std::vector< int > hitCnt
static constexpr unsigned int MaxGaussians
static constexpr float Sqrt2Pi
std::vector< short > fVWireRange
std::vector< short > fVTickRange
virtual Data_t ChiSquare() const override
Returns the of the original fit.
unsigned short fMaxXtraHits
void MakeCrudeHit(unsigned short npt, float *ticks, float *signl)
std::vector< float > hitRMS
#define MF_LOG_WARNING(category)
std::vector< int > RATCnt
second_as<> second
Type of time stored in seconds, in double precision.
FitStats_t FinalFitStats
counts of the good fits
std::vector< double > parerr
std::vector< float > fMinRMS
std::vector< float > bumpChi
exchange data about the originating wire
std::vector< double > par
art framework interface to geometry description