77 auto const*
sce = lar::providerFrom<spacecharge::SpaceChargeService>();
80 auto caloPtr = std::make_unique<std::vector<anab::Calorimetry>>();
81 auto& caloVector = *caloPtr;
85 std::vector<size_t> assnShowerCaloVector;
86 auto associationPtr = std::make_unique<art::Assns<recob::Shower, anab::Calorimetry>>();
91 std::vector<art::Ptr<recob::Shower>> recoShowers;
95 art::FindManyP<recob::Hit> findHitsFromShowers(showerHandle, e,
fShowerTag);
98 for (
size_t i = 0; i < recoShowers.size(); ++i) {
99 auto&
shower = recoShowers[i];
101 int shower_index =
shower.key();
102 MF_LOG_INFO(
"ShowerCalorimetry") <<
"Getting Calorimetry info for " << shower_index <<
"\n";
105 float shower_length =
shower->Length();
107 auto const& hits = findHitsFromShowers.at(shower_index);
109 art::FindManyP<recob::SpacePoint> spFromShowerHits(hits, e,
fSpacePointTag);
113 std::vector<std::vector<size_t>> hit_indices_per_plane(geom->
Nplanes());
114 for (
size_t j = 0; j < hits.size(); ++j) {
115 hit_indices_per_plane[hits[j]->WireID().Plane].push_back(j);
119 for (
size_t j = 0; j < geom->
Nplanes(); ++j) {
121 size_t hits_in_plane = hit_indices_per_plane[j].size();
124 std::vector<float>
dEdx(hits_in_plane);
125 std::vector<float> dQdx(hits_in_plane);
126 std::vector<float> pitch(hits_in_plane);
129 std::vector<float> resRange(hits_in_plane, 0.);
130 std::vector<TVector3> xyz(hits_in_plane, TVector3(0., 0., 0.));
131 std::vector<float> deadwires(hits_in_plane, 0.);
132 std::vector<size_t> hitIndex(hits_in_plane);
136 float kineticEnergy = 0.;
138 for (
size_t k = 0;
k < hits_in_plane; ++
k) {
140 size_t hit_index = hit_indices_per_plane[j][
k];
141 auto& theHit = hits[hit_index];
142 if (!planeID.
isValid) { planeID = theHit->WireID(); }
143 hitIndex[
k] = theHit.key();
144 float wire_pitch = geom->
WirePitch(theHit->View());
146 float theHit_Xpos = -999.;
147 float theHit_Ypos = -999.;
148 float theHit_Zpos = -999.;
150 auto const& sp = spFromShowerHits.at(hit_index);
153 theHit_Xpos = sp[0]->XYZ()[0];
154 theHit_Ypos = sp[0]->XYZ()[1];
155 theHit_Zpos = sp[0]->XYZ()[2];
159 <<
"no sp associated w/this hit ... we will skip this hit";
163 TVector3
pos(theHit_Xpos, theHit_Ypos, theHit_Zpos);
165 float this_pitch = wire_pitch;
167 angleToVert -= 0.5 * ::util::pi<>();
169 float cosgamma =
std::abs(sin(angleToVert) *
shower->Direction().Y() +
170 cos(angleToVert) *
shower->Direction().Z());
171 if (cosgamma > 0) this_pitch /= cosgamma;
172 if (this_pitch < wire_pitch) this_pitch = wire_pitch;
173 pitch[
k] = this_pitch;
176 if (
fSCE &&
sce->EnableCalSpatialSCE()) {
186 pos.Y() + this_pitch *
shower->Direction().Y(),
187 pos.Z() + this_pitch *
shower->Direction().Z()},
188 theHit->WireID().TPC);
190 TVector3 dir_corr = {
191 this_pitch *
shower->Direction().X() - dirOffsets.X() + posOffsets.X(),
192 this_pitch *
shower->Direction().Y() + dirOffsets.Y() - posOffsets.Y(),
193 this_pitch *
shower->Direction().Z() + dirOffsets.Z() - posOffsets.Z()};
195 pitch[
k] = dir_corr.Mag();
197 theHit_Xpos -= posOffsets.X();
198 theHit_Ypos += posOffsets.Y();
199 theHit_Zpos += posOffsets.Z();
201 xyz[
k].SetXYZ(theHit_Xpos, theHit_Ypos, theHit_Zpos);
203 theHit_Ypos -
shower->ShowerStart().Y(),
204 theHit_Zpos -
shower->ShowerStart().Z());
206 dQdx[
k] = theHit->Integral() / pitch[
k];
212 kineticEnergy +=
dEdx[
k];
216 caloVector.emplace_back(kineticEnergy,
228 assnShowerCaloVector.emplace_back(shower_index);
233 for (
size_t i = 0; i < assnShowerCaloVector.size(); i++) {
The data type to uniquely identify a Plane.
bool isValid
Whether this ID points to a valid element.
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
std::enable_if_t< std::is_arithmetic_v< T >, T > hypot(T x, T y)
double dEdx(float dqdx, float Efield)
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
bool CreateAssn(PRODUCER const &prod, art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t indx=UINT_MAX)
Creates a single one-to-one association.
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
static int max(int a, int b)
#define MF_LOG_INFO(category)
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
double dEdx_AREA(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
art::InputTag fSpacePointTag
double WireAngleToVertical(geo::View_t view, geo::TPCID const &tpcid) const
Returns the angle of the wires in the specified view from vertical.