328 std::ofstream
output(
"areaFit.out");
333 art::ServiceHandle<geo::Geometry> geom;
349 art::Handle< std::vector<recob::Wire> > wireVecHandle;
355 art::FindOneP<raw::RawDigit> RawDigits
362 std::vector<float> holder;
363 std::vector<short> rawadc;
365 std::vector<float> startTimes;
366 std::vector<float> maxTimes;
367 std::vector<float> endTimes;
368 std::vector<float> peakHeight;
369 std::vector<float> hitrms;
370 std::vector<double> charge;
373 std::stringstream numConv;
375 unsigned int hwC(0),lwC(5728);
376 unsigned int hwI2(0),lwI2(5728);
377 unsigned int hwI1(0),lwI1(2112);
379 unsigned int nw1hitC(0),nw1hitI2(0),nw1hitI1(0);
381 for(
int jw=0;jw<5600;jw++)
384 for(
int jw=0;jw<5600;jw++)
387 for(
int jw=0;jw<5600;jw++)
391 unsigned int nhitsC(0),nhitsI1(0),nhitsI2(0);
393 unsigned int nnhitsC(0),nnhitsI1(0),nnhitsI2(0);
396 unsigned int minWireC,maxWireC;
397 if(
fThetaAngle==45) {minWireC=2539; maxWireC=3142;}
399 if(
fThetaAngle==20) {minWireC=2539; maxWireC=4000;}
400 if(
fThetaAngle==40) {minWireC=2539; maxWireC=3190;}
401 if(
fThetaAngle==60) {minWireC=2539; maxWireC=2905;}
402 if(
fThetaAngle==70) {minWireC=2539; maxWireC=2805;}
403 if(
fThetaAngle==80) {minWireC=2539; maxWireC=2740;}
407 = art::ServiceHandle<lariov::ChannelStatusService>()->GetProvider();
412 unsigned int minWireI2=2539;
413 unsigned int maxWireI2=4700;
414 unsigned int minDrift=850;
415 unsigned int maxDrift=1500;
419 for(
size_t wireIter = 0; wireIter < wireVecHandle->size(); wireIter++)
424 art::Ptr<recob::Wire> wire(wireVecHandle, wireIter);
425 art::Ptr<raw::RawDigit> rawdigits = RawDigits.at(wireIter);
428 channel = wire->Channel();
430 std::vector<float> signal(wire->Signal());
434 std::vector<geo::WireID> wids = geom->ChannelToWire(channel);
448 channel = rawdigits->Channel();
451 std::vector<geo::WireID> widVec = geom->ChannelToWire(channel);
452 size_t iwire = widVec[0].Wire;
461 int pedestal = (int)rawdigits->GetPedestal();
462 raw::Uncompress(rawdigits->ADCs(), rawadc, pedestal, rawdigits->Compression());
468 mf::LogDebug(
"ICARUSHitFinder") <<
" pedestal " <<rawdigits->GetPedestal() << std::endl;
475 if(plane == 0) holder[
bin]=-holder[
bin];
480 if(plane==0&&iwire<lwI1) lwI1=iwire;
481 if(plane==0&&iwire>hwI1) hwI1=iwire;
482 if(plane==1&&iwire<lwI2) lwI2=iwire;
483 if(plane==1&&iwire>hwI2) hwI2=iwire;
484 if(plane==2&&iwire<lwC) lwC=iwire;
485 if(plane==2&&iwire>hwC) hwC=iwire;
498 bool channelSwitch =
false;
500 for(
auto it = BadChannels.begin(); it != BadChannels.end(); it++)
504 channelSwitch =
true;
509 if(channelSwitch==
false)
527 std::vector<geo::WireID> wids = geom->ChannelToWire(channel);
529 std::vector<float> tempVec = holder;
530 recob::Wire::RegionsOfInterest_t::datarange_t rangeData(
size_t(0),std::move(tempVec));
532 fHitFinderTool->findHitCandidates(rangeData, 0,channel,0,hitCandidateVec);
534 for(
auto& hitCand : hitCandidateVec) {
535 expandHit(hitCand,holder,hitCandidateVec);
541 fHitFinderTool->MergeHitCandidates(rangeData, hitCandidateVec, mergedCandidateHitVec);
556 for(
auto& mergedCands : mergedCandidateHitVec)
567 mf::LogDebug(
"ICARUSHitFinder") <<
" adding localmean " << mean << std::endl;
571 if (endT - startT < 5)
continue;
578 int nGausForFit = mergedCands.size();
583 double chi2PerNDF(0.);
587 const int npk=mergedCands.size();
589 peakParamsVec.clear();
600 if (!(chi2PerNDF < std::numeric_limits<double>::infinity()))
610 if (chi2PerNDF < 10.)
611 fChi2->Fill(chi2PerNDF);
615 peakParamsLong.clear();
621 if(chi2Long<chi2PerNDF&&chi2Long>0.1) {
622 fChi2->Fill(chi2Long);
623 peakParamsVec=peakParamsLong;
625 else {
fChi2->Fill(chi2PerNDF);
635 for(
unsigned int jhit=0;jhit<mergedCands.size(); jhit++)
642 if(jhit>=1&&startInt<mergedCands[jhit-1].stopTick) startInt=mergedCands[jhit-1].stopTick;
643 if(jhit<mergedCands.size()-1&&endInt>mergedCands[jhit+1].startTick) endInt=mergedCands[jhit+1].startTick;
646 float peakAmp, peakMean, peakLeft, peakRight,
peakBaseline;
648 float peakMeanErr, peakAmpErr;
653 Func.SetParameter(0, mergedCands.size());
660 peakAmp = peakParams.peakAmplitude;
661 peakMean = peakParams.peakCenter;
663 peakLeft = peakParams.peakTauLeft;
664 peakRight = peakParams.peakTauRight;
665 peakBaseline = peakParams.peakBaseline;
669 if (std::isnan(peakAmp))
676 peakAmpErr = peakParams.peakAmplitudeError;
677 peakMeanErr = peakParams.peakCenterError;
679 Func.SetParameter(1+5*jhit,peakBaseline);
680 Func.SetParameter(2+5*jhit,peakAmp);
681 Func.SetParameter(3+5*jhit,peakMean);
682 Func.SetParameter(4+5*jhit,peakRight);
683 Func.SetParameter(5+5*jhit,peakLeft);
685 intBaseline+=(endInt-startInt)*peakBaseline;
689 fitCharge=Func.Integral(startInt,endInt)-(endInt-startInt)*
localmeans[jhit];
693 mf::LogWarning(
"ICARUSHitFinder") <<
"Icarus numerical 32 failed";
694 fitCharge=std::accumulate(holder.begin() + (int) startInt, holder.begin() + (int) endInt, 0.)-(endInt-startInt)*
localmeans[jhit];
701 FuncLong.SetParameter(0, mergedCands.size());
709 peakAmp = peakParams.peakAmplitude;
710 peakMean = peakParams.peakCenter;
712 peakLeft = peakParams.peakTauLeft;
713 peakRight = peakParams.peakTauRight;
714 peakBaseline = peakParams.peakBaseline;
715 intBaseline+=(endInt-startInt)*peakBaseline;
716 peakAmpErr = peakParams.peakAmplitudeError;
717 peakMeanErr = peakParams.peakCenterError;
719 FuncLong.SetParameter(1+7*jhit,peakBaseline);
720 FuncLong.SetParameter(2+7*jhit,peakAmp);
721 FuncLong.SetParameter(3+7*jhit,peakMean);
722 FuncLong.SetParameter(4+7*jhit,peakRight);
723 FuncLong.SetParameter(5+7*jhit,peakLeft);
725 FuncLong.SetParameter(7+7*jhit,peakSlope);
729 { fitCharge=FuncLong.Integral(startInt,endInt)-(endInt-startInt)*
localmeans[jhit];
733 {mf::LogWarning(
"ICARUSHitFinder") <<
"Icarus numerical integration failed";
734 fitCharge=std::accumulate(holder.begin() + (int) startInt, holder.begin() + (int) endInt, 0.)-(endInt-startInt)*
localmeans[jhit];
737 if(isnan(fitCharge)&&!islong) fitCharge=std::accumulate(holder.begin() + (int) startInt, holder.begin() + (int) endInt, 0.);
738 if(isnan(fitCharge)&&islong) fitCharge=std::accumulate(holder.begin() + (int) startInt, holder.begin() + (int) endInt, 0.);
743 float totSig=std::accumulate(holder.begin()+ (int) startInt, holder.begin()+ (int) endInt, 0.)-(endInt-startInt)*
localmeans[jhit];
756 (peakLeft+peakRight)/2.,
770 mf::LogDebug(
"ICARUSHitFinder") <<
" fitcharge " << fitCharge <<
" totSig " << totSig << std::endl;
774 hcol.emplace_back(
hit.move(), wire, rawdigits);
776 bool wireWindowC=iWire>=minWireC&&iWire<=maxWireC;
777 bool outWireWindowC=iWire<=minWireC||iWire>=maxWireC;
781 if(nghC==1) nw1hitC++;
786 wCharge[
iWire]+=fitCharge;
794 if(tpc==0&&cryostat==0&&outWireWindowC) {
795 nnhitsC++;
fNoiseC->Fill(peakAmp); }
799 if(cryostat==0&&tpc==0)
801 fAreaC->Fill(wCharge[iwire]);
802 if(cryostat==0&&tpc==0)
805 if(cryostat==0&&tpc==0)
807 output << iwire <<
" " <<wInt[iwire] << std::endl;
808 if(wInt[iwire]>0&&cryostat==0&&tpc==0)
809 fAreaInt->Fill(wCharge[iwire]/wInt[iwire]);
818 if(plane==0||plane==1) {
819 for(
auto& mergedCands : mergedCandidateHitVec)
821 for(
size_t jh=0;jh<mergedCands.size();jh++) {
827 mergedCands[jh].startTick,
828 mergedCands[jh].stopTick,
829 mergedCands[jh].hitSigma,
830 mergedCands[jh].hitCenter,
832 mergedCands[jh].hitHeight,
834 std::accumulate(holder.begin() + (int) mergedCands[jh].startTick, holder.begin() + (int) mergedCands[jh].stopTick, 0.),
836 std::accumulate(holder.begin() + (int) mergedCands[jh].startTick, holder.begin() + (int) mergedCands[jh].stopTick, 0.),
840 int(mergedCands[jh].stopTick-mergedCands[jh].startTick+1)
842 hcol.emplace_back(
hit.move(), wire, rawdigits);
844 bool driftWindow=(mergedCands[jh].hitCenter)>=minDrift&&(mergedCands[jh].hitCenter)<=maxDrift;
845 bool wireWindowI2=iWire>=minWireI2&&iWire<=maxWireI2;
846 bool outDriftWindow=(mergedCands[jh].hitCenter)<=minDrift||(mergedCands[jh].hitCenter)>=maxDrift;
847 bool outWireWindowI2=iWire<=minWireI2||iWire>=maxWireI2;
850 if(plane==0&&driftWindow) {
853 if(nghI1==1) nw1hitI1++;
856 fHeightI1->Fill(mergedCands[jh].hitHeight);
857 fWidthI1->Fill(mergedCands[jh].hitSigma);
859 if(plane==1&&wireWindowI2) {
862 if(nghI2==1) nw1hitI2++;
864 fHeightI2->Fill(mergedCands[jh].hitHeight);
865 fWidthI2->Fill(mergedCands[jh].hitSigma);
869 if(plane==0) nhitsI1++;
870 if(plane==1) nhitsI2++;
872 if(plane==0&&tpc==0&&cryostat==0&&outDriftWindow) {nnhitsI1++;
fNoiseI1->Fill(mergedCands[jh].hitHeight); }
873 if(plane==1&&tpc==0&&cryostat==0&&outWireWindowI2) { nnhitsI2++;
fNoiseI2->Fill(mergedCands[jh].hitHeight); }
886 for(
unsigned int jw=minWireC;jw<maxWireC;jw++)
887 fnhwC->Fill(nhWire[jw]);
size_t fMaxMultiHit
maximum hits for multi fit
std::set< raw::ChannelID_t > ChannelSet_t
Type of set of channel IDs.
unsigned int PlaneID_t
Type for the ID number.
struct ICARUSPeakFitParams{float peakCenter ICARUSPeakFitParams_t
CryostatID_t Cryostat
Index of cryostat.
WireID_t Wire
Index of the wire within its plane.
virtual ChannelSet_t BadChannels() const =0
Returns a copy of set of bad channel IDs for the current run.
void computeBestLocalMean(std::vector< reco_tool::ICandidateHitFinder::HitCandidate > h, std::vector< float > holder, reco_tool::ICandidateHitFinder::MergeHitCandidateVec how, float &localmean)
constexpr details::BinObj< T > bin(T value)
Returns a wrapper to print the specified data in binary format.
double ComputeNullChiSquare(std::vector< float >) const
ICARUSlongHitFitCache fLongFitCache
Cached functions for long hits.
Class managing the creation of a new recob::Hit object.
constexpr ChannelID_t InvalidChannelID
ID of an invalid channel.
A class handling a collection of hits and its associations.
int fFittingRange
semi-width of interval where to fit hit
std::unique_ptr< reco_tool::ICandidateHitFinder > fHitFinderTool
For finding candidate hits.
std::vector< double > localmeans
Class providing information about the quality of channels.
virtual TF1 * Get(size_t nFunc)
Returns a function sum of nFunc base functions.
PlaneID_t Plane
Index of the plane within its TPC.
ICARUShitFitCache fFitCache
Cached functions for multi-peak fits.
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
std::vector< ICARUSPeakFitParams_t > ICARUSPeakParamsVec
BEGIN_PROLOG sequence::SlidingWindowTriggerPatternsOppositeWindows END_PROLOG simSlidingORM6O6 effSlidingORW output
void expandHit(reco_tool::ICandidateHitFinder::HitCandidate &h, std::vector< float > holder, std::vector< reco_tool::ICandidateHitFinder::HitCandidate > how)
unsigned int ChannelID_t
Type representing the ID of a readout channel.
void findLongPeakParameters(const std::vector< float > &, const reco_tool::ICandidateHitFinder::HitCandidateVec &, ICARUSPeakParamsVec &, double &, int &, int) const
TPCID_t TPC
Index of the TPC within its cryostat.
void Uncompress(const std::vector< short > &adc, std::vector< short > &uncompressed, raw::Compress_t compress)
Uncompresses a raw data buffer.
int fIntegratingRange
semi-width of interval where to integrate fitting function
void findMultiPeakParameters(const std::vector< float > &, const reco_tool::ICandidateHitFinder::HitCandidateVec &, ICARUSPeakParamsVec &, double &, int &, int) const
std::string fCalDataModuleLabel