LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
HitFinderAna_module.cc
Go to the documentation of this file.
1 //
3 // HitFinderAna class
4 //
5 // echurch@fnal.gov
6 //
7 // This algorithm is designed to analyze hits on wires after deconvolution
9 // Framework includes
11 
12 // C++ includes
13 #include <utility>
14 #include <vector>
15 
16 // Framework includes
20 #include "art_root_io/TFileService.h"
22 #include "fhiclcpp/ParameterSet.h"
24 
25 // LArSoft includes
33 
34 #include "TTree.h"
35 
37 
38 #include <string>
39 
40 namespace geo {
41  class Geometry;
42 }
43 
45 namespace hit {
46 
48  class HitFinderAna : public art::EDAnalyzer {
49  public:
50  explicit HitFinderAna(fhicl::ParameterSet const& pset);
51 
52  private:
54  void analyze(const art::Event& evt);
55  void beginJob();
56 
58  std::string fLArG4ModuleLabel;
59 
60  TTree* fHTree;
61  Int_t fRun;
62  Int_t fEvt;
63  Int_t fNp0;
64  Int_t fNp1;
65  Int_t fNp2;
66  Int_t fN3p0;
67  Int_t fN3p1;
68  Int_t fN3p2;
69  Float_t* fTimep0;
70  Float_t* fTimep1;
71  Float_t* fTimep2;
72  Int_t* fWirep0;
73  Int_t* fWirep1;
74  Int_t* fWirep2;
75  Float_t* fChgp0;
76  Float_t* fChgp1;
77  Float_t* fChgp2;
78  Float_t* fXYZp0;
79  Float_t* fXYZp1;
80  Float_t* fXYZp2;
81 
82  Int_t* fMCPdg0;
83  Int_t* fMCTId0;
84  Float_t* fMCE0;
85  Int_t* fMCPdg1;
86  Int_t* fMCTId1;
87  Float_t* fMCE1;
88  Int_t* fMCPdg2;
89  Int_t* fMCTId2;
90  Float_t* fMCE2;
91 
92  }; // class HitFinderAna
93 
94 }
95 
96 namespace hit {
97 
98  //-------------------------------------------------
99  HitFinderAna::HitFinderAna(fhicl::ParameterSet const& pset) : EDAnalyzer(pset)
100  {
101  fFFTHitFinderModuleLabel = pset.get<std::string>("HitsModuleLabel");
102  fLArG4ModuleLabel = pset.get<std::string>("LArGeantModuleLabel");
103  }
104 
105  //-------------------------------------------------
107  {
108  // get access to the TFile service
110  fNp0 = 9000;
111  fNp1 = 9000;
112  fNp2 = 9000;
113 
114  fHTree = tfs->make<TTree>("HTree", "HTree");
115  fTimep0 = new Float_t[fNp0];
116  fTimep1 = new Float_t[fNp1];
117  fTimep2 = new Float_t[fNp2];
118  fWirep0 = new Int_t[fNp0];
119  fWirep1 = new Int_t[fNp1];
120  fWirep2 = new Int_t[fNp2];
121  fChgp0 = new Float_t[fNp0];
122  fChgp1 = new Float_t[fNp1];
123  fChgp2 = new Float_t[fNp2];
124  fXYZp0 = new Float_t[fNp0 * 3];
125  fXYZp1 = new Float_t[fNp1 * 3];
126  fXYZp2 = new Float_t[fNp2 * 3];
127 
128  fMCPdg0 = new Int_t[fNp0];
129  fMCPdg1 = new Int_t[fNp1];
130  fMCPdg2 = new Int_t[fNp2];
131  fMCTId0 = new Int_t[fNp0];
132  fMCTId1 = new Int_t[fNp1];
133  fMCTId2 = new Int_t[fNp2];
134  fMCE0 = new Float_t[fNp0];
135  fMCE1 = new Float_t[fNp1];
136  fMCE2 = new Float_t[fNp2];
137 
138  fHTree->Branch("HEvt", &fEvt, "HEvt/I");
139  fHTree->Branch("HRun", &fRun, "HRun/I");
140  fHTree->Branch("HNp0", &fNp0, "HNp0/I");
141  fHTree->Branch("HNp1", &fNp1, "HNp1/I");
142  fHTree->Branch("HNp2", &fNp2, "HNp2/I");
143  fHTree->Branch("HN3p0", &fN3p0, "HN3p0/I");
144  fHTree->Branch("HN3p1", &fN3p1, "HN3p1/I");
145  fHTree->Branch("HN3p2", &fN3p2, "HN3p2/I");
146  fHTree->Branch("Htp0", fTimep0, "Htp0[HNp0]/F");
147  fHTree->Branch("Htp1", fTimep1, "Htp1[HNp1]/F");
148  fHTree->Branch("Htp2", fTimep2, "Htp2[HNp2]/F");
149  fHTree->Branch("Hwp0", fWirep0, "Hwp0[HNp0]/I");
150  fHTree->Branch("Hwp1", fWirep1, "Hwp1[HNp1]/I");
151  fHTree->Branch("Hwp2", fWirep2, "Hwp2[HNp2]/I");
152  fHTree->Branch("Hchgp0", fChgp0, "Hchgp0[HNp0]/F");
153  fHTree->Branch("Hchgp1", fChgp1, "Hchgp1[HNp1]/F");
154  fHTree->Branch("Hchgp2", fChgp2, "Hchgp2[HNp2]/F");
155  fHTree->Branch("HMCXYZp0", fXYZp0, "HMCXYZp0[HN3p0]/F");
156  fHTree->Branch("HMCXYZp1", fXYZp1, "HMCXYZp1[HN3p1]/F");
157  fHTree->Branch("HMCXYZp2", fXYZp2, "HMCXYZp2[HN3p2]/F");
158  fHTree->Branch("HMCPdgp0", fMCPdg0, "HMCPdgp0[HNp0]/I");
159  fHTree->Branch("HMCPdgp1", fMCPdg1, "HMCPdgp1[HNp1]/I");
160  fHTree->Branch("HMCPdgp2", fMCPdg2, "HMCPdgp2[HNp2]/I");
161  fHTree->Branch("HMCTIdp0", fMCTId0, "HMCTIdp0[HNp0]/I");
162  fHTree->Branch("HMCTIdp1", fMCTId1, "HMCTIdp1[HNp1]/I");
163  fHTree->Branch("HMCTIdp2", fMCTId2, "HMCTIdp2[HNp2]/I");
164  fHTree->Branch("HMCEp0", fMCE0, "HMCEp0[HNp0]/F");
165  fHTree->Branch("HMCEp1", fMCE1, "HMCEp1[HNp1]/F");
166  fHTree->Branch("HMCEp2", fMCE2, "HMCEp2[HNp2]/F");
167 
168  return;
169  }
170 
171  //-------------------------------------------------
173  {
174 
175  if (evt.isRealData()) {
176  throw cet::exception("HitFinderAna: ") << "Not for use on Data yet...\n";
177  }
178 
180  evt.getByLabel(fFFTHitFinderModuleLabel, hitHandle);
181 
184  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
185 
186  sim::ParticleList const& _particleList = pi_serv->ParticleList();
187 
188  MF_LOG_VERBATIM("HitFinderAna") << _particleList;
189 
191 
192  std::map<geo::PlaneID, std::vector<art::Ptr<recob::Hit>>> planeIDToHits;
193  for (size_t h = 0; h < hitHandle->size(); ++h)
194  planeIDToHits[hitHandle->at(h).WireID().planeID()].push_back(
195  art::Ptr<recob::Hit>(hitHandle, h));
196 
197  for (auto mapitr : planeIDToHits) {
198  fNp0 = 0;
199  fN3p0 = 0;
200  fNp1 = 0;
201  fN3p1 = 0;
202  fNp2 = 0;
203  fN3p2 = 0;
204 
205  geo::PlaneID pid = mapitr.first;
206  auto itr = mapitr.second.begin();
207  while (itr != mapitr.second.end()) {
208 
209  fRun = evt.run();
210  fEvt = evt.id().event();
211 
212  std::vector<sim::TrackIDE> trackides = bt_serv->HitToTrackIDEs(clockData, *itr);
213  std::vector<sim::TrackIDE>::iterator idesitr = trackides.begin();
214  std::vector<double> xyz = bt_serv->HitToXYZ(clockData, *itr);
215 
216  if (pid.Plane == 0 && fNp0 < 9000) {
217  fTimep0[fNp0] = (*itr)->PeakTime();
218  fWirep0[fNp0] = (*itr)->WireID().Wire;
219  fChgp0[fNp0] = (*itr)->Integral();
220 
221  for (unsigned int kk = 0; kk < 3; ++kk) {
222  fXYZp0[fNp0 * 3 + kk] = xyz[kk];
223  }
224 
225  while (idesitr != trackides.end()) {
226  fMCTId0[fNp0] = (*idesitr).trackID;
227  if (_particleList.find((*idesitr).trackID) != _particleList.end()) {
228  const simb::MCParticle* particle = _particleList.at((*idesitr).trackID);
229  fMCPdg0[fNp0] = particle->PdgCode();
230  fMCE0[fNp0] = particle->E();
231  }
232  idesitr++;
233  }
234 
235  ++fNp0;
236  }
237 
238  else if (pid.Plane == 1 && fNp1 < 9000) {
239  fTimep1[fNp1] = (*itr)->PeakTime();
240  fWirep1[fNp1] = (*itr)->WireID().Wire;
241  fChgp1[fNp1] = (*itr)->Integral();
242 
243  for (unsigned int kk = 0; kk < 3; ++kk) {
244  fXYZp1[fNp1 * 3 + kk] = xyz[kk];
245  }
246 
247  while (idesitr != trackides.end()) {
248  fMCTId1[fNp1] = (*idesitr).trackID;
249  if (_particleList.find((*idesitr).trackID) != _particleList.end()) {
250  const simb::MCParticle* particle = _particleList.at((*idesitr).trackID);
251  fMCPdg1[fNp1] = particle->PdgCode();
252  fMCE1[fNp1] = particle->E();
253  }
254  idesitr++;
255  }
256  ++fNp1;
257  }
258 
259  else if (pid.Plane == 2 && fNp2 < 9000) {
260  fTimep2[fNp2] = (*itr)->PeakTime();
261  fWirep2[fNp2] = (*itr)->WireID().Wire;
262  fChgp2[fNp2] = (*itr)->Integral();
263 
264  for (unsigned int kk = 0; kk < 3; ++kk) {
265  fXYZp2[fNp2 * 3 + kk] = xyz[kk];
266  }
267 
268  while (idesitr != trackides.end()) {
269  fMCTId2[fNp2] = (*idesitr).trackID;
270  if (_particleList.find((*idesitr).trackID) != _particleList.end()) {
271  const simb::MCParticle* particle = _particleList.at((*idesitr).trackID);
272  fMCPdg2[fNp2] = particle->PdgCode();
273  fMCE2[fNp2] = particle->E();
274  }
275  idesitr++;
276  }
277  ++fNp2;
278  }
279 
280  fN3p0 = 3 * fNp0;
281  fN3p1 = 3 * fNp1;
282  fN3p2 = 3 * fNp2;
283 
284  fHTree->Fill();
285  itr++;
286  } // loop on Hits
287  } // loop on map
288 
289  return;
290  } // end analyze method
291 
292 } // end namespace
293 
double E(const int i=0) const
Definition: MCParticle.h:234
intermediate_table::iterator iterator
int PdgCode() const
Definition: MCParticle.h:213
std::vector< sim::TrackIDE > HitToTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
Declaration of signal hit object.
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
mapped_type at(const key_type &key)
Definition: ParticleList.h:330
std::string fFFTHitFinderModuleLabel
bool isRealData() const
Definition: Event.cc:53
std::string fLArG4ModuleLabel
std::vector< double > HitToXYZ(detinfo::DetectorClocksData const &clockData, const recob::Hit &hit) const
Particle class.
Base class for creation of raw signals on wires.
iterator find(const key_type &key)
Definition: ParticleList.h:318
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
void beginJob()
Definition: Breakpoints.cc:14
T get(std::string const &key) const
Definition: ParameterSet.h:314
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
Detector simulation of raw signals on wires.
const sim::ParticleList & ParticleList() const
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
#define MF_LOG_VERBATIM(category)
void analyze(const art::Event &evt)
read/write access to event
EventNumber_t event() const
Definition: EventID.h:116
TCEvent evt
Definition: DataStructs.cxx:8
Particle list in DetSim contains Monte Carlo particle information.
RunNumber_t run() const
Definition: Event.cc:29
Namespace collecting geometry-related classes utilities.
EventID id() const
Definition: Event.cc:23
art framework interface to geometry description
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33