81 auto const*
sce = lar::providerFrom<spacecharge::SpaceChargeService>();
84 auto caloPtr = std::make_unique<std::vector<anab::Calorimetry>>();
85 auto& caloVector = *caloPtr;
89 std::vector<size_t> assnShowerCaloVector;
90 auto associationPtr = std::make_unique<art::Assns<recob::Shower, anab::Calorimetry>>();
95 std::vector<art::Ptr<recob::Shower>> recoShowers;
102 for (
size_t i = 0; i < recoShowers.size(); ++i) {
103 auto&
shower = recoShowers[i];
105 int shower_index =
shower.key();
106 MF_LOG_INFO(
"ShowerCalorimetry") <<
"Getting Calorimetry info for " << shower_index <<
"\n";
109 float shower_length =
shower->Length();
111 auto const&
hits = findHitsFromShowers.at(shower_index);
117 std::vector<std::vector<size_t>> hit_indices_per_plane(geom->
Nplanes());
118 for (
size_t j = 0; j <
hits.size(); ++j) {
119 hit_indices_per_plane[
hits[j]->WireID().Plane].push_back(j);
123 for (
size_t j = 0; j < geom->
Nplanes(); ++j) {
125 size_t hits_in_plane = hit_indices_per_plane[j].size();
128 std::vector<float>
dEdx(hits_in_plane);
129 std::vector<float> dQdx(hits_in_plane);
130 std::vector<float> pitch(hits_in_plane);
133 std::vector<float> resRange(hits_in_plane, 0.);
134 std::vector<TVector3> xyz(hits_in_plane, TVector3(0., 0., 0.));
135 std::vector<float> deadwires(hits_in_plane, 0.);
136 std::vector<size_t> hitIndex(hits_in_plane);
140 float kineticEnergy = 0.;
142 for (
size_t k = 0; k < hits_in_plane; ++k) {
144 size_t hit_index = hit_indices_per_plane[j][k];
145 auto& theHit =
hits[hit_index];
146 if (!planeID.
isValid) { planeID = theHit->WireID(); }
147 hitIndex[k] = theHit.key();
148 float wire_pitch = geom->
WirePitch(theHit->View());
150 float theHit_Xpos = -999.;
151 float theHit_Ypos = -999.;
152 float theHit_Zpos = -999.;
154 auto const& sp = spFromShowerHits.at(hit_index);
157 theHit_Xpos = sp[0]->XYZ()[0];
158 theHit_Ypos = sp[0]->XYZ()[1];
159 theHit_Zpos = sp[0]->XYZ()[2];
163 <<
"no sp associated w/this hit ... we will skip this hit";
167 TVector3 pos(theHit_Xpos, theHit_Ypos, theHit_Zpos);
169 float this_pitch = wire_pitch;
171 angleToVert -= 0.5 * ::util::pi<>();
173 float cosgamma =
std::abs(sin(angleToVert) *
shower->Direction().Y() +
174 cos(angleToVert) *
shower->Direction().Z());
175 if (cosgamma > 0) this_pitch /= cosgamma;
176 if (this_pitch < wire_pitch) this_pitch = wire_pitch;
177 pitch[k] = this_pitch;
180 if (
fSCE &&
sce->EnableCalSpatialSCE()) {
185 posOffsets =
sce->GetCalPosOffsets(
geo::Point_t(pos), theHit->WireID().TPC);
190 pos.Y() + this_pitch *
shower->Direction().Y(),
191 pos.Z() + this_pitch *
shower->Direction().Z()},
192 theHit->WireID().TPC);
194 TVector3 dir_corr = {
195 this_pitch *
shower->Direction().X() - dirOffsets.X() + posOffsets.X(),
196 this_pitch *
shower->Direction().Y() + dirOffsets.Y() - posOffsets.Y(),
197 this_pitch *
shower->Direction().Z() + dirOffsets.Z() - posOffsets.Z()};
199 pitch[k] = dir_corr.Mag();
201 theHit_Xpos -= posOffsets.X();
202 theHit_Ypos += posOffsets.Y();
203 theHit_Zpos += posOffsets.Z();
205 xyz[k].SetXYZ(theHit_Xpos, theHit_Ypos, theHit_Zpos);
206 resRange[k] = std::hypot(theHit_Xpos -
shower->ShowerStart().X(),
207 theHit_Ypos -
shower->ShowerStart().Y(),
208 theHit_Zpos -
shower->ShowerStart().Z());
210 dQdx[k] = theHit->Integral() / pitch[k];
216 kineticEnergy +=
dEdx[k];
220 caloVector.emplace_back(kineticEnergy,
232 assnShowerCaloVector.emplace_back(shower_index);
237 for (
size_t i = 0; i < assnShowerCaloVector.size(); i++) {
238 if (assnShowerCaloVector[i] == std::numeric_limits<size_t>::max())
continue;
245 e.
put(std::move(caloPtr));
246 e.
put(std::move(associationPtr));
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
The data type to uniquely identify a Plane.
bool isValid
Whether this ID points to a valid element.
constexpr auto abs(T v)
Returns the absolute value of the argument.
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
double WireAngleToVertical(View_t view, TPCID const &tpcid) const
Returns the angle of the wires in the specified view from vertical.
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
#define MF_LOG_INFO(category)
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
bool CreateAssn(art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
double dEdx_AREA(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
art::InputTag fSpacePointTag
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
1.8.11 
