LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
ShowerCheater_module.cc
Go to the documentation of this file.
1 //
3 // ShowerCheater module
4 //
5 // brebel@fnal.gov
6 //
8 
9 #include <string>
10 
11 // LArSoft includes
21 
22 // Framework includes
29 #include "fhiclcpp/ParameterSet.h"
31 
32 namespace shwf {
33  class ShowerCheater : public art::EDProducer {
34  public:
35  explicit ShowerCheater(fhicl::ParameterSet const& pset);
36 
37  private:
38  void produce(art::Event& evt) override;
39 
40  std::string const fCheatedClusterLabel;
41  std::string const fG4ModuleLabel;
43  };
44 }
45 
46 namespace shwf {
47 
48  //--------------------------------------------------------------------
50  : EDProducer{pset}
51  , fCheatedClusterLabel{pset.get<std::string>("CheatedClusterLabel", "cluster")}
52  , fG4ModuleLabel{pset.get<std::string>("G4ModuleLabel", "largeant")}
53  {
54  produces<std::vector<recob::Shower>>();
55  produces<std::vector<recob::SpacePoint>>();
56  produces<art::Assns<recob::Shower, recob::Cluster>>();
57  produces<art::Assns<recob::Shower, recob::SpacePoint>>();
58  produces<art::Assns<recob::Shower, recob::Hit>>();
59  produces<art::Assns<recob::Hit, recob::SpacePoint>>();
60  }
61 
62  //--------------------------------------------------------------------
64  {
68 
69  // grab the clusters that have been reconstructed
71  evt.getByLabel(fCheatedClusterLabel, clustercol);
72 
74 
75  // make a vector of them - we aren't writing anything out to a file
76  // so no need for a art::PtrVector here
77  std::vector<art::Ptr<recob::Cluster>> clusters;
78  art::fill_ptr_vector(clusters, clustercol);
79 
80  // make a map of vectors of art::Ptrs keyed by eveID values
81  std::map<int, std::vector<std::pair<size_t, art::Ptr<recob::Cluster>>>> eveClusterMap;
82 
83  // loop over all clusters and fill in the map
84  for (size_t c = 0; c < clusters.size(); ++c) {
85 
86  // in the ClusterCheater module we set the cluster ID to be
87  // the eve particle track ID*1000 + plane*100 + tpc*10 + cryostat number. The
88  // floor function on the cluster ID / 1000 will give us
89  // the eve track ID
90  int eveID = floor(clusters[c]->ID() / 1000.);
91 
92  std::pair<size_t, art::Ptr<recob::Cluster>> indexPtr(c, clusters[c]);
93 
94  eveClusterMap[eveID].push_back(indexPtr);
95 
96  } // end loop over clusters
97 
98  // loop over the map and make prongs
99  std::unique_ptr<std::vector<recob::Shower>> showercol(new std::vector<recob::Shower>);
100  std::unique_ptr<std::vector<recob::SpacePoint>> spcol(new std::vector<recob::SpacePoint>);
101  std::unique_ptr<art::Assns<recob::Shower, recob::Cluster>> scassn(
103  std::unique_ptr<art::Assns<recob::Shower, recob::Hit>> shassn(
105  std::unique_ptr<art::Assns<recob::Shower, recob::SpacePoint>> sspassn(
107  std::unique_ptr<art::Assns<recob::Hit, recob::SpacePoint>> sphassn(
109 
110  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
111 
112  for (auto const& clusterMapItr : eveClusterMap) {
113 
114  // separate out the hits for each particle into the different views
115  std::vector<std::pair<size_t, art::Ptr<recob::Cluster>>> const& eveClusters =
116  clusterMapItr.second;
117 
118  size_t startSPIndx = spcol->size();
119 
120  std::vector<art::Ptr<recob::Cluster>> ptrvs;
121  std::vector<size_t> idxs;
122 
123  for (auto const& idxPtr : eveClusters) {
124  idxs.push_back(idxPtr.first);
125  ptrvs.push_back(idxPtr.second);
126 
127  // need to make the space points for this prong
128  // loop over the hits for this cluster and make
129  // a space point for each one
130  // set the SpacePoint ID to be the cluster ID*10000
131  // + the hit index in the cluster PtrVector of hits
132  std::vector<art::Ptr<recob::Hit>> const& hits = fmh.at(idxPtr.first);
133 
134  for (size_t h = 0; h < hits.size(); ++h) {
135  art::Ptr<recob::Hit> hit = hits[h];
136  // add up the charge from the hits on the collection plane
137  std::vector<double> xyz = bt_serv->HitToXYZ(clockData, hit);
138  double sperr[6] = {0.01, 0.01, 0.1, 0.001, 0.001, 0.001};
139 
140  // make the space point and set its ID and XYZ
141  // the errors and chi^2 are set to "good" values as we know the information perfectly
142  recob::SpacePoint sp(&xyz[0], sperr, 0.9, idxPtr.second->ID() * 10000 + h);
143  spcol->push_back(sp);
144 
145  // associate the space point to the hit
146  util::CreateAssn(evt, *spcol, hit, *sphassn);
147 
148  } // end loop over hits
149  } // end loop over pairs of index values and cluster Ptrs
150 
151  size_t endSPIndx = spcol->size();
152 
153  // is this prong electro-magnetic in nature or
154  // hadronic/muonic? EM --> shower, everything else is a track
155  if (std::abs(pi_serv->ParticleList()[clusterMapItr.first]->PdgCode()) == 11 ||
156  std::abs(pi_serv->ParticleList()[clusterMapItr.first]->PdgCode()) == 22 ||
157  std::abs(pi_serv->ParticleList()[clusterMapItr.first]->PdgCode()) == 111) {
158 
159  mf::LogInfo("ShowerCheater")
160  << "prong of " << clusterMapItr.first << " is a shower with pdg code "
161  << pi_serv->ParticleList()[clusterMapItr.first]->PdgCode();
162 
163  // get the direction cosine for the eve ID particle
164  // just use the same for both the start and end of the prong
165  const TLorentzVector initmom = pi_serv->ParticleList()[clusterMapItr.first]->Momentum();
166  TVector3 dcos(
167  initmom.Px() / initmom.Mag(), initmom.Py() / initmom.Mag(), initmom.Pz() / initmom.Mag());
168  TVector3 dcosErr(1.e-3, 1.e-3, 1.e-3);
169 
171  //double maxTransWidth[2] = { util::kBogusD, util::kBogusD };
172  //double distanceMaxWidth = util::kBogusD;
173 
174  // add a prong to the collection. Make the prong
175  // ID be the same as the track ID for the eve particle
176  recob::Shower s;
177  s.set_id(showercol->size());
178  s.set_direction(dcos);
179  s.set_direction_err(dcosErr);
180  showercol->push_back(s);
181  // associate the shower with its clusters
182  util::CreateAssn(evt, *showercol, ptrvs, *scassn);
183 
184  // get the hits associated with each cluster and associate those with the shower
185  for (size_t i = 0; i < idxs.size(); ++i) {
186  std::vector<art::Ptr<recob::Hit>> hits = fmh.at(i);
187  util::CreateAssn(evt, *showercol, hits, *shassn);
188  }
189 
190  // associate the shower with the space points
191  util::CreateAssn(evt, *showercol, *spcol, *sspassn, startSPIndx, endSPIndx);
192 
193  mf::LogInfo("ShowerCheater") << "adding shower: \n"
194  << showercol->back() << "\nto collection.";
195 
196  } // end if this is a shower
197  } // end loop over the map
198 
199  evt.put(std::move(showercol));
200  evt.put(std::move(spcol));
201  evt.put(std::move(scassn));
202  evt.put(std::move(shassn));
203  evt.put(std::move(sspassn));
204  evt.put(std::move(sphassn));
205  } // end produce
206 
207 } // end namespace
208 
void produce(art::Event &evt) override
void set_direction_err(const TVector3 &dir_e)
Definition: Shower.h:135
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Declaration of signal hit object.
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
constexpr auto abs(T v)
Returns the absolute value of the argument.
shower finding
ShowerCheater(fhicl::ParameterSet const &pset)
std::vector< double > HitToXYZ(detinfo::DetectorClocksData const &clockData, const recob::Hit &hit) const
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
void set_id(const int id)
Definition: Shower.h:127
void hits()
Definition: readHits.C:15
void set_direction(const TVector3 &dir)
Definition: Shower.h:134
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
Declaration of cluster object.
std::string const fCheatedClusterLabel
Detector simulation of raw signals on wires.
const sim::ParticleList & ParticleList() const
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.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Utility object to perform functions of association.
std::string const fG4ModuleLabel
label for module running G4 and making particles, etc
TCEvent evt
Definition: DataStructs.cxx:8
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
Float_t e
Definition: plot.C:35
Namespace collecting geometry-related classes utilities.
art framework interface to geometry description