31 #include "TDecompSVD.h" 71 fAreaNorms = pset.get<std::vector<double>>(
"AreaNorms");
97 std::vector<int> startTimes;
98 std::vector<int> maxTimes;
99 std::vector<int> endTimes;
101 int minTimeHolder = 0;
103 bool maxFound =
false;
104 double threshold = 0.;
105 double fitWidth = 0.;
106 double minWidth = 0.;
107 std::string eqn =
"gaus(0)";
109 std::stringstream numConv;
112 for (
unsigned int wireIter = 0; wireIter < wireVecHandle->size(); wireIter++) {
118 std::vector<float> signal(wire->
Signal());
124 sigType = wireReadoutGeom.SignalType(channel);
138 for (timeIter = signal.begin(); timeIter + 2 < signal.end(); timeIter++) {
140 if (*timeIter > *(timeIter + 1) && *(timeIter + 1) < *(timeIter + 2)) {
143 endTimes.push_back(time + 1);
146 minTimeHolder = time + 2;
149 minTimeHolder = time + 1;
153 else if (*timeIter < *(timeIter + 1) && *(timeIter + 1) > *(timeIter + 2) &&
154 *(timeIter + 1) > threshold) {
156 maxTimes.push_back(time + 1);
157 startTimes.push_back(minTimeHolder);
163 while (maxTimes.size() > endTimes.size())
164 endTimes.push_back(signal.size() - 1);
165 if (startTimes.size() == 0)
continue;
176 double amplitude(0), position(0), width(0);
177 double amplitudeErr(0), positionErr(0), widthErr(0);
178 double goodnessOfFit(0), chargeErr(0);
179 double minPeakHeight(0);
183 std::vector<double> hitSig;
186 while (hitIndex < (
signed)startTimes.size()) {
190 minPeakHeight = signal[maxTimes[hitIndex]];
197 while (numHits <
fMaxMultiHit && numHits + hitIndex < (
signed)endTimes.size() &&
198 signal[endTimes[hitIndex + numHits - 1]] > threshold / 2.0 &&
199 startTimes[hitIndex + numHits] - endTimes[hitIndex + numHits - 1] < 2) {
201 if (signal[maxTimes[hitIndex + numHits]] < minPeakHeight)
202 minPeakHeight = signal[maxTimes[hitIndex + numHits]];
208 startT = startTimes[hitIndex];
210 while (signal[(
int)startT] < minPeakHeight / 2.0)
214 endT = endTimes[hitIndex + numHits - 1];
216 while (signal[(
int)endT] < minPeakHeight / 2.0)
218 size = (int)(endT - startT);
219 TH1D hitSignal(
"hitSignal",
"", size, startT, endT);
220 for (
int i = (
int)startT; i < (int)endT; ++i)
221 hitSignal.Fill(i, signal[i]);
226 for (
int i = 3; i < numHits * 3; i += 3) {
227 eqn.append(
"+gaus(");
230 eqn.append(numConv.str());
234 TF1 gSum(
"gSum", eqn.c_str(), 0,
size);
237 TArrayD data(numHits * numHits);
238 TVectorD amps(numHits);
239 for (
int i = 0; i < numHits; ++i) {
240 amps[i] = signal[maxTimes[hitIndex + i]];
241 for (
int j = 0; j < numHits; j++)
242 data[i + numHits * j] =
243 TMath::Gaus(maxTimes[hitIndex + j], maxTimes[hitIndex + i], fitWidth);
249 TMatrixD h(numHits, numHits);
250 h.Use(numHits, numHits, data.GetArray());
260 for (
int i = 0; i < numHits; ++i) {
267 amplitude = 0.5 * (threshold + signal[maxTimes[hitIndex + i]]);
268 gSum.SetParameter(3 * i, amplitude);
269 gSum.SetParameter(1 + 3 * i, maxTimes[hitIndex + i]);
270 gSum.SetParameter(2 + 3 * i, fitWidth);
271 gSum.SetParLimits(3 * i, 0.0, 3.0 * amplitude);
272 gSum.SetParLimits(1 + 3 * i, startT, endT);
273 gSum.SetParLimits(2 + 3 * i, 0.0, 10.0 * fitWidth);
277 gSum.SetParameters(signal[maxTimes[hitIndex]], maxTimes[hitIndex], fitWidth);
278 gSum.SetParLimits(0, 0.0, 1.5 * signal[maxTimes[hitIndex]]);
279 gSum.SetParLimits(1, startT, endT);
280 gSum.SetParLimits(2, 0.0, 10.0 * fitWidth);
284 hitSignal.Fit(&gSum,
"QNRW",
"", startT, endT);
285 for (
int hitNumber = 0; hitNumber < numHits; ++hitNumber) {
287 if (gSum.GetParameter(3 * hitNumber) > threshold / 2.0 &&
288 gSum.GetParameter(3 * hitNumber + 2) > minWidth) {
289 amplitude = gSum.GetParameter(3 * hitNumber);
290 position = gSum.GetParameter(3 * hitNumber + 1);
291 width = gSum.GetParameter(3 * hitNumber + 2);
292 amplitudeErr = gSum.GetParError(3 * hitNumber);
293 positionErr = gSum.GetParError(3 * hitNumber + 1);
294 widthErr = gSum.GetParError(3 * hitNumber + 2);
295 goodnessOfFit = gSum.GetChisquare() / (double)gSum.GetNDF();
296 int DoF = gSum.GetNDF();
299 chargeErr = std::sqrt(TMath::Pi()) * (amplitudeErr * width + widthErr * amplitude);
303 for (
int sigPos = 0; sigPos <
size; ++sigPos) {
304 hitSig[sigPos] = amplitude * TMath::Gaus(sigPos + startT, position, width);
305 totSig += hitSig[(int)sigPos];
309 totSig = std::sqrt(2 * TMath::Pi()) * amplitude * width /
fAreaNorms[(size_t)sigType];
312 std::vector<geo::WireID> wids = wireReadoutGeom.ChannelToWire(channel);
330 (signal.begin() + (int)startT, signal.begin() + (int)endT, 0.),
double fIndMinWidth
Minimum induction hit width.
FFTHitFinder(fhicl::ParameterSet const &pset)
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
double fMinSigInd
Induction signal height threshold.
EDProducer(fhicl::ParameterSet const &pset)
double fIndWidth
Initial width for induction fit.
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
int fMaxMultiHit
maximum hits for multi fit
double fColWidth
Initial width for collection fit.
static void declare_products(art::ProducesCollector &collector, std::string instance_name="", bool doWireAssns=true, bool doRawDigitAssns=true)
Declares the hit products we are going to fill.
std::vector< double > fAreaNorms
factors for converting area to same units as peak height
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Class managing the creation of a new recob::Hit object.
Helper functions to create a hit.
std::string fCalDataModuleLabel
constexpr ChannelID_t InvalidChannelID
ID of an invalid channel.
A class handling a collection of hits and its associations.
#define DEFINE_ART_MODULE(klass)
Signal from induction planes.
raw::ChannelID_t Channel() const
Returns the ID of the channel (or InvalidChannelID)
enum geo::_plane_sigtype SigType_t
Enumerate the possible plane projections.
double fColMinWidth
Minimum collection hit width.
void emplace_back(recob::Hit &&hit, art::Ptr< recob::Wire > const &wire=art::Ptr< recob::Wire >(), art::Ptr< raw::RawDigit > const &digits=art::Ptr< raw::RawDigit >())
Adds the specified hit to the data collection.
void produce(art::Event &evt) override
Definition of data types for geometry description.
void put_into(art::Event &)
Moves the data into an event.
Detector simulation of raw signals on wires.
std::vector< float > Signal() const
Return a zero-padded full length vector filled with RoI signal.
ProducesCollector & producesCollector() noexcept
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
double fMinSigCol
Collection signal height threshold.
Declaration of basic channel signal object.
int fAreaMethod
Type of area calculation.
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Signal from collection planes.