ShowerCalorimetry_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // 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 //
11 ////////////////////////////////////////////////////////////////////////
12 
19 #include "canvas/Persistency/Common/FindManyP.h"
21 #include "fhiclcpp/ParameterSet.h"
23 
36 
37 #include <TVector3.h>
38 #include <memory>
39 
40 namespace calo {
41  class ShowerCalorimetry;
42 }
43 
45 public:
46  explicit ShowerCalorimetry(fhicl::ParameterSet const& p);
47 
48 private:
49  void produce(art::Event& e) override;
50 
53  bool fSCE;
55 };
56 
58  : EDProducer{p}
59  , fShowerTag(p.get<art::InputTag>("ShowerTag"))
60  , fSpacePointTag(p.get<art::InputTag>("SpacePointTag"))
61  , fSCE(p.get<bool>("CorrectSCE"))
62  , caloAlg(p.get<fhicl::ParameterSet>("CalorimetryAlg"))
63 {
64  produces<std::vector<anab::Calorimetry>>();
65  produces<art::Assns<recob::Shower, anab::Calorimetry>>();
66 }
67 
68 void
70 {
71 
73 
74  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(e);
75  auto const detProp =
77  auto const* sce = lar::providerFrom<spacecharge::SpaceChargeService>();
78 
79  //Make the container for the calo product to put onto the event.
80  auto caloPtr = std::make_unique<std::vector<anab::Calorimetry>>();
81  auto& caloVector = *caloPtr;
82 
83  //Make a container for the track<-->calo associations.
84  //One entry per track, with entry equal to index in calorimetry collection of associated object.
85  std::vector<size_t> assnShowerCaloVector;
86  auto associationPtr = std::make_unique<art::Assns<recob::Shower, anab::Calorimetry>>();
87 
88  auto showerHandle = e.getValidHandle<std::vector<recob::Shower>>(fShowerTag);
89 
90  //Turn it into a vector of art pointers
91  std::vector<art::Ptr<recob::Shower>> recoShowers;
92  art::fill_ptr_vector(recoShowers, showerHandle);
93 
94  //Also get the hits from all the showers
95  art::FindManyP<recob::Hit> findHitsFromShowers(showerHandle, e, fShowerTag);
96 
97  //Go through all of the reconstructed showers in the event
98  for (size_t i = 0; i < recoShowers.size(); ++i) {
99  auto& shower = recoShowers[i];
100 
101  int shower_index = shower.key();
102  MF_LOG_INFO("ShowerCalorimetry") << "Getting Calorimetry info for " << shower_index << "\n";
103 
104  //This wil be used in the calorimetry object later
105  float shower_length = shower->Length();
106  //Get the hits from this shower
107  auto const& hits = findHitsFromShowers.at(shower_index);
108 
109  art::FindManyP<recob::SpacePoint> spFromShowerHits(hits, e, fSpacePointTag);
110 
111  //Sort the hits by their plane
112  //This vector stores the index of each hit on each plane
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);
116  }
117 
118  //Go through each plane and make calorimetry objects
119  for (size_t j = 0; j < geom->Nplanes(); ++j) {
120 
121  size_t hits_in_plane = hit_indices_per_plane[j].size();
122 
123  //Reserve vectors for each part of the calorimetry object
124  std::vector<float> dEdx(hits_in_plane);
125  std::vector<float> dQdx(hits_in_plane);
126  std::vector<float> pitch(hits_in_plane);
127 
128  //residual range, xyz, and deadwire default for now
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);
133 
134  geo::PlaneID planeID;
135 
136  float kineticEnergy = 0.;
137 
138  for (size_t k = 0; k < hits_in_plane; ++k) {
139 
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());
145 
146  float theHit_Xpos = -999.;
147  float theHit_Ypos = -999.;
148  float theHit_Zpos = -999.;
149 
150  auto const& sp = spFromShowerHits.at(hit_index);
151  if (!sp.empty()) {
152  //only use first space point
153  theHit_Xpos = sp[0]->XYZ()[0];
154  theHit_Ypos = sp[0]->XYZ()[1];
155  theHit_Zpos = sp[0]->XYZ()[2];
156  }
157  else {
158  MF_LOG_INFO("ShowerCalorimetry")
159  << "no sp associated w/this hit ... we will skip this hit";
160  continue;
161  }
162 
163  TVector3 pos(theHit_Xpos, theHit_Ypos, theHit_Zpos);
164  //correct pitch for hit direction
165  float this_pitch = wire_pitch;
166  float angleToVert = geom->WireAngleToVertical(theHit->View(), theHit->WireID());
167  angleToVert -= 0.5 * ::util::pi<>();
168 
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;
174 
175  //Correct for SCE
176  if (fSCE && sce->EnableCalSpatialSCE()) {
177 
178  geo::Vector_t posOffsets = {0., 0., 0.};
179  geo::Vector_t dirOffsets = {0., 0., 0.};
180 
181  posOffsets = sce->GetCalPosOffsets(geo::Point_t(pos), theHit->WireID().TPC);
182 
183  //For now, use the shower direction from Pandora...a better idea?
184  dirOffsets =
185  sce->GetCalPosOffsets(geo::Point_t{pos.X() + this_pitch * shower->Direction().X(),
186  pos.Y() + this_pitch * shower->Direction().Y(),
187  pos.Z() + this_pitch * shower->Direction().Z()},
188  theHit->WireID().TPC);
189 
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()};
194 
195  pitch[k] = dir_corr.Mag();
196  //correct position for SCE
197  theHit_Xpos -= posOffsets.X();
198  theHit_Ypos += posOffsets.Y();
199  theHit_Zpos += posOffsets.Z();
200  }
201  xyz[k].SetXYZ(theHit_Xpos, theHit_Ypos, theHit_Zpos);
202  resRange[k] = std::hypot(theHit_Xpos - shower->ShowerStart().X(),
203  theHit_Ypos - shower->ShowerStart().Y(),
204  theHit_Zpos - shower->ShowerStart().Z());
205 
206  dQdx[k] = theHit->Integral() / pitch[k];
207 
208  // Just for now, use dQdx for dEdx
209  // dEdx[k] = theHit->Integral() / this_pitch;
210  dEdx[k] = caloAlg.dEdx_AREA(clockData, detProp, *theHit, pitch[k]),
211 
212  kineticEnergy += dEdx[k];
213  }
214 
215  //Make a calo object in the vector
216  caloVector.emplace_back(kineticEnergy,
217  dEdx,
218  dQdx,
219  resRange,
220  deadwires,
221  shower_length,
222  pitch,
224  hitIndex,
225  planeID);
226 
227  //Place the shower index in the association object
228  assnShowerCaloVector.emplace_back(shower_index);
229  }
230  }
231 
232  //Make the associations for ART
233  for (size_t i = 0; i < assnShowerCaloVector.size(); i++) {
234  if (assnShowerCaloVector[i] == std::numeric_limits<size_t>::max()) continue;
235 
236  art::Ptr<recob::Shower> shower_ptr(showerHandle, assnShowerCaloVector[i]);
237  util::CreateAssn(e, caloVector, shower_ptr, *associationPtr, i);
238  }
239 
240  //Finish up: Put the objects into the event
241  e.put(std::move(caloPtr));
242  e.put(std::move(associationPtr));
243 }
244 
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
bool isValid
Whether this ID points to a valid element.
Definition: geo_types.h:211
ShowerCalorimetry(fhicl::ParameterSet const &p)
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:164
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
art framework interface to geometry description
std::enable_if_t< std::is_arithmetic_v< T >, T > hypot(T x, T y)
Definition: hypot.h:60
void produce(art::Event &e) override
T abs(T value)
double dEdx(float dqdx, float Efield)
Definition: doAna.cpp:21
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
const double e
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
def move(depos, offset)
Definition: depos.py:107
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:441
p
Definition: test.py:223
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: DataViewImpl.h:686
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)
Definition: TrackingTypes.h:55
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.
Definition: geo_vectors.h:184
Declaration of signal hit object.
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)
Definition: Ptr.h:297
double WireAngleToVertical(geo::View_t view, geo::TPCID const &tpcid) const
Returns the angle of the wires in the specified view from vertical.
calorimetry