All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HitFinderAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // HitFinderAna class
4 //
5 // echurch@fnal.gov
6 //
7 // This algorithm is designed to analyze hits on wires after deconvolution
8 ////////////////////////////////////////////////////////////////////////
9 // Framework includes
10 #include "art/Framework/Core/ModuleMacros.h"
11 
12 // C++ includes
13 #include <utility>
14 #include <vector>
15 
16 // Framework includes
17 #include "art/Framework/Principal/Event.h"
18 #include "art/Framework/Principal/Handle.h"
19 #include "art/Framework/Services/Registry/ServiceHandle.h"
20 #include "art_root_io/TFileService.h"
21 #include "canvas/Persistency/Common/Ptr.h"
22 #include "fhiclcpp/ParameterSet.h"
23 #include "messagefacility/MessageLogger/MessageLogger.h"
24 
25 // LArSoft includes
26 #include "larcore/Geometry/Geometry.h"
27 #include "lardata/DetectorInfoServices/DetectorClocksService.h"
28 #include "lardataobj/RecoBase/Hit.h"
29 #include "larsim/MCCheater/BackTrackerService.h"
30 #include "larsim/MCCheater/ParticleInventoryService.h"
31 #include "nug4/ParticleNavigation/ParticleList.h"
32 #include "nusimdata/SimulationBase/MCParticle.h"
33 
34 #include "TTree.h"
35 
36 #include "art/Framework/Core/EDAnalyzer.h"
37 
38 #include <string>
39 
40 namespace geo {
41  class Geometry;
42 }
43 
44 ///Detector simulation of raw signals on wires
45 namespace hit {
46 
47  /// Base class for creation of raw signals on wires.
48  class HitFinderAna : public art::EDAnalyzer {
49  public:
50  explicit HitFinderAna(fhicl::ParameterSet const& pset);
51 
52  private:
53  /// read/write access to event
54  void analyze(const art::Event& evt);
55  void beginJob();
56 
57  std::string fFFTHitFinderModuleLabel;
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  //-------------------------------------------------
106  void
107  HitFinderAna::beginJob()
108  {
109  // get access to the TFile service
110  art::ServiceHandle<art::TFileService const> tfs;
111  fNp0 = 9000;
112  fNp1 = 9000;
113  fNp2 = 9000;
114 
115  fHTree = tfs->make<TTree>("HTree", "HTree");
116  fTimep0 = new Float_t[fNp0];
117  fTimep1 = new Float_t[fNp1];
118  fTimep2 = new Float_t[fNp2];
119  fWirep0 = new Int_t[fNp0];
120  fWirep1 = new Int_t[fNp1];
121  fWirep2 = new Int_t[fNp2];
122  fChgp0 = new Float_t[fNp0];
123  fChgp1 = new Float_t[fNp1];
124  fChgp2 = new Float_t[fNp2];
125  fXYZp0 = new Float_t[fNp0 * 3];
126  fXYZp1 = new Float_t[fNp1 * 3];
127  fXYZp2 = new Float_t[fNp2 * 3];
128 
129  fMCPdg0 = new Int_t[fNp0];
130  fMCPdg1 = new Int_t[fNp1];
131  fMCPdg2 = new Int_t[fNp2];
132  fMCTId0 = new Int_t[fNp0];
133  fMCTId1 = new Int_t[fNp1];
134  fMCTId2 = new Int_t[fNp2];
135  fMCE0 = new Float_t[fNp0];
136  fMCE1 = new Float_t[fNp1];
137  fMCE2 = new Float_t[fNp2];
138 
139  fHTree->Branch("HEvt", &fEvt, "HEvt/I");
140  fHTree->Branch("HRun", &fRun, "HRun/I");
141  fHTree->Branch("HNp0", &fNp0, "HNp0/I");
142  fHTree->Branch("HNp1", &fNp1, "HNp1/I");
143  fHTree->Branch("HNp2", &fNp2, "HNp2/I");
144  fHTree->Branch("HN3p0", &fN3p0, "HN3p0/I");
145  fHTree->Branch("HN3p1", &fN3p1, "HN3p1/I");
146  fHTree->Branch("HN3p2", &fN3p2, "HN3p2/I");
147  fHTree->Branch("Htp0", fTimep0, "Htp0[HNp0]/F");
148  fHTree->Branch("Htp1", fTimep1, "Htp1[HNp1]/F");
149  fHTree->Branch("Htp2", fTimep2, "Htp2[HNp2]/F");
150  fHTree->Branch("Hwp0", fWirep0, "Hwp0[HNp0]/I");
151  fHTree->Branch("Hwp1", fWirep1, "Hwp1[HNp1]/I");
152  fHTree->Branch("Hwp2", fWirep2, "Hwp2[HNp2]/I");
153  fHTree->Branch("Hchgp0", fChgp0, "Hchgp0[HNp0]/F");
154  fHTree->Branch("Hchgp1", fChgp1, "Hchgp1[HNp1]/F");
155  fHTree->Branch("Hchgp2", fChgp2, "Hchgp2[HNp2]/F");
156  fHTree->Branch("HMCXYZp0", fXYZp0, "HMCXYZp0[HN3p0]/F");
157  fHTree->Branch("HMCXYZp1", fXYZp1, "HMCXYZp1[HN3p1]/F");
158  fHTree->Branch("HMCXYZp2", fXYZp2, "HMCXYZp2[HN3p2]/F");
159  fHTree->Branch("HMCPdgp0", fMCPdg0, "HMCPdgp0[HNp0]/I");
160  fHTree->Branch("HMCPdgp1", fMCPdg1, "HMCPdgp1[HNp1]/I");
161  fHTree->Branch("HMCPdgp2", fMCPdg2, "HMCPdgp2[HNp2]/I");
162  fHTree->Branch("HMCTIdp0", fMCTId0, "HMCTIdp0[HNp0]/I");
163  fHTree->Branch("HMCTIdp1", fMCTId1, "HMCTIdp1[HNp1]/I");
164  fHTree->Branch("HMCTIdp2", fMCTId2, "HMCTIdp2[HNp2]/I");
165  fHTree->Branch("HMCEp0", fMCE0, "HMCEp0[HNp0]/F");
166  fHTree->Branch("HMCEp1", fMCE1, "HMCEp1[HNp1]/F");
167  fHTree->Branch("HMCEp2", fMCE2, "HMCEp2[HNp2]/F");
168 
169  return;
170  }
171 
172  //-------------------------------------------------
173  void
174  HitFinderAna::analyze(const art::Event& evt)
175  {
176 
177  if (evt.isRealData()) {
178  throw cet::exception("HitFinderAna: ") << "Not for use on Data yet...\n";
179  }
180 
181  art::Handle<std::vector<recob::Hit>> hitHandle;
182  evt.getByLabel(fFFTHitFinderModuleLabel, hitHandle);
183 
184  art::ServiceHandle<cheat::BackTrackerService const> bt_serv;
185  art::ServiceHandle<cheat::ParticleInventoryService const> pi_serv;
186  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
187 
188  sim::ParticleList const& _particleList = pi_serv->ParticleList();
189 
190  MF_LOG_VERBATIM("HitFinderAna") << _particleList;
191 
192  art::ServiceHandle<geo::Geometry const> geom;
193 
194  std::map<geo::PlaneID, std::vector<art::Ptr<recob::Hit>>> planeIDToHits;
195  for (size_t h = 0; h < hitHandle->size(); ++h)
196  planeIDToHits[hitHandle->at(h).WireID().planeID()].push_back(
197  art::Ptr<recob::Hit>(hitHandle, h));
198 
199  for (auto mapitr : planeIDToHits) {
200  fNp0 = 0;
201  fN3p0 = 0;
202  fNp1 = 0;
203  fN3p1 = 0;
204  fNp2 = 0;
205  fN3p2 = 0;
206 
207  geo::PlaneID pid = mapitr.first;
208  auto itr = mapitr.second.begin();
209  while (itr != mapitr.second.end()) {
210 
211  fRun = evt.run();
212  fEvt = evt.id().event();
213 
214  std::vector<sim::TrackIDE> trackides = bt_serv->HitToTrackIDEs(clockData, *itr);
215  std::vector<sim::TrackIDE>::iterator idesitr = trackides.begin();
216  std::vector<double> xyz = bt_serv->HitToXYZ(clockData, *itr);
217 
218  if (pid.Plane == 0 && fNp0 < 9000) {
219  fTimep0[fNp0] = (*itr)->PeakTime();
220  fWirep0[fNp0] = (*itr)->WireID().Wire;
221  fChgp0[fNp0] = (*itr)->Integral();
222 
223  for (unsigned int kk = 0; kk < 3; ++kk) {
224  fXYZp0[fNp0 * 3 + kk] = xyz[kk];
225  }
226 
227  while (idesitr != trackides.end()) {
228  fMCTId0[fNp0] = (*idesitr).trackID;
229  if (_particleList.find((*idesitr).trackID) != _particleList.end()) {
230  const simb::MCParticle* particle = _particleList.at((*idesitr).trackID);
231  fMCPdg0[fNp0] = particle->PdgCode();
232  fMCE0[fNp0] = particle->E();
233  }
234  idesitr++;
235  }
236 
237  ++fNp0;
238  }
239 
240  else if (pid.Plane == 1 && fNp1 < 9000) {
241  fTimep1[fNp1] = (*itr)->PeakTime();
242  fWirep1[fNp1] = (*itr)->WireID().Wire;
243  fChgp1[fNp1] = (*itr)->Integral();
244 
245  for (unsigned int kk = 0; kk < 3; ++kk) {
246  fXYZp1[fNp1 * 3 + kk] = xyz[kk];
247  }
248 
249  while (idesitr != trackides.end()) {
250  fMCTId1[fNp1] = (*idesitr).trackID;
251  if (_particleList.find((*idesitr).trackID) != _particleList.end()) {
252  const simb::MCParticle* particle = _particleList.at((*idesitr).trackID);
253  fMCPdg1[fNp1] = particle->PdgCode();
254  fMCE1[fNp1] = particle->E();
255  }
256  idesitr++;
257  }
258  ++fNp1;
259  }
260 
261  else if (pid.Plane == 2 && fNp2 < 9000) {
262  fTimep2[fNp2] = (*itr)->PeakTime();
263  fWirep2[fNp2] = (*itr)->WireID().Wire;
264  fChgp2[fNp2] = (*itr)->Integral();
265 
266  for (unsigned int kk = 0; kk < 3; ++kk) {
267  fXYZp2[fNp2 * 3 + kk] = xyz[kk];
268  }
269 
270  while (idesitr != trackides.end()) {
271  fMCTId2[fNp2] = (*idesitr).trackID;
272  if (_particleList.find((*idesitr).trackID) != _particleList.end()) {
273  const simb::MCParticle* particle = _particleList.at((*idesitr).trackID);
274  fMCPdg2[fNp2] = particle->PdgCode();
275  fMCE2[fNp2] = particle->E();
276  }
277  idesitr++;
278  }
279  ++fNp2;
280  }
281 
282  fN3p0 = 3 * fNp0;
283  fN3p1 = 3 * fNp1;
284  fN3p2 = 3 * fNp2;
285 
286  fHTree->Fill();
287  itr++;
288  } // loop on Hits
289  } // loop on map
290 
291  return;
292  } // end analyze method
293 
294 } // end namespace
295 
296 DEFINE_ART_MODULE(hit::HitFinderAna)
Hit.h
Declaration of signal hit object.
geo::PlaneID
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
hit::HitFinderAna::fFFTHitFinderModuleLabel
std::string fFFTHitFinderModuleLabel
Definition: HitFinderAna_module.cc:57
hit
process_name hit
Definition: cheaterreco.fcl:51
hit::HitFinderAna::fLArG4ModuleLabel
std::string fLArG4ModuleLabel
Definition: HitFinderAna_module.cc:58
hit::HitFinderAna::HitFinderAna
HitFinderAna(fhicl::ParameterSet const &pset)
Definition: HitFinderAna_module.cc:99
hit::HitFinderAna
Base class for creation of raw signals on wires.
Definition: HitFinderAna_module.cc:48
h
while getopts h
Definition: getStatistics.sh:13
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
hit::HitFinderAna::analyze
void analyze(const art::Event &evt)
read/write access to event
Definition: HitFinderAna_module.cc:174
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
Geometry.h
art framework interface to geometry description
Generated by   doxygen 1.8.5