LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
HitCheater_module.cc
Go to the documentation of this file.
1 // Class: HitCheater
3 // Module Type: producer
4 // File: HitCheater_module.cc
5 //
6 // Generated at Tue Nov 8 09:41:20 2011 by Brian Rebel using artmod
7 // from art v1_00_02.
9 #ifndef HitCheater_h
10 #define HitCheater_h
11 
12 // C/C++ standard libraries
13 #include <cmath> // std::sqrt()
14 #include <string>
15 #include <map>
16 #include <vector>
17 
25 #include "cetlib_except/exception.h"
26 
28 #include "larcoreobj/SimpleTypesAndConstants/RawTypes.h" // raw::ChannelID_t
40 
41 namespace hit {
42  class HitCheater;
43 }
44 
46 public:
47  explicit HitCheater(fhicl::ParameterSet const & p);
48  virtual ~HitCheater();
49 
50  virtual void produce(art::Event & e);
51 
52  virtual void beginJob();
53  virtual void reconfigure(fhicl::ParameterSet const & p);
54 
55 private:
56 
57  void FindHitsOnChannel(
58  const sim::SimChannel* sc,
59  std::vector<recob::Hit>& hits,
60  int spill
61  );
62 
63 
64  std::string fG4ModuleLabel;
65  std::string fWireModuleLabel;
66  double fMinCharge;
67  double fElectronsToADC;
71  double fSamplingRate;
74 };
75 
76 //-------------------------------------------------------------------
78 {
79  this->reconfigure(p);
80 
81  // let HitCollectionCreator declare that we are going to produce
82  // hits and associations with wires and raw digits
83  // (with no particular product label)
85 }
86 
87 //-------------------------------------------------------------------
89 {
90 }
91 
92 //-------------------------------------------------------------------
94 {
95  // this object contains the hit collection
96  // and its associations to wires and raw digits:
98 
99  // Read in the wire List object(s).
102  e.getValidHandle<std::vector<recob::Wire>>(WireInputTag);
103 
104  int whatSpill = 1;
105  if( !fCalDataProductInstanceName.empty() ) {
106  if( fCalDataProductInstanceName.find("ost") != std::string::npos) whatSpill=2;
107  else whatSpill=0;
108  }
109 
110  // also get the raw digits associated with the wires;
111  // we assume they have been created by the same module as the wires
112  art::FindOneP<raw::RawDigit> WireToRawDigits(wHandle, e, WireInputTag);
113 
114  // make a map of wires to channel numbers
115  std::map<raw::ChannelID_t, art::Ptr<recob::Wire>> wireMap;
116 
117  for(size_t wc = 0; wc < wHandle->size(); ++wc){
118  art::Ptr<recob::Wire> wire(wHandle, wc);
119  wireMap[wire->Channel()] = wire;
120  }
121 
122  // get the sim::SimChannels out of the event
123  std::vector<const sim::SimChannel*> sccol;
124  e.getView(fG4ModuleLabel, sccol);
125 
126  // find the hits on each channel
127  for(sim::SimChannel const* sc: sccol) {
128  std::vector<recob::Hit> hits_on_channel;
129 
130  FindHitsOnChannel(sc, hits_on_channel, whatSpill);
131 
132  art::Ptr<recob::Wire> const& wire = wireMap[sc->Channel()];
133  art::Ptr<raw::RawDigit> rawdigits; // null by default
134  if (wire.isNonnull()) rawdigits = WireToRawDigits.at(wire.key());
135 
136  // add all the hits found on this channel to the data product,
137  // all associated to the same hit and wire
138  for (recob::Hit& hit: hits_on_channel)
139  hits.emplace_back(std::move(hit), wire, rawdigits);
140 
141  }// end loop over SimChannels
142 
143  // put the cheated hits into the event
144  LOG_DEBUG("HitCheater") << "putting " << hits.size() << " hits into the event";
145  hits.put_into(e);
146 
147  return;
148 }
149 
150 //-------------------------------------------------------------------
152  std::vector<recob::Hit>& hits,
153  int spill)
154 {
156 
157  raw::ChannelID_t channel = sc->Channel();
158  geo::SigType_t signal_type = geo->SignalType(channel);
159  geo::View_t view = geo->View(channel);
160 
161 
162  // determine the possible geo::WireIDs for this particular channel
163  // then make a map of tdc to electrons for each one of those geo::WireIDs
164  // then find hits on each geo::WireID
165  std::vector<geo::WireID> wireids = geo->ChannelToWire(channel);
166 
167  std::map<geo::WireID, std::map< unsigned int, double> > wireIDSignals;
168 
169  auto const& idemap = sc->TDCIDEMap();
170 
171  for(auto const& mapitr : idemap){
172  unsigned short tdc = mapitr.first;
173 
175  if( tdc < spill*fReadOutWindowSize ||
176  tdc > (spill+1)*fReadOutWindowSize ) continue;
177  } else {
178  if ( tdc < 0 || tdc > fReadOutWindowSize) continue;
179  }
180 
181  // figure out which TPC we are in for each voxel
182 
183  for(auto const& ideitr : mapitr.second){
184 
185  const float edep = ideitr.numElectrons;
186 
187  std::array<double, 3> pos;
188  pos[0] = ideitr.x;
189  pos[1] = ideitr.y;
190  pos[2] = ideitr.z;
191 
192  geo::TPCID tpcID = geo->FindTPCAtPosition(pos.data());
193  if (!tpcID.isValid) {
194  mf::LogWarning("HitCheater") << "TPC for position ( "
195  << pos[0] << " ; " << pos[1] << " ; " << pos[2] << " )"
196  << " in no TPC; move on to the next sim::IDE";
197  continue;
198  }
199  const unsigned int tpc = tpcID.TPC;
200  const unsigned int cstat = tpcID.Cryostat;
201 
202  for( auto const& wid : wireids){
203  if(wid.TPC == tpc && wid.Cryostat == cstat){
204  // in the right TPC, now figure out which wire we want
205  // this works because there is only one plane option for
206  // each WireID in each TPC
207  if(wid.Wire == geo->NearestWire(pos.data(), wid.Plane, wid.TPC, wid.Cryostat))
208  wireIDSignals[wid][tdc] += edep;
209  }// end if in the correct TPC and Cryostat
210  }// end loop over wireids for this channel
211  }// end loop over ides for this channel
212  }// end loop over tdcs for this channel
213 
214  // now loop over each wire ID and determine where the hits are
215  for( auto const& widitr : wireIDSignals){
216 
217  // get the first tdc in the
218  unsigned short prev = widitr.second.begin()->first;
219  unsigned short startTime = prev;
220  double totCharge = 0.;
221  double maxCharge = -1.;
222  double peakTime = 0.;
223  int multiplicity = 1 ;
224  lar::util::StatCollector<double> time; // reduce rounding errors
225 
226  // loop over all the tdcs for this geo::WireID
227  for( auto tdcitr : widitr.second){
228  unsigned short tdc = tdcitr.first;
229  if (tdc < prev) {
231  << "SimChannel TDCs going backward!";
232  }
233 
234  // more than a one tdc gap between times with
235  // signal, start a new hit
236  if(tdc - prev > fNewHitTDCGap){
237 
238  if(totCharge > fMinCharge){
239  hits.emplace_back(
240  channel, // channel
241  (raw::TDCtick_t) startTime, // start_tick
242  (raw::TDCtick_t) prev, // end_tick
243  peakTime, // peak_time
244  1., // sigma_peak_time
245  time.RMS(), // RMS
246  maxCharge, // peak_amplitude
247  std::sqrt(maxCharge), // sigma_peak_amplitude
248  totCharge, // summedADC
249  totCharge, // hit_integral
250  std::sqrt(totCharge), // hit_sigma_integral
251  multiplicity, // multiplicity
252  0, // local_index
253  1., // goodness_of_fit
254  0., // dof
255  view, // view
256  signal_type, // signal_type
257  widitr.first // wireID
258  );
259 
260  LOG_DEBUG("HitCheater") << "new hit is " << hits.back();
261 
262  }// end if charge is large enough
263 
264  // reset the variables for each hit
265  startTime = tdc;
266  peakTime = tdc;
267  totCharge = 0.;
268  maxCharge = -1.;
269  time.clear();
270 
271  }// end if need to start a new hit
272 
273  const double adc = fElectronsToADC*tdcitr.second;
274 
275  totCharge += adc;
276  // use ADC as weight; the shift reduces the precision needed;
277  // average would need to be shifted back: time.Averatge() + startTime
278  // but RMS is not affected
279  time.add(tdc - startTime, adc);
280  if(adc > maxCharge){
281  maxCharge = adc;
282  peakTime = tdc;
283  }
284 
285  prev = tdc;
286 
287  }// end loop over tdc values for the current geo::WireID
288 
289 
290  // We might have missed the last hit, so do it now
291  if(totCharge > fMinCharge){
292  hits.emplace_back(
293  channel, // channel
294  (raw::TDCtick_t) startTime, // start_tick
295  (raw::TDCtick_t) prev + 1, // end_tick; prev included in the hit
296  peakTime, // peak_time
297  1., // sigma_peak_time
298  time.RMS(), // RMS
299  maxCharge, // peak_amplitude
300  std::sqrt(maxCharge), // sigma_peak_amplitude
301  totCharge, // summedADC
302  totCharge, // hit_integral
303  std::sqrt(totCharge), // hit_sigma_integral
304  multiplicity, // multiplicity
305  0, // local_index
306  1., // goodness_of_fit
307  0., // dof
308  view, // view
309  signal_type, // signal_type
310  widitr.first // wireID
311  );
312 
313  LOG_DEBUG("HitCheater") << "last hit is " << hits.back();
314 
315  }// end if charge is large enough
316 
317  }// end loop over map of geo::WireID to map<tdc,electrons>
318 
319 
320  return;
321 }
322 
323 //-------------------------------------------------------------------
325 {
326  return;
327 }
328 
329 //-------------------------------------------------------------------
331 {
332  fG4ModuleLabel = p.get< std::string >("G4ModuleLabel", "largeant");
333  fWireModuleLabel = p.get< std::string >("WireModuleLabel", "caldata" );
334  fMinCharge = p.get< double >("MinimumCharge", 5. );
335  fNewHitTDCGap = p.get< int >("NewHitTDCGap", 1 );
336 
337  const detinfo::DetectorProperties* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
338  fElectronsToADC = detprop->ElectronsToADC();
339  fSamplingRate = detprop->SamplingRate();
340  fTriggerOffset = detprop->TriggerOffset();
341 
343  size_t pos = fWireModuleLabel.find(":");
344  if( pos!=std::string::npos ) {
345  fCalDataProductInstanceName = fWireModuleLabel.substr( pos+1 );
346  fWireModuleLabel = fWireModuleLabel.substr( 0, pos );
347  }
348 
349  fReadOutWindowSize = detprop->ReadOutWindowSize();
350  fNumberTimeSamples = detprop->NumberTimeSamples();
351 
352  return;
353 }
354 
355 #endif /* HitCheater_h */
356 
key_type key() const
Definition: Ptr.h:356
Store parameters for running LArG4.
int fTriggerOffset
from detinfo::DetectorPropertiesService
double fSamplingRate
from detinfo::DetectorPropertiesService
std::string fWireModuleLabel
label name for module making recob::Wires
Energy deposited on a readout channel by simulated tracks.
Definition: SimChannel.h:143
bool isNonnull() const
Definition: Ptr.h:335
int fNumberTimeSamples
Number of total time samples (N*readoutwindowsize)
virtual void beginJob()
int fNewHitTDCGap
gap allowed in tdcs without charge before making a new hit
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
void FindHitsOnChannel(const sim::SimChannel *sc, std::vector< recob::Hit > &hits, int spill)
Declaration of signal hit object.
bool isValid
Whether this ID points to a valid element.
Definition: geo_types.h:129
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:130
Definition of basic raw digits.
SigType_t SignalType(geo::PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
std::string fCalDataProductInstanceName
label name for module making recob::Wires
static void declare_products(ModuleType &producer, std::string instance_name="", bool doWireAssns=true, bool doRawDigitAssns=true)
Declares the hit products we are going to fill.
Definition: HitCreator.h:1117
Classes gathering simple statistics.
Weight_t RMS() const
Returns the root mean square.
double fMinCharge
Minimum charge required to make a hit.
int TDCtick_t
Type representing a TDC tick.
Definition: RawTypes.h:24
geo::TPCID FindTPCAtPosition(double const worldLoc[3]) const
Returns the ID of the TPC at specified location.
std::size_t getView(std::string const &moduleLabel, std::string const &productInstanceName, std::vector< ELEMENT const * > &result) const
Definition: DataViewImpl.h:474
Helper functions to create a hit.
void hits()
Definition: readHits.C:15
A class handling a collection of hits and its associations.
Definition: HitCreator.h:513
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:42
std::string fG4ModuleLabel
label name for module making sim::SimChannels
raw::ChannelID_t Channel() const
Returns the ID of the channel (or InvalidChannelID)
Definition: Wire.h:164
enum geo::_plane_sigtype SigType_t
Enumerate the possible plane projections.
T get(std::string const &key) const
Definition: ParameterSet.h:231
geo::WireID::WireID_t NearestWire(geo::Point_t const &point, geo::PlaneID const &planeid) const
Returns the index of wire closest to position in the specified TPC.
void emplace_back(recob::Hit &&hit, art::Ptr< recob::Wire > const &wire=art::Ptr< recob::Wire >(), art::Ptr< raw::RawDigit > const &digits=art::Ptr< raw::RawDigit >())
Adds the specified hit to the data collection.
Definition: HitCreator.cxx:245
The data type to uniquely identify a TPC.
Definition: geo_types.h:195
virtual void reconfigure(fhicl::ParameterSet const &p)
Double_t edep
Definition: macro.C:13
View_t View(geo::PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
Definition of data types for geometry description.
void put_into(art::Event &)
Moves the data into an event.
Definition: HitCreator.h:663
Detector simulation of raw signals on wires.
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
raw::ChannelID_t Channel() const
Returns the readout channel this object describes.
Definition: SimChannel.h:332
size_t size() const
Returns the number of hits currently in the collection.
Definition: HitCreator.h:647
Encapsulate the construction of a single detector plane.
int fReadOutWindowSize
Number of samples in a readout window; NOT total samples.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
object containing MC truth information necessary for making RawDigits and doing back tracking ...
HitCheater(fhicl::ParameterSet const &p)
#define LOG_DEBUG(id)
Declaration of basic channel signal object.
virtual void produce(art::Event &e)
TDCIDEs_t const & TDCIDEMap() const
Returns all the deposited energy information as stored.
Definition: SimChannel.h:331
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:49
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:27
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
double fElectronsToADC
Conversion factor of electrons to ADC counts.
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:203
Float_t e
Definition: plot.C:34
void clear()
Clears all the statistics.
Namespace collecting geometry-related classes utilities.
Collects statistics on a single quantity (weighted)
Definition: fwd.h:25
art framework interface to geometry description
Encapsulate the construction of a single detector plane.
void add(Data_t value, Weight_t weight=Weight_t(1.0))
Adds one entry with specified value and weight.