20 #include "canvas/Persistency/Common/FindOneP.h" 21 #include "canvas/Persistency/Common/FindManyP.h" 27 #include "art_root_io/TFileService.h" 32 #include "nurandom/RandomUtils/NuRandomService.h" 72 #include "CLHEP/Random/RandomEngine.h" 82 class RobustHitFinder;
119 bool ValidTrigger(std::vector<unsigned int> evtTriggers,
unsigned int & c1arg,
unsigned int & c2arg,
unsigned int & trignumarg);
120 float TimeToDriftDist(
float thistime,
unsigned int thistpc);
121 float TimeToDisplacement(
float thistime);
122 float TimeToX(
float thistime,
unsigned int thistpc);
123 float hitGeomDist(TVector3 hitloc, TVector3 trigloc1, TVector3 trigloc2);
278 fClks(lar::providerFrom<detinfo::DetectorClocksService>()),
279 fDetProp(lar::providerFrom<detinfo::DetectorPropertiesService>()),
299 mf::LogError(
"RobustHitFinder") <<
"No recob::Wires found when expected";
311 mf::LogError(
"RobustHitFinder") <<
"No T0s found in this event";
324 for (
size_t iwire = 0; iwire < wireHandle->size(); ++iwire)
342 for (
size_t i_t0 = 0; i_t0 < t0Handle->size(); i_t0++)
348 mf::LogVerbatim(
"RobustHitFinder") <<
"Found t0 at TPC tick = " << tick0;
350 std::vector<art::Ptr<raw::ExternalTrigger> > trigvec = triggers.at(i_t0);
352 std::vector<unsigned int> evtTriggers;
353 for (
auto const &trig : trigvec) evtTriggers.push_back(trig->GetTrigID());
382 std::cout <<
"AuxDet " <<
c1 <<
" not found. Aborting..." <<
std::endl;
396 std::cout <<
"AuxDet " <<
c2 <<
" not found. Aborting..." <<
std::endl;
401 int triggercode = 100*
c1+
c2;
438 for (
size_t iwire = 0; iwire < wireHandle->size(); ++iwire)
463 chan.
chanz =
static_cast<float>(wirexyz[2]);
478 chanMap.emplace(std::make_pair(chan.
channelID,chan));
481 if (chanMap.size() == 0)
487 for (
auto & chanmapitr : chanMap)
491 if (hitVec.size() == 0)
496 std::vector<dune::HitLineFitAlg::HitLineFitData> fitdata;
497 for (
auto &
hit : hitVec)
519 fitdata.push_back(hlfd);
525 for (
size_t i_h = 0; i_h < fitdata.size(); ++i_h)
527 hitVec[i_h].fitrealhit = fitdata[i_h].hitREAL;
532 if (retval == 1 || retval == 0)
540 for (
auto &
hit : hitVec)
565 hTracksByTrigger =
fTfs->make<TH1I>(
"hTracksByTrigger",
"Num reconstructed tracks by trigger",3716,0,3716);
568 fAllTree->Branch(
"run",&
run,
"run/I");
569 fAllTree->Branch(
"subrun",&
subrun,
"subrun/I");
570 fAllTree->Branch(
"event",&
event,
"event/I");
571 fAllTree->Branch(
"c1",&
c1,
"c1/i");
572 fAllTree->Branch(
"c2",&
c2,
"c2/i");
573 fAllTree->Branch(
"trignum",&
trignum,
"trignum/i");
583 fTree =
fTfs->make<TTree>(
"RobustHitFinder",
"RobustHitFinder");
587 fTree->Branch(
"t0",&
t0,
"t0/D");
588 fTree->Branch(
"c1",&
c1,
"c1/i");
589 fTree->Branch(
"c2",&
c2,
"c2/i");
671 fEfield = p.
get<std::vector<float> >(
"Efield");
691 for (
int i_tpc = 0; i_tpc<8; ++i_tpc)
693 float chanbeginz = 99999;
694 float chanendz = -99999;
695 int chanbeginid = 99999;
696 int chanendid = -99999;
697 for (
auto & chanitr : chanMap)
699 if (chanitr.second.tpcNum != i_tpc)
continue;
700 chanitr.second.nGoodHits = 0;
701 for (
auto &
hit : hitVec)
703 if (
hit.channelID == chanitr.first &&
hit.fitrealhit)
705 chanitr.second.goodHitStartTick =
hit.hitBeginTick;
706 chanitr.second.goodHitEndTick =
hit.hitEndTick;
707 chanitr.second.nGoodHits++;
710 if (chanitr.second.nGoodHits > 0 && chanitr.second.chanz < chanbeginz)
712 chanbeginz = chanitr.second.chanz;
713 chanbeginid = chanitr.first;
715 if (chanitr.second.nGoodHits > 0 && chanitr.second.chanz > chanendz)
717 chanendz = chanitr.second.chanz;
718 chanendid = chanitr.first;
722 if (chanendz < chanbeginz)
724 mf::LogError(
"RobustHitFinder") <<
"Problem: chanbeginz=" << chanbeginz <<
" is greater than chanendz=" << chanendz;
728 if (chanbeginid == 99999 || chanendid == -99999)
734 if (chanMap[chanbeginid].tpcNum % 2 != chanMap[chanendid].tpcNum % 2)
736 mf::LogError(
"RobustHitFinder") <<
"Track crosses APA. Don't know how to deal with this yet.";
740 std::vector<std::pair<float,int> > trackchans;
741 for (
auto const & chanitr : chanMap)
743 if (chanitr.second.tpcNum != i_tpc)
continue;
744 if (chanitr.second.chanz >= chanbeginz && chanitr.second.chanz <= chanendz && chanitr.second.tpcNum % 2 == chanMap[chanbeginid].tpcNum % 2)
746 trackchans.push_back(std::make_pair(chanitr.second.chanz,chanitr.first));
749 std::sort(trackchans.begin(),trackchans.end());
751 for (
size_t i_tc = 0; i_tc < trackchans.size(); ++i_tc)
762 chnearlow = chanMap[trackchans[i_tc-sub].second];
770 chnearhigh = chanMap[trackchans[i_tc+
add].second];
779 double clz = chnearlow.
chanz;
780 double chz = chnearhigh.
chanz;
782 int chstart = cls + (((chs - cls) / (chz - clz)) * (ch.
chanz - clz));
783 int chend = cle + (((che - cle) / (chz - clz)) * (ch.
chanz - clz));
791 if (chz-clz < 0)
continue;
793 ch.
pulse_ends.push_back(std::make_pair(chstart,chend));
864 float tempsegmentlength = 0.449;
869 double thetayz = TMath::ATan2(model->Eval(hit.
hitz+1)-model->Eval(hit.
hitz-1),2);
870 double tan2thetayz = TMath::Power(TMath::Tan(thetayz),2);
872 double projL = sqrt(1+tan2thetayz+y2z2);
873 tempsegmentlength *=
static_cast<float>(projL);
877 double thetayx = TMath::ATan2(2,model->Eval(hit.
hitx+1)-model->Eval(hit.
hitx-1));
878 double tan2thetayx = TMath::Power(TMath::Tan(thetayx),2);
880 double projL = sqrt(1+tan2thetayx*(1+y2x2));
881 tempsegmentlength *=
static_cast<float>(projL);
892 for (
size_t i_hit = 0; i_hit < chan.
pulse_ends.size(); i_hit++)
896 int begin_index,end_index;
901 float peak = std::distance(chan.
signalVec.begin(),std::max_element(bi,ei));
914 std::vector<float> pulse(beginitr,enditr);
917 std::vector<float> pulseFilter(fbeginitr,fenditr);
920 hit.
hitWidth = end_index-begin_index;
924 hit.
hitSigmaIntegral = TMath::Sqrt(pulse.size())*TMath::RMS(pulse.size(),pulse.data());
925 hit.
hitSigmaIntegralFilter = TMath::Sqrt(pulseFilter.size())*TMath::RMS(pulseFilter.size(),pulseFilter.data());
973 float tempsegmentlength = 0.449;
976 double thetayz = TMath::ATan2((hit.
hitx+1)-(hit.
hitx-1),2);
977 double tan2thetayz = TMath::Power(TMath::Tan(thetayz),2);
979 double projL = sqrt(1+tan2thetayz+y2z2);
980 tempsegmentlength *=
static_cast<float>(projL);
984 double thetayx = TMath::ATan2(2,(hit.
hitz+1)-(hit.
hitz-1));
985 double tan2thetayx = TMath::Power(TMath::Tan(thetayx),2);
987 double projL = sqrt(1+tan2thetayx*(1+y2x2));
988 tempsegmentlength *=
static_cast<float>(projL);
1009 hit.
artHit = temphit.move();
1016 std::vector<MCHit> MCHitVec;
1025 bool channelexists =
false;
1026 for (
size_t i_wire = 0; i_wire < wireHandle->size(); ++i_wire)
1029 if (pwire->
Channel() == sc->Channel())
1031 channelexists =
true;
1035 if (!channelexists)
continue;
1038 if (fCSP->
IsBad(sc->Channel()))
continue;
1043 auto const& tdcidemap = sc->TDCIDEMap();
1045 for (
auto const& tdcIt : tdcidemap)
1048 auto const& ideVec = tdcIt.second;
1051 unsigned short tdc = tdcIt.first;
1054 for (
auto const& ideIt : ideVec)
1059 if (
abs(ideIt.trackID) != 1)
continue;
1061 if (part->
Mother() != 0)
continue;
1064 mchit.
idevec.push_back(ideIt);
1067 if (mchit.
idevec.size() > 1 && mchit.
channel != sc->Channel())
1068 throw cet::exception(
"RobustHitFinder") <<
"Channel IDs don't match!";
1071 mchit.
channel = sc->Channel();
1073 throw cet::exception(
"RobustHitFinder") <<
"Pdg from previous IDE (" << mchit.
pdg <<
") doesn't equal this IDE Pdg (" << part->
PdgCode() <<
")";
1079 mchit.
tdcvec.push_back(tdc);
1085 if (mchit.
idevec.size()==0)
continue;
1086 MCHitVec.push_back(mchit);
1089 for (
auto mchit : MCHitVec)
1092 int mcchannel = mchit.channel;
1094 double meanSimTime = TMath::Mean(mchit.tdcvec.begin(),mchit.tdcvec.end(),mchit.chargevec.begin());
1095 double sigmaSimTime = TMath::RMS(mchit.tdcvec.begin(),mchit.tdcvec.end(),mchit.chargevec.begin());
1096 double totalSimCharge = std::accumulate(mchit.chargevec.begin(),mchit.chargevec.end(),0);
1097 totalSimCharge *= 1.602e-4;
1098 totalSimCharge *= 14.0*7.615;
1099 totalSimCharge *= 2.808;
1104 wire.push_back(chanMap.at(mcchannel).wireID);
1105 tpc.push_back(chanMap.at(mcchannel).tpcNum);
1109 baseline.push_back(chanMap.at(mcchannel).baseline);
1110 rms.push_back(chanMap.at(mcchannel).rms);
1112 rmsFilter.push_back(chanMap.at(mcchannel).rmsFilter);
1113 integral.push_back(totalSimCharge);
1126 width.push_back(sigmaSimTime);
1127 hitx.push_back(
TimeToX(meanSimTime,chanMap.at(mcchannel).tpcNum));
1129 hitz.push_back(chanMap.at(mcchannel).chanz);
1137 hitt.push_back(meanSimTime);
1145 float tempsegmentlength = 0.449;
1149 double tan2thetayz = TMath::Power(TMath::Tan(thetayz),2);
1151 double projL = sqrt(1+tan2thetayz+y2z2);
1152 tempsegmentlength *=
static_cast<float>(projL);
1157 double tan2thetayx = TMath::Power(TMath::Tan(thetayx),2);
1159 double projL = sqrt(1+tan2thetayx*(1+y2x2));
1160 tempsegmentlength *=
static_cast<float>(projL);
1168 c1arg=999; c2arg=999;
1169 int contains_111 = 0, contains_112 = 0, contains_113 = 0;
1170 int contains_Ntrigs = 0, contains_NU = 0, contains_NL = 0, contains_SU = 0, contains_SL = 0;
1171 int contains_EL = 0, contains_WU = 0, contains_TEL = 0;
1172 for (
size_t i_c = 0; i_c < evtTriggers.size(); i_c++)
1174 unsigned int trigID = evtTriggers[i_c];
1177 if (trigID <= 5) contains_SL++;
1178 if (trigID >= 6 && trigID <= 15) contains_EL++;
1179 if (trigID >= 16 && trigID <= 21) contains_NL++;
1180 if (trigID >= 22 && trigID <= 27) contains_NU++;
1181 if (trigID >= 28 && trigID <= 37) contains_WU++;
1182 if (trigID >= 38 && trigID <= 43) contains_SU++;
1183 if (trigID >= 44 && trigID <= 92) contains_TEL++;
1184 if (trigID == 111) contains_111++;
1185 if (trigID == 112) contains_112++;
1186 if (trigID == 113) contains_113++;
1189 if (contains_111 + contains_112 + contains_113 != 1)
return false;
1191 (contains_NU || contains_NL || contains_SU || contains_SL || contains_EL || contains_WU))
return false;
1192 if (contains_Ntrigs != 3)
return false;
1193 if (contains_111 && (contains_NU || contains_NL || contains_SU || contains_SL))
return false;
1194 if (contains_112 && (contains_EL || contains_WU || contains_SU || contains_NL))
return false;
1195 if (contains_113 && (contains_EL || contains_WU || contains_NU || contains_SL))
return false;
1196 if (contains_111 && (!contains_EL || !contains_WU))
return false;
1197 if (contains_112 && (!contains_NU || !contains_SL))
return false;
1198 if (contains_113 && (!contains_SU || !contains_NL))
return false;
1200 std::vector<unsigned int> counterIDs;
1202 for (
size_t i_c = 0; i_c < evtTriggers.size(); i_c++)
1204 unsigned int trigID = evtTriggers[i_c];
1205 if (trigID >= 44 && trigID <= 100)
continue;
1206 if (trigID >= 111 && trigID <= 113)
1208 trignumarg = trigID;
1211 counterIDs.push_back(trigID);
1213 if (counterIDs.size() != 2)
return false;
1214 if (trignumarg == 0)
return false;
1216 if (trignumarg == 112 || trignumarg == 113)
1220 c1arg = counterIDs[0];
1221 c2arg = counterIDs[1];
1225 c1arg = counterIDs[1];
1226 c2arg = counterIDs[0];
1229 else if (trignumarg == 111)
1233 c1arg = counterIDs[0];
1234 c2arg = counterIDs[1];
1238 c1arg = counterIDs[1];
1239 c2arg = counterIDs[0];
1242 if (c1arg == c2arg)
return false;
1243 if (c1arg == 999 || c2arg == 999)
return false;
1255 float ld = -1,
l1 = -1, l2 = -1, l3 = -1;
1259 if (thistpc == 1 || thistpc == 3 || thistpc == 5 || thistpc == 7)
1266 else if (thistpc == 0 || thistpc == 2 || thistpc == 4 || thistpc == 6)
1280 if (ld < 0 ||
l1 < 0 || l2 < 0 || l3 < 0)
return -99998;
1282 float t3max = l3/v3;
1283 float t2max = t3max+(l2/v2);
1284 float t1max = t2max+(
l1/v1);
1285 float tdmax = t1max+(ld/vd);
1287 if (0 <= thistime && thistime < t3max)
return v3*thistime;
1288 else if (t3max <= thistime && thistime < t2max)
return l3+v2*(thistime-t3max);
1289 else if (t2max <= thistime && thistime < t1max)
return l3+l2+v1*(thistime-t2max);
1290 else if (t1max <= thistime && thistime < tdmax)
return l3+l2+
l1+vd*(thistime-t1max);
1303 if (driftdistance < -89999)
return driftdistance;
1306 if (thistpc == 1 || thistpc == 3 || thistpc == 5 || thistpc == 7)
1310 else if (thistpc == 0 || thistpc == 2 || thistpc == 4 || thistpc == 6)
1324 return (((hitloc-trigloc1).Cross(hitloc-trigloc2)).Mag()/(trigloc2-trigloc1).Mag());
virtual bool IsBad(raw::ChannelID_t channel) const =0
Returns whether the specified channel is bad in the current run.
Pedestal provider class for DUNE.
code to link reconstructed objects back to the MC truth information
virtual float PedRms(raw::ChannelID_t ch) const =0
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
RobustHitFinder(fhicl::ParameterSet const &p)
std::vector< float > hiterrzhi
std::vector< float > hiterrxhi
void reconfigure(fhicl::ParameterSet const &p)
EventNumber_t event() const
void SetTreeVariablesMCTruth(const ChanMap_t &chanMap, art::Event &e)
CLHEP::HepRandomEngine & fEngine
void FilterWaveform(std::vector< float > wf, std::vector< float > &fwf)
std::vector< std::vector< float > > signal
std::vector< float > hiterrxlo
std::vector< float > peaktime
std::vector< float > pedmean
std::vector< int > numGoodHitsChan
std::vector< float > hiterrylo
const double & Time() const
std::vector< bool > ismctruth
float TimeToDisplacement(float thistime)
const simb::MCParticle * TrackIdToParticle_P(int id) const
std::vector< float > driftdist
void RobustRMSBase(std::vector< float > wf, float &bl, float &r)
std::vector< float > integralFilter
Declaration of signal hit object.
std::vector< float > hitz
std::vector< float > perpdist
double MinX() const
Returns the world x coordinate of the start of the box.
Coord add(Coord c1, Coord c2)
std::vector< int > channel
std::string fCounterT0ModuleLabel
bool fUseMeasuredCounterPositions
int FitLine(std::vector< HitLineFitData > &data, HitLineFitResults &bestfit)
std::vector< float > peaktick
std::vector< float > amplitudeFilter
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
Definition of basic raw digits.
Planes which measure Z direction.
std::vector< float > integral
double MaxX() const
Returns the world x coordinate of the end of the box.
std::vector< float > rmsFilter
bool ValidTrigger(std::vector< unsigned int > evtTriggers, unsigned int &c1arg, unsigned int &c2arg, unsigned int &trignumarg)
SigType_t SignalType(geo::PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
int fMissedBufferTicksLow
std::vector< float > baseline
std::vector< float > hiterrzlo
virtual double TPCG4Time2Tick(double g4time) const =0
Converts simulation time into a TPC electronics time tick.
std::vector< float > hity
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.
float TimeToDriftDist(float thistime, unsigned int thistpc)
void SetSeed(UInt_t seed)
std::vector< sim::IDE > idevec
std::map< int, float > bestValError
detinfo::DetectorClocks const * fClks
geo::View_t View() const
Returns the view the channel belongs to.
Class managing the creation of a new recob::Hit object.
std::vector< bool > countercut
Helper functions to create a hit.
std::vector< float > peaktickFilter
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
A class handling a collection of hits and its associations.
#define DEFINE_ART_MODULE(klass)
void produce(art::Event &e) override
std::string fWireModuleLabel
raw::ChannelID_t Channel() const
Returns the ID of the channel (or InvalidChannelID)
void FindHits(dune::ChannelInformation &chan)
float TimeToX(float thistime, unsigned int thistpc)
std::vector< float > pedrms
T get(std::string const &key) const
std::vector< std::vector< float > > signalFilter
void SetParameter(int i, double startValue, double minValue, double maxValue)
void SetTreeVariables(const dune::ChannelInformation &chan, const dune::HitInformation &hit, const dune::HitLineFitAlg::HitLineFitResults &fitresult)
SubRunNumber_t subRun() const
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.
Class providing information about the quality of channels.
std::vector< float > fEfield
detinfo::DetectorProperties const * fDetProp
std::vector< float > baselineFilter
std::vector< float > sigmaintegralFilter
std::vector< unsigned short > tdcvec
Definition of data types for geometry description.
void put_into(art::Event &)
Moves the data into an event.
std::vector< float > hitt
AuxDetGeo const & AuxDet(unsigned int const ad=0) const
Returns the specified auxiliary detector.
std::vector< float > peaktimeFilter
Detector simulation of raw signals on wires.
Conversion of times between different formats and references.
std::vector< float > Signal() const
Return a zero-padded full length vector filled with RoI signal.
std::vector< int > begintick
std::vector< const simb::MCParticle * > particlevec
Encapsulate the geometry of an auxiliary detector.
int fMissedBufferTicksHigh
ProducesCollector & producesCollector() noexcept
std::vector< int > endtick
std::map< int, float > bestVal
WireGeo const & Wire(geo::WireID const &wireid) const
Returns the specified wire.
std::vector< HitInformation > HitVec_t
art::ServiceHandle< geo::Geometry > fGeom
std::vector< bool > assumedhit
dune::RMSHitFinderAlg fHitFinderAlg
virtual double DriftVelocity(double efield=0., double temperature=0.) const =0
std::vector< float > amplitude
std::vector< float > hiterryhi
void FillHitInformation(dune::ChannelInformation &chan, dune::HitVec_t &hitVec, bool assumedHit)
object containing MC truth information necessary for making RawDigits and doing back tracking ...
Interface for experiment-specific channel quality info provider.
art::ServiceHandle< art::TFileService > fTfs
virtual float PedMean(raw::ChannelID_t ch) const =0
Retrieve pedestal information.
void SetHorizVertRanges(float hmin, float hmax, float vmin, float vmax)
const std::vector< art::Ptr< sim::SimChannel > > & SimChannels() const
std::map< int, ChannelInformation > ChanMap_t
Declaration of basic channel signal object.
void MakeCounterPositionMap(std::string CounterDir, std::string CounterFile, std::map< unsigned int, std::pair< TVector3, std::vector< TVector3 > > > &CounterPositionMap, double fExtendCountersX=0, double fExtendCountersY=0, double fExtendCountersZ=0)
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
virtual double TPCTick2TrigTime(double tick) const =0
Converts a TPC time (in ticks) into a trigger time [µs].
void MakeupMissedHits(ChanMap_t &chanMap, HitVec_t &hitVec)
std::vector< float > sumADC
float fHitGeomDistanceCut
void GetCenter(double *xyz, double localz=0.0) const
Fills the world coordinate of a point on the wire.
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Interface for experiment-specific service for channel quality info.
float hitGeomDist(TVector3 hitloc, TVector3 trigloc1, TVector3 trigloc2)
dune::HitLineFitAlg fFitAlg
std::vector< float > sigmaintegral
void SetSearchTicks(int s, int e)
Tools and modules for checking out the basics of the Monte Carlo.
std::vector< double > chargevec
std::map< unsigned int, std::pair< TVector3, std::vector< TVector3 > > > fCounterPositionMap
std::vector< float > hitx
const double * PlaneLocation(unsigned int p) const
std::vector< bool > fitrealhit
cet::coded_exception< error, detail::translate > exception
std::vector< float > segmentlength
QTextStream & endl(QTextStream &s)
Event finding and building.
void GetCenter(double *xyz, double localz=0.0) const
Return the center position of an AuxDet.
std::vector< int > allchannels
std::vector< int > signalsize
Signal from collection planes.
const lariov::DetPedestalProvider & fPedestalRetrievalAlg