LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
MCHitFinder_module.cc
Go to the documentation of this file.
1 // Class: MCHitFinder
3 // Module Type: producer
4 // File: MCHitFinder_module.cc
5 //
6 // Generated at Tue Jun 24 07:30:08 2014 by Kazuhiro Terao using artmod
7 // from cetpkgsupport v1_05_04.
9 
16 #include "fhiclcpp/ParameterSet.h"
17 
18 #include <memory>
19 
21 
26 
27 namespace hit {
28 
29  class MCHitFinder;
30 
31  class MCHitFinder : public art::EDProducer {
32  public:
33  explicit MCHitFinder(fhicl::ParameterSet const& p);
34 
35  private:
36  void produce(art::Event& e) override;
37 
38  std::string fLArG4ModuleName;
39  bool fVerbose;
41  };
42 
44  {
45  fLArG4ModuleName = p.get<std::string>("LArG4ModuleName");
46  fMakeMCWire = p.get<bool>("MakeMCWire");
47  fVerbose = p.get<bool>("Verbose");
48  produces<std::vector<sim::MCHitCollection>>();
49  if (fMakeMCWire) produces<std::vector<sim::MCWireCollection>>();
50  }
51 
53  {
54  const unsigned int nch = art::ServiceHandle<geo::WireReadout const>()->Get().Nchannels();
55 
56  std::unique_ptr<std::vector<sim::MCHitCollection>> hits_v(
57  new std::vector<sim::MCHitCollection>());
58  std::unique_ptr<std::vector<sim::MCWireCollection>> wires_v(
59  new std::vector<sim::MCWireCollection>());
60 
61  hits_v->reserve(nch);
62  wires_v->reserve(nch);
63  for (size_t ch = 0; ch < nch; ++ch) {
64 
65  hits_v->push_back(sim::MCHitCollection(ch));
66  wires_v->push_back(sim::MCWireCollection(ch));
67  }
68 
70  e.getByLabel(fLArG4ModuleName, simchArray);
71  if (!simchArray.isValid())
72  throw cet::exception(__PRETTY_FUNCTION__)
73  << "Did not find sim::SimChannel with a label: " << fLArG4ModuleName.c_str() << std::endl;
74 
75  // Loop over SimChannel
76  for (size_t simch_index = 0; simch_index < simchArray->size(); ++simch_index) {
77 
78  const art::Ptr<sim::SimChannel> simch_ptr(simchArray, simch_index);
79 
80  size_t ch = simch_ptr->Channel();
81 
82  if (ch >= hits_v->size())
83 
84  throw cet::exception(__PRETTY_FUNCTION__)
85  << "Channel number " << ch << " exceeds total # of channels: " << nch << std::endl;
86 
87  auto& mchits = hits_v->at(ch);
88  auto& mcwires = wires_v->at(ch);
89 
90  std::map<sim::MCEnDep, sim::MCWire> edep_wire_map;
91  auto tdc_ide_map = simch_ptr->TDCIDEMap();
92 
94  // Loop over stored data for a particular sim::SimChannel //
96 
97  //
98  // Read & Convert
99  //
100  if (fVerbose) std::cout << std::endl << "Processing Ch: " << ch << std::endl;
101 
102  for (auto const& tdc_ide_pair : tdc_ide_map) {
103 
104  auto const& tdc = tdc_ide_pair.first;
105  auto const& ide_v = tdc_ide_pair.second;
106 
107  for (auto const& ide : ide_v) {
109  edep.SetVertex(ide.x, ide.y, ide.z);
110  edep.SetEnergy(ide.energy);
111  edep.SetTrackId(ide.trackID);
112 
113  sim::MCWire wire;
114  wire.SetStartTDC(tdc);
115 
116  auto edep_iter = edep_wire_map.insert(std::make_pair(edep, wire));
117 
118  auto last_tdc =
119  (edep_iter).first->second.StartTDC() + (edep_iter).first->second.size() - 1;
120 
121  if (fVerbose) {
122 
123  if (edep_iter.second) std::cout << std::endl;
124 
125  std::cout << " Track: " << ide.trackID << " Vtx: " << ide.x << " " << ide.y << " "
126  << ide.z << " " << ide.energy << " ... @ TDC = " << tdc << " ... "
127  << ide.numElectrons << std::endl;
128  }
129 
130  if (!(edep_iter.second)) {
131  if (last_tdc + 1 != tdc) { continue; }
132  }
133  (edep_iter).first->second.push_back(ide.numElectrons);
134  } // Done looping over IDEs in a particular TDC
135  } // Done looping over TDC
136 
137  //
138  // Store
139  //
140  for (auto const& edep_wire_pair : edep_wire_map) {
141 
142  auto const& edep = edep_wire_pair.first;
143  auto const& wire = edep_wire_pair.second;
144 
145  // Create MCHit
147  float vtx[3] = {float(edep.Vertex()[0]), float(edep.Vertex()[1]), float(edep.Vertex()[2])};
148  hit.SetParticleInfo(vtx, edep.Energy(), edep.TrackId());
149 
150  double max_time = 0;
151  double qsum = 0;
152  double qmax = 0;
153  for (size_t wire_index = 0; wire_index < wire.size(); ++wire_index) {
154 
155  auto q = wire.at(wire_index);
156 
157  qsum += q;
158 
159  if (q > qmax) {
160  qmax = q;
161  max_time = wire.StartTDC() + wire_index;
162  }
163  }
164  hit.SetCharge(qsum, qmax);
165  hit.SetTime(max_time, 0);
166 
167  mchits.push_back(hit);
168 
169  if (!fMakeMCWire) continue;
170 
171  mcwires.push_back(wire);
172  } // End looping over all MCEnDep-MCWire pairs
173  } // End looping over all SimChannels
174 
175  simchArray.removeProduct();
176 
177  std::sort((*hits_v).begin(), (*hits_v).end());
178  e.put(std::move(hits_v));
179 
180  if (fMakeMCWire) {
181  std::sort((*wires_v).begin(), (*wires_v).end());
182  e.put(std::move(wires_v));
183  }
184  }
185 }
186 
void SetCharge(float qsum, float amp)
Setter function for charge/amplitude.
Definition: MCHit.h:47
void SetEnergy(float e)
Definition: MCDataHolder.h:36
void SetParticleInfo(const float vtx[], const float energy, const int trackId)
Setter function for partile info.
Definition: MCHit.h:61
void SetTime(const float peak, const float width)
Setter function for time.
Definition: MCHit.h:54
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
void SetStartTDC(const unsigned int start)
Setter function for time.
Definition: MCWire.h:34
void produce(art::Event &e) override
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
void SetTrackId(unsigned int id)
Definition: MCDataHolder.h:38
bool isValid() const noexcept
Definition: Handle.h:203
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
void SetVertex(float x, float y, float z)
Definition: MCDataHolder.h:29
std::string fLArG4ModuleName
bool removeProduct()
Definition: Handle.h:263
Double_t edep
Definition: macro.C:13
Detector simulation of raw signals on wires.
raw::ChannelID_t Channel() const
Returns the readout channel this object describes.
Definition: SimChannel.h:323
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
object containing MC truth information necessary for making RawDigits and doing back tracking ...
TDCIDEs_t const & TDCIDEMap() const
Returns all the deposited energy information as stored.
Definition: SimChannel.h:319
Float_t e
Definition: plot.C:35
MCHitFinder(fhicl::ParameterSet const &p)
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33