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