LArSoft  v08_35_01
Liquid Argon Software toolkit - http://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 
14 #include "fhiclcpp/ParameterSet.h"
15 
16 #include <memory>
17 
19 
24 
25 namespace hit {
26 
27  class MCHitFinder;
28 
29  class MCHitFinder : public art::EDProducer {
30  public:
31  explicit MCHitFinder(fhicl::ParameterSet const & p);
32 
33  private:
34  void produce(art::Event & e) override;
35 
36  std::string fLArG4ModuleName;
37  bool fVerbose;
39 
40  };
41 
42 
44  : EDProducer{p}
45  {
46  fLArG4ModuleName = p.get<std::string>("LArG4ModuleName");
47  fMakeMCWire = p.get<bool>("MakeMCWire");
48  fVerbose = p.get<bool>("Verbose");
49  produces< std::vector< sim::MCHitCollection> >();
50  if(fMakeMCWire)
51  produces< std::vector< sim::MCWireCollection> >();
52  }
53 
55  {
56 
58 
59  const unsigned int nch = geo->Nchannels();
60 
61  std::unique_ptr< std::vector<sim::MCHitCollection> > hits_v ( new std::vector<sim::MCHitCollection>() );
62  std::unique_ptr< std::vector<sim::MCWireCollection> > wires_v ( new std::vector<sim::MCWireCollection>() );
63 
64  hits_v->reserve(nch);
65  wires_v->reserve(nch);
66  for(size_t ch=0; ch<nch; ++ch) {
67 
68  hits_v->push_back(sim::MCHitCollection(ch));
69  wires_v->push_back(sim::MCWireCollection(ch));
70 
71  }
72 
74  e.getByLabel(fLArG4ModuleName,simchArray);
75  if(!simchArray.isValid())
76  throw cet::exception(__PRETTY_FUNCTION__)
77  << "Did not find sim::SimChannel with a label: " << fLArG4ModuleName.c_str() << std::endl;
78 
79  // Loop over SimChannel
80  for(size_t simch_index=0; simch_index<simchArray->size(); ++simch_index) {
81 
82  const art::Ptr<sim::SimChannel> simch_ptr(simchArray,simch_index);
83 
84  size_t ch = simch_ptr->Channel();
85 
86  if(ch >= hits_v->size())
87 
88  throw cet::exception(__PRETTY_FUNCTION__)
89  << "Channel number " << ch << " exceeds total # of channels: " << nch << std::endl;
90 
91  auto &mchits = hits_v->at(ch);
92  auto &mcwires = wires_v->at(ch);
93 
94  std::map<sim::MCEnDep,sim::MCWire> edep_wire_map;
95  auto tdc_ide_map = simch_ptr->TDCIDEMap();
96 
98  // Loop over stored data for a particular sim::SimChannel //
100 
101  //
102  // Read & Convert
103  //
104  if(fVerbose)
105  std::cout<<std::endl<<"Processing Ch: "<<ch<<std::endl;
106 
107  for(auto const& tdc_ide_pair : tdc_ide_map) {
108 
109  auto const& tdc = tdc_ide_pair.first;
110  auto const& ide_v = tdc_ide_pair.second;
111 
112  for(auto const& ide : ide_v) {
114  edep.SetVertex(ide.x,ide.y,ide.z);
115  edep.SetEnergy(ide.energy);
116  edep.SetTrackId(ide.trackID);
117 
118  sim::MCWire wire;
119  wire.SetStartTDC(tdc);
120 
121  auto edep_iter = edep_wire_map.insert(std::make_pair(edep,wire));
122 
123  //auto edep_iter = edep_wire_map.find(edep);
124 
125  auto last_tdc = (edep_iter).first->second.StartTDC() + (edep_iter).first->second.size() - 1;
126 
127  if(fVerbose) {
128 
129  if( edep_iter.second ) std::cout<<std::endl;
130 
131  std::cout<<" Track: "<<ide.trackID
132  <<" Vtx: " <<ide.x << " " << ide.y << " " <<ide.z << " " <<ide.energy
133  << " ... @ TDC = "<<tdc<<" ... "<<ide.numElectrons << std::endl;
134  }
135 
136  if( !(edep_iter.second) ) {
137 
138  if( last_tdc+1 != tdc ) {
139  /*
140  std::cerr
141  << "Found discontinuous TDC! "
142  << " Last (ADC @ TDC): " << (*(edep_iter.first->second.rbegin())) << " @ " << last_tdc
143  << " while current (ADC @ TDC): " << (ide.numElectrons * detp->ElectronsToADC()) << " @ " << tdc
144  << " ... skipping to store!"
145  << std::endl;
146  */
147  continue;
148  }
149  }
150  (edep_iter).first->second.push_back(ide.numElectrons);
151  }// Done looping over IDEs in a particular TDC
152  }// Done looping over TDC
153 
154  //
155  // Store
156  //
157  for(auto const& edep_wire_pair : edep_wire_map) {
158 
159  auto const& edep = edep_wire_pair.first;
160  auto const& wire = edep_wire_pair.second;
161 
162  // Create MCHit
164  float vtx[3] = {float(edep.Vertex()[0]),
165  float(edep.Vertex()[1]),
166  float(edep.Vertex()[2])};
167  hit.SetParticleInfo(vtx, edep.Energy(), edep.TrackId());
168 
169  double max_time = 0;
170  double qsum = 0;
171  double qmax = 0;
172  for(size_t wire_index=0; wire_index < wire.size(); ++wire_index) {
173 
174  auto q = wire.at(wire_index);
175 
176  qsum += q;
177 
178  if( q > qmax) { qmax = q; max_time = wire.StartTDC() + wire_index; }
179 
180  }
181  hit.SetCharge(qsum,qmax);
182  hit.SetTime(max_time,0);
183 
184  mchits.push_back(hit);
185 
186  if(!fMakeMCWire) continue;
187 
188  mcwires.push_back(wire);
189  } // End looping over all MCEnDep-MCWire pairs
190  } // End looping over all SimChannels
191 
192  e.removeCachedProduct(simchArray);
193 
194  std::sort((*hits_v).begin(),(*hits_v).end());
195  e.put(std::move(hits_v));
196 
197  if(fMakeMCWire) {
198  std::sort((*wires_v).begin(),(*wires_v).end());
199  e.put(std::move(wires_v));
200  }
201 
202  }
203 }
void SetCharge(float qsum, float amp)
Setter function for charge/amplitude.
Definition: MCHit.h:54
void SetEnergy(float e)
Definition: MCDataHolder.h:39
void SetParticleInfo(const float vtx[], const float energy, const int trackId)
Setter function for partile info.
Definition: MCHit.h:64
void SetTime(const float peak, const float width)
Setter function for time.
Definition: MCHit.h:57
EDProducer()=default
void SetStartTDC(const unsigned int start)
Setter function for time.
Definition: MCWire.h:42
void produce(art::Event &e) override
void SetTrackId(unsigned int id)
Definition: MCDataHolder.h:41
unsigned int Nchannels() const
Returns the number of TPC readout channels in the detector.
bool isValid() const
Definition: Handle.h:183
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:437
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:56
void SetVertex(float x, float y, float z)
Definition: MCDataHolder.h:32
std::string fLArG4ModuleName
Double_t edep
Definition: macro.C:13
bool removeCachedProduct(Handle< PROD > &) const
Definition: DataViewImpl.h:933
Detector simulation of raw signals on wires.
raw::ChannelID_t Channel() const
Returns the readout channel this object describes.
Definition: SimChannel.h:329
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:328
Float_t e
Definition: plot.C:34
MCHitFinder(fhicl::ParameterSet const &p)
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
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33