CheckT0_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: CheckT0
3 // Plugin Type: analyzer (art v3_02_06)
4 // File: CheckT0_module.cc
5 //
6 // Generated at Tue Jul 30 21:43:10 2019 by Tingjun Yang using cetskelgen
7 // from cetlib version v3_07_02.
8 ////////////////////////////////////////////////////////////////////////
9 
17 #include "canvas/Persistency/Common/FindManyP.h"
19 #include "fhiclcpp/ParameterSet.h"
21 #include "art_root_io/TFileService.h"
22 
34 
35 #include "TTree.h"
36 
37 #include <vector>
38 
39 namespace pdsp {
40  class CheckT0;
41 }
42 
43 
45 public:
46  explicit CheckT0(fhicl::ParameterSet const& p);
47  // The compiler-generated destructor is fine for non-base
48  // classes without bare pointers or other resource use.
49 
50  // Plugins should not be copied or assigned.
51  CheckT0(CheckT0 const&) = delete;
52  CheckT0(CheckT0&&) = delete;
53  CheckT0& operator=(CheckT0 const&) = delete;
54  CheckT0& operator=(CheckT0&&) = delete;
55  void beginJob() override;
56 
57  // Required functions.
58  void analyze(art::Event const& e) override;
59 
60 private:
61 
62  // Declare member data here.
63  TTree *t0tree;
64  int run;
65  int event;
66  std::vector<int> trackid;
67  std::vector<double> t0crt2;
68  std::vector<double> t0crt1;
69  std::vector<double> t0pandora;
70  std::vector<double> t0anodep;
71  std::vector<double> t0truth;
72  std::vector<double> trackstartx;
73  std::vector<double> trackstarty;
74  std::vector<double> trackstartz;
75  std::vector<double> trackendx;
76  std::vector<double> trackendy;
77  std::vector<double> trackendz;
78  std::vector<double> trackstartx_sce;
79  std::vector<double> trackstarty_sce;
80  std::vector<double> trackstartz_sce;
81  std::vector<double> trackendx_sce;
82  std::vector<double> trackendy_sce;
83  std::vector<double> trackendz_sce;
84  std::vector<double> crt1x;
85  std::vector<double> crt1y;
86  std::vector<double> crt1z;
87  std::vector<double> crt2x0;
88  std::vector<double> crt2y0;
89  std::vector<double> crt2z0;
90  std::vector<double> crt2x1;
91  std::vector<double> crt2y1;
92  std::vector<double> crt2z1;
94  std::vector<double> rdtimestamp_digits;
95 
96 };
97 
98 
100  : EDAnalyzer{p} // ,
101  // More initializers here.
102 {
103 
104 }
105 
107 {
108  // Implementation of required member function here.
109  run = e.run();
110  event = e.id().event();
111  trackid.clear();
112  t0crt2.clear();
113  t0crt1.clear();
114  t0pandora.clear();
115  t0anodep.clear();
116  t0truth.clear();
117  trackstartx.clear();
118  trackstarty.clear();
119  trackstartz.clear();
120  trackendx.clear();
121  trackendy.clear();
122  trackendz.clear();
123  trackstartx_sce.clear();
124  trackstarty_sce.clear();
125  trackstartz_sce.clear();
126  trackendx_sce.clear();
127  trackendy_sce.clear();
128  trackendz_sce.clear();
129  crt1x.clear();
130  crt1y.clear();
131  crt1z.clear();
132  crt2x0.clear();
133  crt2y0.clear();
134  crt2z0.clear();
135  crt2x1.clear();
136  crt2y1.clear();
137  crt2z1.clear();
138  rdtimestamp_evt = -1;
139  rdtimestamp_digits.clear();
140 
141  //Services
144  //Space charge service
145  auto const* SCE = lar::providerFrom<spacecharge::SpaceChargeService>();
146 
147  // Reconstruciton information
148  std::vector < art::Ptr < recob::Track > > trackList;
149  auto trackListHandle = e.getHandle< std::vector < recob::Track > >("pandoraTrack");
150  if (trackListHandle) {
151  art::fill_ptr_vector(trackList, trackListHandle);
152  }
153  else return;
154 
155  //Get hits associated with track
156  art::FindManyP < recob::Hit > hitsFromTrack(trackListHandle, e, "pandoraTrack");
157 
158  //Get PFParticles
159  auto pfpListHandle = e.getHandle< std::vector<recob::PFParticle> >("pandora");
160 
161  //Get Track-PFParticle association
162  art::FindManyP<recob::PFParticle> fmpfp(trackListHandle, e, "pandoraTrack");
163 
164  //Get Pandora T0-PFParticle association
165  art::FindManyP<anab::T0> fmt0pandora(pfpListHandle, e, "pandora");
166 
167  //Get Anode Piercing T0
168  art::FindManyP<anab::T0> fmt0anodepiercer(pfpListHandle, e, "anodepiercerst0");
169 
170  //Get 2-CRT T0
171  art::FindManyP<anab::T0> fmt0crt2(trackListHandle, e, "crtreco");
172  art::FindManyP<anab::CosmicTag> fmctcrt2(trackListHandle, e, "crtreco");
173 
174  //Get 1-CRT T0
175  art::FindManyP<anab::T0> fmt0crt1(trackListHandle, e, "crttag");
176  art::FindManyP<anab::CosmicTag> fmctcrt1(trackListHandle, e, "crttag");
177 
178  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService>()->DataFor(e);
179  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService>()->DataFor(e, clockData);
180 
181  for (auto const& track : trackList){
182  int this_trackid = track.key();
183  double this_t0crt2 = -DBL_MAX;
184  double this_t0crt1 = -DBL_MAX;
185  double this_t0pandora = -DBL_MAX;
186  double this_t0anodep = -DBL_MAX;
187  double this_t0truth = -DBL_MAX;
188  double this_trackstartx = -DBL_MAX;
189  double this_trackstarty = -DBL_MAX;
190  double this_trackstartz = -DBL_MAX;
191  double this_trackendx = -DBL_MAX;
192  double this_trackendy = -DBL_MAX;
193  double this_trackendz = -DBL_MAX;
194  double this_trackstartx_sce = -DBL_MAX;
195  double this_trackstarty_sce = -DBL_MAX;
196  double this_trackstartz_sce = -DBL_MAX;
197  double this_trackendx_sce = -DBL_MAX;
198  double this_trackendy_sce = -DBL_MAX;
199  double this_trackendz_sce = -DBL_MAX;
200  double this_crt1x = -DBL_MAX;
201  double this_crt1y = -DBL_MAX;
202  double this_crt1z = -DBL_MAX;
203  double this_crt2x0 = -DBL_MAX;
204  double this_crt2y0 = -DBL_MAX;
205  double this_crt2z0 = -DBL_MAX;
206  double this_crt2x1 = -DBL_MAX;
207  double this_crt2y1 = -DBL_MAX;
208  double this_crt2z1 = -DBL_MAX;
209 
210  if (fmt0crt2.isValid()){
211  auto const& vt0crt2 = fmt0crt2.at(track.key());
212  if (!vt0crt2.empty()) this_t0crt2 = vt0crt2[0]->Time();
213  }
214 
215  if (fmt0crt1.isValid()){
216  auto const& vt0crt1 = fmt0crt1.at(track.key());
217  if (!vt0crt1.empty()) this_t0crt1 = vt0crt1[0]->Time();
218  }
219 
220  if (fmctcrt2.isValid()){
221  auto const& vctcrt2 = fmctcrt2.at(track.key());
222  if (!vctcrt2.empty()){
223  this_crt2x0 = vctcrt2[0]->EndPoint1()[0];
224  this_crt2y0 = vctcrt2[0]->EndPoint1()[1];
225  this_crt2z0 = vctcrt2[0]->EndPoint1()[2];
226  this_crt2x1 = vctcrt2[0]->EndPoint2()[0];
227  this_crt2y1 = vctcrt2[0]->EndPoint2()[1];
228  this_crt2z1 = vctcrt2[0]->EndPoint2()[2];
229  }
230  }
231 
232  if (fmctcrt1.isValid()){
233  auto const& vctcrt1 = fmctcrt1.at(track.key());
234  if (!vctcrt1.empty()){
235  this_crt1x = vctcrt1[0]->EndPoint1()[0];
236  this_crt1y = vctcrt1[0]->EndPoint1()[1];
237  this_crt1z = vctcrt1[0]->EndPoint1()[2];
238  }
239  }
240 
241  if (fmpfp.isValid()){
242  auto const &pfps = fmpfp.at(track.key());
243  if(!pfps.empty()){
244  //Find T0 for PFParticle
245  if (fmt0pandora.isValid()){
246  auto const &t0s = fmt0pandora.at(pfps[0].key());
247  if(!t0s.empty()){
248  this_t0pandora = t0s[0]->Time();
249  }
250  }
251  if (fmt0anodepiercer.isValid()){
252  auto const &t0aps = fmt0anodepiercer.at(pfps[0].key());
253  if (!t0aps.empty()){
254  this_t0anodep = t0aps[0]->Time();
255  }
256  }
257  }
258  }
259  if (this_t0crt2 > -DBL_MAX ||
260  this_t0crt1 > -DBL_MAX ||
261  this_t0pandora > -DBL_MAX ||
262  this_t0anodep > -DBL_MAX){
263 
264  auto const & allHits = hitsFromTrack.at(track.key());
265 
266  if (!e.isRealData()){
267  // Find true particle for reconstructed track
268  int TrackID = 0;
269  std::map<int,double> trkide;
270  for(size_t h = 0; h < allHits.size(); ++h){
271  art::Ptr<recob::Hit> hit = allHits[h];
272  std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
273  for(size_t e = 0; e < TrackIDs.size(); ++e){
274  trkide[TrackIDs[e].trackID] += TrackIDs[e].energy;
275  }
276  }
277  // Work out which IDE despoited the most charge in the hit if there was more than one.
278  double maxe = -1;
279  double tote = 0;
280  for (std::map<int,double>::iterator ii = trkide.begin(); ii!=trkide.end(); ++ii){
281  tote += ii->second;
282  if ((ii->second)>maxe){
283  maxe = ii->second;
284  TrackID = ii->first;
285  }
286  }
287  // Now have trackID, so get PdG code and T0 etc.
288  const simb::MCParticle *particle = pi_serv->TrackIdToParticle_P(TrackID);
289  if (particle){
290  this_t0truth = particle->T();
291  }
292  }
293 
294  this_trackstartx = track->Vertex().X();
295  this_trackstarty = track->Vertex().Y();
296  this_trackstartz = track->Vertex().Z();
297  this_trackendx = track->End().X();
298  this_trackendy = track->End().Y();
299  this_trackendz = track->End().Z();
300 
301  if (std::abs(this_t0pandora+DBL_MAX)<1e-10){
302  //no pandora t0 found, correct for t0
303  double ticksOffset = 0;
304  if (this_t0crt2 > -DBL_MAX) ticksOffset = this_t0crt2/500.+detProp.GetXTicksOffset(allHits[0]->WireID());
305  else if (this_t0crt1 > -DBL_MAX) ticksOffset = this_t0crt1/500.+detProp.GetXTicksOffset(allHits[0]->WireID());
306  else if (this_t0anodep > -DBL_MAX) ticksOffset = this_t0anodep/500.+detProp.GetXTicksOffset(allHits[0]->WireID());
307  double xOffset = detProp.ConvertTicksToX(ticksOffset,allHits[0]->WireID());
308  this_trackstartx -= xOffset;
309  this_trackendx -= xOffset;
310  }
311 
312  auto const & posOffsets = SCE->GetCalPosOffsets(geo::Point_t(this_trackstartx, this_trackstarty, this_trackstartz), allHits[0]->WireID().TPC);
313  this_trackstartx_sce = this_trackstartx - posOffsets.X();
314  this_trackstarty_sce = this_trackstarty + posOffsets.Y();
315  this_trackstartz_sce = this_trackstartz + posOffsets.Z();
316  this_trackendx_sce = this_trackendx - posOffsets.X();
317  this_trackendy_sce = this_trackendy + posOffsets.Y();
318  this_trackendz_sce = this_trackendz + posOffsets.Z();
319 
320  trackid.push_back(this_trackid);
321  t0crt2.push_back(this_t0crt2);
322  t0crt1.push_back(this_t0crt1);
323  t0pandora.push_back(this_t0pandora);
324  t0anodep.push_back(this_t0anodep);
325  t0truth.push_back(this_t0truth);
326  trackstartx.push_back(this_trackstartx);
327  trackstarty.push_back(this_trackstarty);
328  trackstartz.push_back(this_trackstartz);
329  trackendx.push_back(this_trackendx);
330  trackendy.push_back(this_trackendy);
331  trackendz.push_back(this_trackendz);
332  trackstartx_sce.push_back(this_trackstartx_sce);
333  trackstarty_sce.push_back(this_trackstarty_sce);
334  trackstartz_sce.push_back(this_trackstartz_sce);
335  trackendx_sce.push_back(this_trackendx_sce);
336  trackendy_sce.push_back(this_trackendy_sce);
337  trackendz_sce.push_back(this_trackendz_sce);
338  crt1x.push_back(this_crt1x);
339  crt1y.push_back(this_crt1y);
340  crt1z.push_back(this_crt1z);
341  crt2x0.push_back(this_crt2x0);
342  crt2y0.push_back(this_crt2y0);
343  crt2z0.push_back(this_crt2z0);
344  crt2x1.push_back(this_crt2x1);
345  crt2y1.push_back(this_crt2y1);
346  crt2z1.push_back(this_crt2z1);
347  }
348  }
349 
350  //RDTimeStamps for the event
351  std::vector < art::Ptr < raw::RDTimeStamp > > rdtsevtList;
352  art::InputTag itag("timingrawdecoder","daq");
353  auto rdts_evt = e.getHandle< std::vector < raw::RDTimeStamp > >(itag);
354  if (rdts_evt) {
355  art::fill_ptr_vector(rdtsevtList, rdts_evt);
356  }
357 
358  if (!rdtsevtList.empty()){
359  rdtimestamp_evt = rdtsevtList[0]->GetTimeStamp();
360  }
361 
362  //RDTimeStamps for each raw digit (channel)
363  art::InputTag itag2("tpcrawdecoder","daq");
364  std::vector < art::Ptr < raw::RDTimeStamp > > rdtsdigitList;
365  auto rdts_digit = e.getHandle < std::vector < raw::RDTimeStamp > >(itag2);
366  if (rdts_digit) {
367  art::fill_ptr_vector(rdtsdigitList, rdts_digit);
368  }
369 
370  for (const auto & rdts : rdtsdigitList){
371  rdtimestamp_digits.push_back(rdts->GetTimeStamp());
372  }
373 
374  if (!trackid.empty()) t0tree->Fill();
375 }
376 
378  art::ServiceHandle<art::TFileService> fileServiceHandle;
379  t0tree = fileServiceHandle->make<TTree>("t0tree", "t0 info");
380  t0tree->Branch("run", &run, "run/I");
381  t0tree->Branch("event", &event, "event/I");
382  t0tree->Branch("trackid", &trackid);
383  t0tree->Branch("t0crt2", &t0crt2);
384  t0tree->Branch("t0crt1", &t0crt1);
385  t0tree->Branch("t0pandora",&t0pandora);
386  t0tree->Branch("t0anodep",&t0anodep);
387  t0tree->Branch("t0truth", &t0truth);
388  t0tree->Branch("trackstartx",&trackstartx);
389  t0tree->Branch("trackstarty",&trackstarty);
390  t0tree->Branch("trackstartz",&trackstartz);
391  t0tree->Branch("trackendx",&trackendx);
392  t0tree->Branch("trackendy",&trackendy);
393  t0tree->Branch("trackendz",&trackendz);
394  t0tree->Branch("trackstartx_sce",&trackstartx_sce);
395  t0tree->Branch("trackstarty_sce",&trackstarty_sce);
396  t0tree->Branch("trackstartz_sce",&trackstartz_sce);
397  t0tree->Branch("trackendx_sce",&trackendx_sce);
398  t0tree->Branch("trackendy_sce",&trackendy_sce);
399  t0tree->Branch("trackendz_sce",&trackendz_sce);
400  t0tree->Branch("crt1x",&crt1x);
401  t0tree->Branch("crt1y",&crt1y);
402  t0tree->Branch("crt1z",&crt1z);
403  t0tree->Branch("crt2x0",&crt2x0);
404  t0tree->Branch("crt2y0",&crt2y0);
405  t0tree->Branch("crt2z0",&crt2z0);
406  t0tree->Branch("crt2x1",&crt2x1);
407  t0tree->Branch("crt2y1",&crt2y1);
408  t0tree->Branch("crt2z1",&crt2z1);
409  t0tree->Branch("rdtimestamp_evt", &rdtimestamp_evt, "rdtimestamp_evt/D");
410  t0tree->Branch("rdtimestamp_digits", &rdtimestamp_digits);
411 }
intermediate_table::iterator iterator
std::vector< double > crt2z1
std::vector< double > crt1y
std::vector< double > trackstartx
std::vector< double > trackendy_sce
CheckT0 & operator=(CheckT0 const &)=delete
std::vector< double > crt2y0
std::vector< sim::TrackIDE > HitToTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
std::vector< double > t0crt2
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
std::vector< TrackID > TrackIDs
const simb::MCParticle * TrackIdToParticle_P(int id) const
void beginJob() override
std::vector< double > crt2y1
std::vector< double > trackendx
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
std::vector< double > crt1x
CheckT0(fhicl::ParameterSet const &p)
std::vector< double > trackendy
std::vector< double > rdtimestamp_digits
std::vector< double > crt1z
bool isRealData() const
std::vector< double > trackstarty_sce
std::vector< double > crt2z0
T abs(T value)
std::vector< double > t0crt1
const double e
double rdtimestamp_evt
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
def key(type, name=None)
Definition: graph.py:13
std::vector< double > trackendz
void analyze(art::Event const &e) override
double T(const int i=0) const
Definition: MCParticle.h:224
std::vector< double > trackendx_sce
p
Definition: test.py:223
std::vector< double > trackendz_sce
RunNumber_t run() const
Definition: DataViewImpl.cc:71
std::vector< double > trackstartz
Detector simulation of raw signals on wires.
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:184
std::vector< double > trackstartx_sce
std::vector< int > trackid
Declaration of signal hit object.
std::vector< double > t0truth
std::vector< double > t0pandora
std::vector< double > trackstartz_sce
Provides recob::Track data product.
EventNumber_t event() const
Definition: EventID.h:116
std::vector< double > crt2x1
std::vector< double > trackstarty
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
std::vector< double > t0anodep
EventID id() const
Definition: Event.cc:34
Event finding and building.
std::vector< double > crt2x0