LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
ShowerCalorimetry_module.cc
Go to the documentation of this file.
1 // Class: ShowerCalorimetry
3 // Plugin Type: producer (art v3_02_06)
4 // File: ShowerCalorimetry_module.cc
5 // Authors: calcuttj@msu.edu
6 // ahiguera@Central.UH.EDU
7 // tjyang@fnal.gov
8 // Generated at Fri Jul 12 14:14:46 2019 by Jacob Calcutt using cetskelgen
9 // from cetlib version v3_07_02.
10 //
12 
22 #include "fhiclcpp/ParameterSet.h"
24 
38 
39 #include <TVector3.h>
40 
41 #include <memory>
42 #include <utility>
43 #include <vector>
44 
45 namespace calo {
46  class ShowerCalorimetry;
47 }
48 
50 public:
51  explicit ShowerCalorimetry(fhicl::ParameterSet const& p);
52 
53 private:
54  void produce(art::Event& e) override;
55 
58  bool fSCE;
60 };
61 
63  : EDProducer{p}
64  , fShowerTag(p.get<art::InputTag>("ShowerTag"))
65  , fSpacePointTag(p.get<art::InputTag>("SpacePointTag"))
66  , fSCE(p.get<bool>("CorrectSCE"))
67  , caloAlg(p.get<fhicl::ParameterSet>("CalorimetryAlg"))
68 {
69  produces<std::vector<anab::Calorimetry>>();
70  produces<art::Assns<recob::Shower, anab::Calorimetry>>();
71 }
72 
74 {
75 
77 
78  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(e);
79  auto const detProp =
81  auto const* sce = lar::providerFrom<spacecharge::SpaceChargeService>();
82 
83  //Make the container for the calo product to put onto the event.
84  auto caloPtr = std::make_unique<std::vector<anab::Calorimetry>>();
85  auto& caloVector = *caloPtr;
86 
87  //Make a container for the track<-->calo associations.
88  //One entry per track, with entry equal to index in calorimetry collection of associated object.
89  std::vector<size_t> assnShowerCaloVector;
90  auto associationPtr = std::make_unique<art::Assns<recob::Shower, anab::Calorimetry>>();
91 
92  auto showerHandle = e.getValidHandle<std::vector<recob::Shower>>(fShowerTag);
93 
94  //Turn it into a vector of art pointers
95  std::vector<art::Ptr<recob::Shower>> recoShowers;
96  art::fill_ptr_vector(recoShowers, showerHandle);
97 
98  //Also get the hits from all the showers
99  art::FindManyP<recob::Hit> findHitsFromShowers(showerHandle, e, fShowerTag);
100 
101  //Go through all of the reconstructed showers in the event
102  for (size_t i = 0; i < recoShowers.size(); ++i) {
103  auto& shower = recoShowers[i];
104 
105  int shower_index = shower.key();
106  MF_LOG_INFO("ShowerCalorimetry") << "Getting Calorimetry info for " << shower_index << "\n";
107 
108  //This wil be used in the calorimetry object later
109  float shower_length = shower->Length();
110  //Get the hits from this shower
111  auto const& hits = findHitsFromShowers.at(shower_index);
112 
114 
115  //Sort the hits by their plane
116  //This vector stores the index of each hit on each plane
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);
120  }
121 
122  //Go through each plane and make calorimetry objects
123  for (size_t j = 0; j < geom->Nplanes(); ++j) {
124 
125  size_t hits_in_plane = hit_indices_per_plane[j].size();
126 
127  //Reserve vectors for each part of the calorimetry object
128  std::vector<float> dEdx(hits_in_plane);
129  std::vector<float> dQdx(hits_in_plane);
130  std::vector<float> pitch(hits_in_plane);
131 
132  //residual range, xyz, and deadwire default for now
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);
137 
138  geo::PlaneID planeID;
139 
140  float kineticEnergy = 0.;
141 
142  for (size_t k = 0; k < hits_in_plane; ++k) {
143 
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());
149 
150  float theHit_Xpos = -999.;
151  float theHit_Ypos = -999.;
152  float theHit_Zpos = -999.;
153 
154  auto const& sp = spFromShowerHits.at(hit_index);
155  if (!sp.empty()) {
156  //only use first space point
157  theHit_Xpos = sp[0]->XYZ()[0];
158  theHit_Ypos = sp[0]->XYZ()[1];
159  theHit_Zpos = sp[0]->XYZ()[2];
160  }
161  else {
162  MF_LOG_INFO("ShowerCalorimetry")
163  << "no sp associated w/this hit ... we will skip this hit";
164  continue;
165  }
166 
167  TVector3 pos(theHit_Xpos, theHit_Ypos, theHit_Zpos);
168  //correct pitch for hit direction
169  float this_pitch = wire_pitch;
170  float angleToVert = geom->WireAngleToVertical(theHit->View(), theHit->WireID());
171  angleToVert -= 0.5 * ::util::pi<>();
172 
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;
178 
179  //Correct for SCE
180  if (fSCE && sce->EnableCalSpatialSCE()) {
181 
182  geo::Vector_t posOffsets = {0., 0., 0.};
183  geo::Vector_t dirOffsets = {0., 0., 0.};
184 
185  posOffsets = sce->GetCalPosOffsets(geo::Point_t(pos), theHit->WireID().TPC);
186 
187  //For now, use the shower direction from Pandora...a better idea?
188  dirOffsets =
189  sce->GetCalPosOffsets(geo::Point_t{pos.X() + this_pitch * shower->Direction().X(),
190  pos.Y() + this_pitch * shower->Direction().Y(),
191  pos.Z() + this_pitch * shower->Direction().Z()},
192  theHit->WireID().TPC);
193 
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()};
198 
199  pitch[k] = dir_corr.Mag();
200  //correct position for SCE
201  theHit_Xpos -= posOffsets.X();
202  theHit_Ypos += posOffsets.Y();
203  theHit_Zpos += posOffsets.Z();
204  }
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());
209 
210  dQdx[k] = theHit->Integral() / pitch[k];
211 
212  // Just for now, use dQdx for dEdx
213  // dEdx[k] = theHit->Integral() / this_pitch;
214  dEdx[k] = caloAlg.dEdx_AREA(clockData, detProp, *theHit, pitch[k]),
215 
216  kineticEnergy += dEdx[k];
217  }
218 
219  //Make a calo object in the vector
220  caloVector.emplace_back(kineticEnergy,
221  dEdx,
222  dQdx,
223  resRange,
224  deadwires,
225  shower_length,
226  pitch,
228  hitIndex,
229  planeID);
230 
231  //Place the shower index in the association object
232  assnShowerCaloVector.emplace_back(shower_index);
233  }
234  }
235 
236  //Make the associations for ART
237  for (size_t i = 0; i < assnShowerCaloVector.size(); i++) {
238  if (assnShowerCaloVector[i] == std::numeric_limits<size_t>::max()) continue;
239 
240  art::Ptr<recob::Shower> shower_ptr(showerHandle, assnShowerCaloVector[i]);
241  util::CreateAssn(e, caloVector, shower_ptr, *associationPtr, i);
242  }
243 
244  //Finish up: Put the objects into the event
245  e.put(std::move(caloPtr));
246  e.put(std::move(associationPtr));
247 }
248 
Utilities related to art service access.
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:160
Declaration of signal hit object.
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
bool isValid
Whether this ID points to a valid element.
Definition: geo_types.h:210
constexpr auto abs(T v)
Returns the absolute value of the argument.
ShowerCalorimetry(fhicl::ParameterSet const &p)
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
void produce(art::Event &e) override
double WireAngleToVertical(View_t view, TPCID const &tpcid) const
Returns the angle of the wires in the specified view from vertical.
void hits()
Definition: readHits.C:15
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
Definition: TrackingTypes.h:68
#define MF_LOG_INFO(category)
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
Definition: PFPUtils.cxx:2675
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.
Definition: geo_vectors.h:180
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Utility object to perform functions of association.
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.
Definition: GeometryCore.h:977
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
Float_t e
Definition: plot.C:35
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
art framework interface to geometry description
calorimetry