24 #include "art_root_io/TFileService.h" 43 #include "range/v3/view.hpp" 46 using ranges::views::transform;
135 ,
fMinX(pset.
get<
double>("MinX", -1.e10))
136 ,
fMaxX(pset.
get<
double>("MaxX", 1.e10))
137 ,
fMinY(pset.
get<
double>("MinY", -1.e10))
138 ,
fMaxY(pset.
get<
double>("MaxY", 1.e10))
139 ,
fMinZ(pset.
get<
double>("MinZ", -1.e10))
140 ,
fMaxZ(pset.
get<
double>("MaxZ", 1.e10))
179 mf::LogInfo(
"SpacePointAna") <<
"SpacePointAna configured with the following parameters:\n" 193 art::TFileDirectory
dir = tfs->mkdir(
"sptana",
"SpacePointAna histograms");
195 unsigned int nwiresU = 0, nwiresV = 0, nwiresW = 0;
202 unsigned int const nwires = pgeom.Nwires();
203 switch (pgeom.View()) {
204 case geo::kU: nwiresU = nwires;
break;
205 case geo::kV: nwiresV = nwires;
break;
206 case geo::kZ: nwiresW = nwires;
break;
213 fHDTUE = dir.make<TH1F>(
"MCDTUE",
"U-Drift Electrons Time Difference", 100, -5., 5.);
214 fHDTVE = dir.make<TH1F>(
"MCDTVE",
"V-Drift Electrons Time Difference", 100, -5., 5.);
215 fHDTWE = dir.make<TH1F>(
"MCDTWE",
"W-Drift Electrons Time Difference", 100, -5., 5.);
216 fHDTUPull = dir.make<TH1F>(
"MCDTUPull",
"U-Drift Electrons Time Pull", 100, -50., 50.);
217 fHDTVPull = dir.make<TH1F>(
"MCDTVPull",
"V-Drift Electrons Time Pull", 100, -50., 50.);
218 fHDTWPull = dir.make<TH1F>(
"MCDTWPull",
"W-Drift Electrons Time Pull", 100, -50., 50.);
221 fHDTUV = dir.make<TH1F>(
"DTUV",
"U-V time difference", 100, -20., 20.);
222 fHDTVW = dir.make<TH1F>(
"DTVW",
"V-W time difference", 100, -20., 20.);
223 fHDTWU = dir.make<TH1F>(
"DTWU",
"W-U time difference", 100, -20., 20.);
225 "DTUVU",
"U-V time difference vs. U", 100, 0., double(nwiresU), 100, -20., 20.);
227 "DTUVV",
"U-V time difference vs. V", 100, 0., double(nwiresV), 100, -20., 20.);
229 "DTVWV",
"V-W time difference vs. V", 100, 0., double(nwiresV), 100, -20., 20.);
231 "DTVWW",
"V-W time difference vs. W", 100, 0., double(nwiresW), 100, -20., 20.);
233 "DTWUW",
"W-U time difference vs. W", 100, 0., double(nwiresW), 100, -20., 20.);
235 "DTWUU",
"W-U time difference vs. U", 100, 0., double(nwiresU), 100, -20., 20.);
236 fHS = dir.make<TH1F>(
"DS",
"Spatial Separation", 100, -2., 2.);
238 fHchisq = dir.make<TH1F>(
"chisq",
"Chisquare", 100, 0., 20.);
240 fHx = dir.make<TH1F>(
"xpos",
"X Position", 100, 0., 2. * geom->
DetHalfWidth());
243 fHz = dir.make<TH1F>(
"zpos",
"Z Position", 100, 0., geom->
DetLength());
244 fHAmpU = dir.make<TH1F>(
"ampU",
"U Hit Amplitude", 50, 0., 50.);
245 fHAmpV = dir.make<TH1F>(
"ampV",
"V Hit Amplitude", 50, 0., 50.);
246 fHAmpW = dir.make<TH1F>(
"ampW",
"W Hit Amplitude", 50, 0., 50.);
247 fHAreaU = dir.make<TH1F>(
"areaU",
"U Hit Area", 100, 0., 500.);
248 fHAreaV = dir.make<TH1F>(
"areaV",
"V Hit Area", 100, 0., 500.);
249 fHAreaW = dir.make<TH1F>(
"areaW",
"W Hit Area", 100, 0., 500.);
250 fHSumU = dir.make<TH1F>(
"sumU",
"U Hit Sum ADC", 100, 0., 500.);
251 fHSumV = dir.make<TH1F>(
"sumV",
"V Hit Sum ADC", 100, 0., 500.);
252 fHSumW = dir.make<TH1F>(
"sumW",
"W Hit Sum ADC", 100, 0., 500.);
254 fHMCdx = dir.make<TH1F>(
"MCdx",
"X MC Residual", 100, -1., 1.);
255 fHMCdy = dir.make<TH1F>(
"MCdy",
"Y MC Residual", 100, -1., 1.);
256 fHMCdz = dir.make<TH1F>(
"MCdz",
"Z MC Residual", 100, -1., 1.);
257 fHMCxpull = dir.make<TH1F>(
"MCxpull",
"X MC Pull", 100, -50., 50.);
258 fHMCypull = dir.make<TH1F>(
"MCypull",
"Y MC Pull", 100, -50., 50.);
259 fHMCzpull = dir.make<TH1F>(
"MCzpull",
"Z MC Pull", 100, -50., 50.);
292 int nclus = clusterh->size();
294 for (
int i = 0; i < nclus; ++i) {
296 std::vector<art::Ptr<recob::Hit>> clushits = fm.at(i);
297 int nhits = clushits.size();
301 ihit != clushits.end();
315 int nhits = hith->size();
318 for (
int i = 0; i < nhits; ++i)
335 for (
auto const& hitPtr : hits) {
344 std::vector<double> hitxyz = bt_serv->
HitToXYZ(clockData, hit);
345 tav = detProp.ConvertXToTicks(
353 fHDTUE->Fill(tpeak - tav);
357 fHDTVE->Fill(tpeak - tav);
361 fHDTWE->Fill(tpeak - tav);
370 std::vector<recob::SpacePoint> spts1;
371 std::vector<recob::SpacePoint> spts2;
372 std::vector<recob::SpacePoint> spts3;
386 <<
"Found " << spts1.size() <<
" space points using special time cut.";
401 <<
"Found " << spts2.size() <<
" space points using special seperation cut.";
413 MF_LOG_DEBUG(
"SpacePointAna") <<
"Found " << spts3.size()
414 <<
" space points using default cuts.";
420 for (
auto const& spt : spts1) {
421 if (spt.XYZ()[0] <
fMinX || spt.XYZ()[0] >
fMaxX || spt.XYZ()[1] <
fMinY ||
431 for (
auto const& hit1 : spthits | transform(
to_element)) {
433 unsigned int tpc1 = hit1WireID.
TPC;
434 unsigned int plane1 = hit1WireID.
Plane;
435 unsigned int wire1 = hit1WireID.
Wire;
440 for (
auto const& hit2 : spthits | transform(
to_element)) {
442 unsigned int tpc2 = hit2WireID.
TPC;
443 unsigned int plane2 = hit2WireID.
Plane;
444 unsigned int wire2 = hit2WireID.
Wire;
448 if (tpc1 == tpc2 && plane1 != plane2) {
456 fHDTUVU->Fill(
double(wire1), t1 - t2);
457 fHDTUVV->Fill(
double(wire2), t1 - t2);
461 fHDTWUW->Fill(
double(wire2), t2 - t1);
462 fHDTWUU->Fill(
double(wire1), t2 - t1);
468 fHDTVWV->Fill(
double(wire1), t1 - t2);
469 fHDTVWW->Fill(
double(wire2), t1 - t2);
473 fHDTUVU->Fill(
double(wire2), t2 - t1);
474 fHDTUVV->Fill(
double(wire1), t2 - t1);
480 fHDTWUW->Fill(
double(wire1), t1 - t2);
481 fHDTWUU->Fill(
double(wire2), t1 - t2);
485 fHDTVWV->Fill(
double(wire2), t2 - t1);
486 fHDTVWW->Fill(
double(wire1), t2 - t1);
496 for (
auto const& spt : spts2) {
497 if (spt.XYZ()[0] <
fMinX || spt.XYZ()[0] >
fMaxX || spt.XYZ()[1] <
fMinY ||
514 for (
auto const& spt : spts3) {
515 if (spt.XYZ()[0] <
fMinX || spt.XYZ()[0] >
fMaxX || spt.XYZ()[1] <
fMinY ||
520 fHx->Fill(spt.XYZ()[0]);
521 fHy->Fill(spt.XYZ()[1]);
522 fHz->Fill(spt.XYZ()[2]);
526 std::vector<art::Ptr<recob::Hit>> spthits;
528 for (
auto const& ptr : av_spthits) {
529 spthits.push_back(ptr);
534 for (
auto const& hitPtr : spthits) {
562 fHMCdx->Fill(spt.XYZ()[0] - mcxyz[0]);
563 fHMCdy->Fill(spt.XYZ()[1] - mcxyz[1]);
564 fHMCdz->Fill(spt.XYZ()[2] - mcxyz[2]);
565 if (spt.ErrXYZ()[0] > 0.)
566 fHMCxpull->Fill((spt.XYZ()[0] - mcxyz[0]) / std::sqrt(spt.ErrXYZ()[0]));
567 if (spt.ErrXYZ()[2] > 0.)
568 fHMCypull->Fill((spt.XYZ()[1] - mcxyz[1]) / std::sqrt(spt.ErrXYZ()[2]));
569 if (spt.ErrXYZ()[5] > 0.)
570 fHMCzpull->Fill((spt.XYZ()[2] - mcxyz[2]) / std::sqrt(spt.ErrXYZ()[5]));
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
void reserve(size_type n)
constexpr to_element_t to_element
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Declaration of signal hit object.
Length_t DetHalfWidth(TPCID const &tpcid=tpc_zero) const
Returns the half width of the active volume of the specified TPC.
CryostatID_t Cryostat
Index of cryostat.
Planes which measure Z direction.
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
geo::View_t View() const
View for the plane of the hit.
WireID_t Wire
Index of the wire within its plane.
geo::WireID const & WireID() const
Initial tdc tick for hit.
std::vector< double > HitToXYZ(detinfo::DetectorClocksData const &clockData, const recob::Hit &hit) const
EDAnalyzer(fhicl::ParameterSet const &pset)
Length_t DetLength(TPCID const &tpcid=tpc_zero) const
Returns the length of the active volume of the specified TPC.
const art::PtrVector< recob::Hit > & getAssociatedHits(const recob::SpacePoint &spt) const
const SpacePointAlg fSptalgSep
float PeakAmplitude() const
The estimated amplitude of the hit at its peak, in ADC units.
void makeMCTruthSpacePoints(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts) const
double separation(const art::PtrVector< recob::Hit > &hits) const
bool isValid() const noexcept
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
#define DEFINE_ART_MODULE(klass)
std::string fHitModuleLabel
void push_back(Ptr< U > const &p)
void makeSpacePoints(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts) const
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
void analyze(const art::Event &evt) override
bool merge() const noexcept
PlaneID_t Plane
Index of the plane within its TPC.
double correctedTime(detinfo::DetectorPropertiesData const &detProp, const recob::Hit &hit) const
Declaration of cluster object.
Detector simulation of raw signals on wires.
float PeakTime() const
Time of the signal peak, in tick units.
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
SpacePointAna(fhicl::ParameterSet const &pset)
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Encapsulate the construction of a single detector plane.
const SpacePointAlg fSptalgTime
std::vector< double > SpacePointHitsToWeightedXYZ(detinfo::DetectorClocksData const &clockData, std::vector< art::Ptr< recob::Hit >> const &hits) const
void bookHistograms(bool mc)
float SummedADC() const
The sum of calibrated ADC counts of the hit (0. by default)
const SpacePointAlg fSptalgDefault
constexpr WireID()=default
Default constructor: an invalid TPC ID.
float SigmaPeakTime() const
Uncertainty for the signal peak, in tick units.
Length_t DetHalfHeight(TPCID const &tpcid=tpc_zero) const
Returns the half height of the active volume of the specified TPC.
2D representation of charge deposited in the TDC/wire plane
TPCID_t TPC
Index of the TPC within its cryostat.
Algorithm for generating space points from hits.
std::string fClusterModuleLabel
art framework interface to geometry description
cet::coded_exception< error, detail::translate > exception
Encapsulate the construction of a single detector plane.