80 auto const*
sce = lar::providerFrom<spacecharge::SpaceChargeService>();
83 auto caloPtr = std::make_unique<std::vector<anab::Calorimetry>>();
84 auto& caloVector = *caloPtr;
88 std::vector<size_t> assnShowerCaloVector;
89 auto associationPtr = std::make_unique<art::Assns<recob::Shower, anab::Calorimetry>>();
94 std::vector<art::Ptr<recob::Shower>> recoShowers;
101 for (
size_t i = 0; i < recoShowers.size(); ++i) {
102 auto&
shower = recoShowers[i];
104 int shower_index =
shower.key();
105 MF_LOG_INFO(
"ShowerCalorimetry") <<
"Getting Calorimetry info for " << shower_index <<
"\n";
108 float shower_length =
shower->Length();
110 auto const&
hits = findHitsFromShowers.at(shower_index);
116 std::vector<std::vector<size_t>> hit_indices_per_plane(wireReadoutGeom.Nplanes());
117 for (
size_t j = 0; j <
hits.size(); ++j) {
118 hit_indices_per_plane[
hits[j]->WireID().Plane].push_back(j);
122 for (
size_t j = 0; j < wireReadoutGeom.Nplanes(); ++j) {
124 size_t hits_in_plane = hit_indices_per_plane[j].size();
127 std::vector<float>
dEdx(hits_in_plane);
128 std::vector<float> dQdx(hits_in_plane);
129 std::vector<float> pitch(hits_in_plane);
132 std::vector<float> resRange(hits_in_plane, 0.);
133 std::vector<TVector3> xyz(hits_in_plane, TVector3(0., 0., 0.));
134 std::vector<float> deadwires(hits_in_plane, 0.);
135 std::vector<size_t> hitIndex(hits_in_plane);
139 float kineticEnergy = 0.;
141 for (
size_t k = 0; k < hits_in_plane; ++k) {
143 size_t hit_index = hit_indices_per_plane[j][k];
144 auto& theHit =
hits[hit_index];
145 if (!planeID.
isValid) { planeID = theHit->WireID(); }
146 hitIndex[k] = theHit.key();
147 float wire_pitch = wireReadoutGeom.
Plane({0, 0, theHit->View()}).WirePitch();
149 float theHit_Xpos = -999.;
150 float theHit_Ypos = -999.;
151 float theHit_Zpos = -999.;
153 auto const& sp = spFromShowerHits.at(hit_index);
156 theHit_Xpos = sp[0]->XYZ()[0];
157 theHit_Ypos = sp[0]->XYZ()[1];
158 theHit_Zpos = sp[0]->XYZ()[2];
162 <<
"no sp associated w/this hit ... we will skip this hit";
166 TVector3 pos(theHit_Xpos, theHit_Ypos, theHit_Zpos);
168 float this_pitch = wire_pitch;
169 float angleToVert = wireReadoutGeom.WireAngleToVertical(theHit->View(), theHit->WireID());
170 angleToVert -= 0.5 * ::util::pi<>();
172 float cosgamma =
std::abs(sin(angleToVert) *
shower->Direction().Y() +
173 cos(angleToVert) *
shower->Direction().Z());
174 if (cosgamma > 0) this_pitch /= cosgamma;
175 if (this_pitch < wire_pitch) this_pitch = wire_pitch;
176 pitch[k] = this_pitch;
179 if (
fSCE &&
sce->EnableCalSpatialSCE()) {
184 posOffsets =
sce->GetCalPosOffsets(
geo::Point_t(pos), theHit->WireID().TPC);
189 pos.Y() + this_pitch *
shower->Direction().Y(),
190 pos.Z() + this_pitch *
shower->Direction().Z()},
191 theHit->WireID().TPC);
193 TVector3 dir_corr = {
194 this_pitch *
shower->Direction().X() - dirOffsets.X() + posOffsets.X(),
195 this_pitch *
shower->Direction().Y() + dirOffsets.Y() - posOffsets.Y(),
196 this_pitch *
shower->Direction().Z() + dirOffsets.Z() - posOffsets.Z()};
198 pitch[k] = dir_corr.Mag();
200 theHit_Xpos -= posOffsets.X();
201 theHit_Ypos += posOffsets.Y();
202 theHit_Zpos += posOffsets.Z();
204 xyz[k].SetXYZ(theHit_Xpos, theHit_Ypos, theHit_Zpos);
205 resRange[k] = std::hypot(theHit_Xpos -
shower->ShowerStart().X(),
206 theHit_Ypos -
shower->ShowerStart().Y(),
207 theHit_Zpos -
shower->ShowerStart().Z());
209 dQdx[k] = theHit->Integral() / pitch[k];
215 kineticEnergy +=
dEdx[k];
219 caloVector.emplace_back(kineticEnergy,
231 assnShowerCaloVector.emplace_back(shower_index);
236 for (
size_t i = 0; i < assnShowerCaloVector.size(); i++) {
237 if (assnShowerCaloVector[i] == std::numeric_limits<size_t>::max())
continue;
244 e.
put(std::move(caloPtr));
245 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.
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
PlaneID_t Plane
Index of the plane within its TPC.
#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
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
art::InputTag fSpacePointTag
1.8.11 
