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