PixelMapProducer.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file PixelMapProducer.h
3 /// \brief PixelMapProducer for CVN
4 /// \author Alexander Radovic - a.radovic@gmail.com
5 //
6 // Modifications to allow unwrapped collection view
7 // - Leigh Whitehead - leigh.howard.whitehead@cern.ch
8 ////////////////////////////////////////////////////////////////////////
9 
10 #include <iostream>
11 #include <ostream>
12 #include <list>
13 #include <algorithm>
14 #include <numeric>
15 
18 #include "TVector2.h"
19 #include "TH2D.h"
22 
28 
29 namespace cvn
30 {
31 
32  PixelMapProducer::PixelMapProducer(unsigned int nWire, unsigned int nTdc, double tRes):
33  fNWire(nWire),
34  fNTdc(nTdc),
35  fTRes(tRes),
36  fUnwrapped(2),
37  fProtoDUNE(false)
38  {
39 
41  if (fGeometry->DetectorName().find("dunevd10kt_3view") != std::string::npos)
43  }
44 
46  {
48  if (fGeometry->DetectorName().find("dunevd10kt_3view") != std::string::npos)
50  }
51 
54  {
55  std::vector<const recob::Hit*> newCluster;
56  for(const art::Ptr<recob::Hit> hit : cluster){
57  newCluster.push_back(hit.get());
58  }
59  return CreateMap(detProp, newCluster);
60  }
61 
63  const std::vector<const recob::Hit* >& cluster)
64  {
65  if (fGeometry->DetectorName().find("dunevd10kt_3view") != std::string::npos){
66  fVDPlane0.clear();
67  fVDPlane1.clear();
69  }
70 
71  Boundary bound = DefineBoundary(detProp, cluster);
72  return CreateMapGivenBoundary(detProp, cluster, bound);
73  }
74 
76  const std::vector<const recob::Hit*>& cluster,
77  const Boundary& bound)
78  {
79 
80  PixelMap pm(fNWire, fNTdc, bound);
81 
82  for(size_t iHit = 0; iHit < cluster.size(); ++iHit)
83  {
84 
85  geo::WireID wireid = cluster[iHit]->WireID();
86  double temptdc = cluster[iHit]->PeakTime();
87  unsigned int tempWire = wireid.Wire;
88  unsigned int tempPlane = wireid.Plane;
89  // Leigh: Simple modification to unwrap the collection view wire plane
90  if(!fProtoDUNE){
91  if(fUnwrapped == 1){
92  // Jeremy: Autodetect geometry for DUNE 10kt module. Is this a bad idea??
93  if (fGeometry->DetectorName() == "dune10kt_v1") {
94  if (wireid.TPC%6 == 0 or wireid.TPC%6 == 5) continue; // Skip dummy TPCs in 10kt module
95  GetDUNE10ktGlobalWireTDC(detProp, wireid.Wire,cluster[iHit]->PeakTime(),
96  wireid.Plane,wireid.TPC,tempWire,tempPlane,temptdc);
97  }
98  else if (fGeometry->DetectorName().find("dunevd10kt_3view") != std::string::npos){
99  GetDUNEVertDrift3ViewGlobalWire(wireid.Wire, wireid.Plane,wireid.TPC,tempWire,tempPlane);
100  }
101  // Default to 1x2x6. Should probably specifically name this function as such
102  else {
103  GetDUNEGlobalWireTDC(detProp, wireid.Wire,cluster[iHit]->PeakTime(),
104  wireid.Plane,wireid.TPC,tempWire,tempPlane,temptdc);
105  }
106  }
107  else if(fUnwrapped == 2){
108  // Old method that has problems with the APA crossers, kept for old times' sake
109  GetDUNEGlobalWire(wireid.Wire,wireid.Plane,wireid.TPC,tempWire,tempPlane);
110  }
111  }
112  else{
113  GetProtoDUNEGlobalWire(wireid.Wire,wireid.Plane,wireid.TPC,tempWire,tempPlane);
114  }
115  const double pe = cluster[iHit]->Integral();
116  const unsigned int wire = tempWire;
117  const unsigned int wirePlane = tempPlane;
118  const double tdc = temptdc;
119  pm.Add(wire, tdc, wirePlane, pe);
120 
121  }
122  return pm;
123  }
124 
125  std::ostream& operator<<(std::ostream& os, const PixelMapProducer& p)
126  {
127  os << "PixelMapProducer: "
128  << p.NTdc() <<" tdcs X " << p.NWire() << " wires";
129  return os;
130  }
131 
133  {
134  const geo::WireGeo* pwire = fGeometry->WirePtr(wireid);
136  double slope = 0.;
137  if(!pwire->isVertical()) slope = pwire->TanThetaZ();
138 
139  double intercept = center.Y() - slope*center.Z();
140  if(wireid.Plane == 2) intercept = 0.;
141 
142  return intercept;
143  }
144 
146 
147  // double spacing = 0.847;
148  for(int plane = 0; plane < 2; plane++){
149 
150  int nCRM_row = 6;
151  for(int diag_tpc = 0; diag_tpc < nCRM_row; diag_tpc++){
152 
153  unsigned int nWiresTPC = fGeometry->Nwires(plane, 0, 0);
154  int tpc_id = plane == 0 ? (nCRM_row+1)*diag_tpc : (nCRM_row-1)*(nCRM_row-diag_tpc);
155  geo::WireID start = geo::WireID(0, tpc_id, plane, 0);
156  geo::WireID end = geo::WireID(0, tpc_id, plane, nWiresTPC-1);
157 
158  double start_intercept = _getIntercept(start);
159  double end_intercept = _getIntercept(end);
160  if(plane == 0){
161  fVDPlane0.push_back(start_intercept);
162  fVDPlane0.push_back(end_intercept);
163  }
164  else{
165  fVDPlane1.push_back(_getIntercept(end));
166  fVDPlane1.push_back(_getIntercept(start));
167  }
168  }
169 
170  }
171  }
172 
174  const std::vector< const recob::Hit*>& cluster)
175  {
176 
177  std::vector<double> time_0;
178  std::vector<double> time_1;
179  std::vector<double> time_2;
180  //std::vector<int> wire;
181 
182  std::vector<int> wire_0;
183  std::vector<int> wire_1;
184  std::vector<int> wire_2;
185 
186 
187  for(size_t iHit = 0; iHit < cluster.size(); ++iHit)
188  {
189  geo::WireID wireid = cluster[iHit]->WireID();
190 
191  unsigned int globalWire = wireid.Wire;
192  unsigned int globalPlane = wireid.Plane;
193  double globalTime = cluster[iHit]->PeakTime();
194  if(!fProtoDUNE){
195  if(fUnwrapped == 1) {
196  if (fGeometry->DetectorName() == "dune10kt_v1") {
197  if (wireid.TPC%6 == 0 or wireid.TPC%6 == 5) continue; // Skip dummy TPCs in 10kt module
198  GetDUNE10ktGlobalWireTDC(detProp, wireid.Wire,cluster[iHit]->PeakTime(),wireid.Plane,wireid.TPC,globalWire,globalPlane,globalTime);
199  }
200  else if (fGeometry->DetectorName().find("dunevd10kt_3view") != std::string::npos){
201  GetDUNEVertDrift3ViewGlobalWire(wireid.Wire, wireid.Plane,wireid.TPC,globalWire,globalPlane);
202  }
203  else {
204  GetDUNEGlobalWireTDC(detProp, wireid.Wire,cluster[iHit]->PeakTime(),wireid.Plane,wireid.TPC,globalWire,globalPlane,globalTime);
205  }
206  }
207  else if(fUnwrapped == 2){
208  GetDUNEGlobalWire(wireid.Wire,wireid.Plane,wireid.TPC,globalWire,globalPlane);
209  }
210  }
211  else{
212  GetProtoDUNEGlobalWire(wireid.Wire,wireid.Plane,wireid.TPC,globalWire,globalPlane);
213  }
214 
215  if(globalPlane==0){
216  time_0.push_back(globalTime);
217  wire_0.push_back(globalWire);
218  }
219  if(globalPlane==1){
220  time_1.push_back(globalTime);
221  wire_1.push_back(globalWire);
222  }
223  if(globalPlane==2){
224  time_2.push_back(globalTime);
225  wire_2.push_back(globalWire);
226  }
227  }
228 
229  double tsum_0 = std::accumulate(time_0.begin(), time_0.end(), 0.0);
230  double tmean_0 = tsum_0 / time_0.size();
231 
232  double tsum_1 = std::accumulate(time_1.begin(), time_1.end(), 0.0);
233  double tmean_1 = tsum_1 / time_1.size();
234 
235  double tsum_2 = std::accumulate(time_2.begin(), time_2.end(), 0.0);
236  double tmean_2 = tsum_2 / time_2.size();
237 
238  std::cout << "Boundary wire vector sizes: " << wire_0.size() << ", " << wire_1.size() << ", " << wire_2.size() << std::endl;
239 
240  auto minwireelement_0= std::min_element(wire_0.begin(), wire_0.end());
241  std::cout<<"minwire 0: "<<*minwireelement_0<<std::endl;
242  auto minwireelement_1= std::min_element(wire_1.begin(), wire_1.end());
243  std::cout<<"minwire 1: "<<*minwireelement_1<<std::endl;
244  auto minwireelement_2= std::min_element(wire_2.begin(), wire_2.end());
245  std::cout<<"minwire 2: "<<*minwireelement_2<<std::endl;
246 
247  //auto maxwireelement_0= std::max_element(wire_0.begin(), wire_0.end());
248  //std::cout<<"maxwire 0: "<<*maxwireelement_0<<std::endl;
249  //auto maxwireelement_1= std::max_element(wire_1.begin(), wire_1.end());
250  //std::cout<<"maxwire 1: "<<*maxwireelement_1<<std::endl;
251  //auto maxwireelement_2= std::max_element(wire_2.begin(), wire_2.end());
252  //std::cout<<"maxwire 2: "<<*maxwireelement_2<<std::endl;
253 
254 
255  int minwire_0 = 0;
256  int minwire_1 = 0;
257  int minwire_2 = 0;
258  if(wire_0.size() > 0) minwire_0 = *minwireelement_0-1;
259  if(wire_1.size() > 0) minwire_1 = *minwireelement_1-1;
260  if(wire_2.size() > 0) minwire_2 = *minwireelement_2-1;
261 
262  Boundary bound(fNWire,fTRes,minwire_0,minwire_1,minwire_2,tmean_0,tmean_1,tmean_2);
263 
264  return bound;
265  }
266 
267  void PixelMapProducer::GetDUNEGlobalWire(unsigned int localWire, unsigned int plane, unsigned int tpc, unsigned int& globalWire, unsigned int& globalPlane) const
268  {
269  unsigned int nWiresTPC = 400;
270 
271  globalWire = localWire;
272  globalPlane = 0;
273 
274  // Collection plane has more wires
275  if(plane == 2){
276  nWiresTPC=480;
277  globalPlane = 2;
278  }
279 
280  // Workspace geometry has two drift regions
281  // |-----|-----| / /
282  // y ^ | 3 | 2 |/ /
283  // | -| z |-----|-----| /
284  // | / | 1 | 0 | /
285  // x <---|/ |-----|-----|/
286  //
287 
288  int tpcMod4 = tpc%4;
289 
290  // Induction views depend on the drift direction
291  if(plane < 2){
292  // For drift in negative x direction keep U and V as defined.
293  if(tpcMod4 == 0 || tpcMod4 == 3){
294  globalPlane = plane;
295  }
296  // For drift in positive x direction, swap U and V.
297  else{
298  if(plane == 0) globalPlane = 1;
299  else globalPlane = 0;
300  }
301  }
302 
303  if(globalPlane != 1){
304  globalWire += (tpc/4)*nWiresTPC;
305  }
306  else{
307  globalWire += ((23-tpc)/4)*nWiresTPC;
308  }
309 
310  }
311 
312  // Based on Robert's code in adcutils
314  unsigned int localWire, double localTDC, unsigned int plane, unsigned int tpc,
315  unsigned int& globalWire, unsigned int& globalPlane, double& globalTDC) const
316  {
317 
318  unsigned int nWiresTPC = 400;
319  unsigned int wireGap = 4;
320  double driftLen = fGeometry->TPC(tpc,0).DriftDistance();
321  double apaLen = fGeometry->TPC(tpc,0).Width() - fGeometry->TPC(tpc,0).ActiveWidth();
322  double driftVel = detProp.DriftVelocity();
323  unsigned int drift_size = (driftLen / driftVel) * 2; // Time in ticks to cross a TPC
324  unsigned int apa_size = 4*(apaLen / driftVel) * 2; // Width of the whole APA in TDC
325 
326  globalWire = 0;
327  globalPlane = 0;
328 // int dir = fGeometry->TPC(tpc,0).DetectDriftDirection();
329 
330  // Collection plane has more wires
331  if(plane == 2){
332  nWiresTPC = 480;
333  wireGap = 5;
334  globalPlane = 2;
335  }
336 
337  bool includeZGap = true;
338  if(includeZGap) nWiresTPC += wireGap;
339 
340  // Workspace geometry has two drift regions
341  // |-----|-----| / /
342  // y ^ | 3 | 2 |/ /
343  // | -| z |-----|-----| /
344  // | / | 1 | 0 | /
345  // x <---|/ |-----|-----|/
346  //
347  int tpcMod4 = tpc%4;
348  // Induction views depend on the drift direction
349  if (plane < 2 and tpc%2 == 1) globalPlane = !plane;
350  else globalPlane = plane;
351 
352  int offset = 752; // Offset between upper and lower modules in induction views, from Robert & Dorota's code
353  // Second induction plane gets offset from the back of the TPC
354  if (globalPlane != 1) globalWire += (tpc/4)*nWiresTPC;
355  else globalWire += ((23-tpc)/4)*nWiresTPC;
356  // Reverse wires and add offset for upper modules in induction views
357  // Nitish : what's the difference between Nwires here and nWiresTPC?
358  if (tpcMod4 > 1 and globalPlane < 2) globalWire += fGeometry->Nwires(globalPlane, tpc, 0) + offset - localWire;
359  else globalWire += localWire;
360 
361  if(tpcMod4 == 0 || tpcMod4 == 2){
362  globalTDC = drift_size - localTDC;
363  }
364  else{
365  globalTDC = localTDC + drift_size + apa_size;
366  }
367  }
368 
370  unsigned int localWire, double localTDC, unsigned int plane, unsigned int tpc,
371  unsigned int& globalWire, unsigned int& globalPlane, double& globalTDC) const
372  {
373  unsigned int nWiresTPC = 400;
374  unsigned int wireGap = 4;
375  double driftLen = fGeometry->TPC(tpc).DriftDistance();
376  double apaLen = fGeometry->TPC(tpc).Width() - fGeometry->TPC(tpc).ActiveWidth();
377  double driftVel = detProp.DriftVelocity();
378  unsigned int drift_size = (driftLen / driftVel) * 2; // Time in ticks to cross a TPC
379  unsigned int apa_size = 4*(apaLen / driftVel) * 2; // Width of the whole APA in TDC
380 
381 
382  globalWire = 0;
383  globalPlane = 0;
384 
385  // Collection plane has more wires
386  if(plane == 2){
387  nWiresTPC = 480;
388  wireGap = 5;
389  globalPlane = 2;
390  }
391 
392  bool includeZGap = true;
393  if(includeZGap) nWiresTPC += wireGap;
394 
395  // 10kt has four real TPCs and two dummies in each slice
396  //
397  // |--|-----|-----|-----|-----|--| / /
398  // y ^ |11| 10 | 9 | 8 | 7 | 6|/ /
399  // | -| z |--|-----|-----|-----|-----|--| /
400  // | / | 5| 4 | 3 | 2 | 1 | 0| /
401  // x <---|/ |--|-----|-----|-----|-----|--|/
402  // ^ wires ^ ^ wires ^
403  //
404  // We already filtered out the dummies, so we can assume 0->3 as follows:
405  //
406  // |-----|-----|-----|-----| / /
407  // y ^ | 7 | 6 | 5 | 4 |/ /
408  // | -| z |-----|-----|-----|-----| /
409  // | / | 3 | 2 | 1 | 0 | /
410  // x <---|/ |-----|-----|-----|-----|/
411  // ^ wires ^ ^ wires ^
412  //
413 
414  size_t tpc_x = (tpc%6) - 1; // x coordinate in 0->4 range
415  size_t tpc_xy = (tpc%12) - 1; // xy coordinate as 0->3 & 6->9 (converted from 1->4, 7->10)
416  if (tpc_xy > 3) tpc_xy -= 2; // now subtract 2 so it's in the 0->7 range
417 
418  // Induction views depend on the drift direction
419  if (plane < 2 and tpc%2 == 1) globalPlane = !plane;
420  else globalPlane = plane;
421 
422  int offset = 752; // Offset between upper and lower modules in induction views, from Robert & Dorota's code
423  // Second induction plane gets offset from the back of the TPC
424  if (globalPlane != 1) globalWire += (tpc/12)*nWiresTPC;
425  else globalWire += ((300-tpc)/12)*nWiresTPC;
426  // Reverse wires and add offset for upper modules in induction views
427  if (tpc_xy > 3 and globalPlane < 2) globalWire += fGeometry->Nwires(globalPlane, tpc, 0) + offset - localWire;
428  else globalWire += localWire;
429 
430  if (tpc_x % 2 == 0) globalTDC = localTDC;
431  else globalTDC = (2*drift_size) - localTDC;
432  if (tpc_x > 1) globalTDC += 2 * (drift_size + apa_size);
433 
434  } // function PixelMapProducer::GetDUNE10ktGlobalWireTDC
435 
436  // Special case for ProtoDUNE where we want to extract single particles to mimic CCQE interactions. The output pixel maps should be the same as the workspace
437  // but we need different treatment of the wire numbering
438  void PixelMapProducer::GetProtoDUNEGlobalWire(unsigned int localWire, unsigned int plane, unsigned int tpc, unsigned int& globalWire, unsigned int& globalPlane) const
439  {
440  unsigned int nWiresTPC = 400;
441 
442  globalWire = localWire;
443  globalPlane = 0;
444 
445  // Collection plane has more wires
446  if(plane == 2){
447  nWiresTPC=480;
448  globalPlane = 2;
449  }
450 
451  // ProtoDUNE has a central CPA so times are fine
452  // It (annoyingly) has two dummy TPCs on the sides
453  //
454  // y ^ |-|-----|-----|-| /
455  // | -| z | | | | | /
456  // | / |3| 2 | 1 |0| /
457  // x <---|/ |-|-----|-----|-|/
458  //
459 
460  int tpcMod4 = tpc%4;
461  // tpcMod4: 1 for -ve drift, 2 for +ve drift
462  // Induction views depend on the drift direction
463  if(plane < 2){
464  // For drift in negative x direction keep U and V as defined.
465  if(tpcMod4 == 1){
466  globalPlane = plane;
467  }
468  // For drift in positive x direction, swap U and V.
469  else{
470  if(plane == 0) globalPlane = 1;
471  else globalPlane = 0;
472  }
473  }
474 
475  if(globalPlane != 1){
476  globalWire += (tpc/4)*nWiresTPC;
477  }
478  else{
479  globalWire += ((12-tpc)/4)*nWiresTPC;
480  }
481 
482  } // function PixelMapProducer::GetProtoDUNEGlobalWire
483 
484  // Special case for ProtoDUNE where we want to extract single particles to mimic CCQE interactions. The output pixel maps should be the same as the workspace
485  // but we need different treatment of the wire numbering
486  void PixelMapProducer::GetProtoDUNEGlobalWireTDC(unsigned int localWire, double localTDC, unsigned int plane, unsigned int tpc,
487  unsigned int& globalWire, double& globalTDC, unsigned int& globalPlane) const
488  {
489  // We can just use the existing function to get the global wire & plane
490  // GetProtoDUNEGlobalWire(localWire, plane, tpc, globalWire, globalPlane);
491  GetDUNEGlobalWire(localWire, plane, tpc, globalWire, globalPlane);
492  // Implement additional mirroring here?
493 
494  } // function GetProtoDUNEGlobalWireTDC
495 
496  void PixelMapProducer::GetDUNEVertDrift3ViewGlobalWire(unsigned int localWire, unsigned int plane, unsigned int tpc, unsigned int& globalWire, unsigned int& globalPlane) const
497  {
498  // Preliminary function for VD Geometries (3 View for now)
499  // 1x6x6 -- single drift volume should make things significantly simpler
500 
501  int nCRM_row = 6;
502  // spacing between y-intercepts of parallel wires in a given plane.
503  // seems its not actually the pitch/cos(theta_z) but rather the pitch/cos(theta_z) - 2*r_wire ??
504  double spacing = 0.847;
505 
506  globalPlane = plane;
507  unsigned int nWiresTPC = fGeometry->Nwires(globalPlane, tpc, 0);
508 
509  if(globalPlane < 2){
510 
511  geo::WireID wire_id = geo::WireID(0, tpc, globalPlane, localWire);
512  double wire_intercept = _getIntercept(wire_id);
513  double low_bound, upper_bound;
514  int start, end, diag_tpc;
515  // double matched_intercept;
516  // get wires on diagonal CRMs and their intercepts which bound the current wire's intercept
517  if(globalPlane == 0){
518  start = std::lower_bound(fVDPlane0.begin(), fVDPlane0.end(), wire_intercept) - fVDPlane0.begin() - 1;
519  end = std::upper_bound(fVDPlane0.begin(), fVDPlane0.end(), wire_intercept) - fVDPlane0.begin();
520  low_bound = fVDPlane0[start];
521  upper_bound = fVDPlane0[end];
522  diag_tpc = (start/2);
523  }
524  else{
525  end = std::lower_bound(fVDPlane1.begin(), fVDPlane1.end(), wire_intercept) - fVDPlane1.begin() - 1;
526  start = std::upper_bound(fVDPlane1.begin(), fVDPlane1.end(), wire_intercept) - fVDPlane1.begin();
527  low_bound = fVDPlane1[end];
528  upper_bound = fVDPlane1[start];
529  diag_tpc = (nCRM_row-(end/2) - 1);
530  }
531  // if the intercept of the wire is in between two diagonal CRMs, assign it to the diagonal CRM its closest to
532  if((start % 2)^globalPlane){
533 
534  int diag_idx = diag_tpc + !globalPlane;
535  globalWire = (wire_intercept > (low_bound+upper_bound)*0.5) ? (nWiresTPC-1)*diag_idx + !globalPlane : (nWiresTPC-1)*diag_idx + globalPlane;
536  }
537  // otherwise assign it to the closest wire within the same CRM
538  else{
539  int diag_idx = diag_tpc;
540  int offset = globalPlane ? std::round((upper_bound - wire_intercept)/spacing) : std::round((wire_intercept-low_bound)/spacing);
541  globalWire = (nWiresTPC-1)*diag_idx + offset + 1;
542 
543  }
544  }
545  else{
546  int tpc_z = tpc/6;
547  globalWire = localWire + tpc_z*nWiresTPC;
548  }
549 
550  }
551 
553  art::Ptr<recob::Hit>& hit, std::vector<int>& pdgs,
554  std::vector<int>& tracks, std::vector<float>& energy, std::vector<std::string>& process) {
555 
556  // BackTracker and ParticleInventory
559 
560  // Get true particle and PDG responsible for this hit
561  std::vector<sim::TrackIDE> IDEs = bt->HitToTrackIDEs(clockData, hit);
562  for (sim::TrackIDE & k : IDEs) {
563  tracks.push_back(k.trackID); // add track ID
564  simb::MCParticle p = pi->TrackIdToParticle(k.trackID);
565  process.push_back(p.Process()); // add G4 process string
566 
567  // Manually check to see if we have a Michel electron here
568  int pdg = p.PdgCode();
569  if (abs(pdg) == 11 && p.Mother() != 0) {
570  int pdgParent = pi->TrackIdToParticle(p.Mother()).PdgCode();
571  if (abs(pdgParent) == 13 && ((pdg>0)==(pdgParent>0))) {
572  // If it's a Michel electron, tag it with an unphysical PDG code
573  if (p.Process() == "muMinusCaptureAtRest" || p.Process() == "muPlusCaptureAtRest") {
574  pdg = 99;
575  }
576  } // If electron parent is muon
577  } // If non-primary electron
578  pdgs.push_back(pdg);
579  energy.push_back(k.energy);
580  } // for trackIDE
581  } // function PixelMapProducer::GetHitTruth
582 
584  detinfo::DetectorClocksData const& clockData,
585  detinfo::DetectorPropertiesData const& detProp,
587  bool usePixelTruth) {
588 
589  // 3-dimensional coordinates (wire, time TPC) and 3 views
590  SparsePixelMap map(3, 3, usePixelTruth);
591 
594 
595  for(size_t iHit = 0; iHit < cluster.size(); ++iHit) {
596 
597  geo::WireID wireid = cluster[iHit]->WireID();
598  double globalTime = cluster[iHit]->PeakTime();
599  unsigned int globalWire = wireid.Wire;
600  unsigned int globalPlane = wireid.Plane;
601 
602  if (fGeometry->DetectorName().find("1x2x6") != std::string::npos) {
603  GetDUNEGlobalWireTDC(detProp, wireid.Wire, cluster[iHit]->PeakTime(),
604  wireid.Plane, wireid.TPC, globalWire, globalPlane, globalTime);
605  }
606  else if (fGeometry->DetectorName() == "dune10kt_v1") {
607  if (wireid.TPC%6 == 0 or wireid.TPC%6 == 5) continue;
608  GetDUNE10ktGlobalWireTDC(detProp, wireid.Wire, cluster[iHit]->PeakTime(),
609  wireid.Plane, wireid.TPC, globalWire, globalPlane, globalTime);
610  }
611  else if (fGeometry->DetectorName().find("protodune") != std::string::npos) {
612  GetProtoDUNEGlobalWireTDC(wireid.Wire, cluster[iHit]->PeakTime(),
613  wireid.Plane, wireid.TPC, globalWire, globalTime, globalPlane);
614  }
616  << "Geometry " << fGeometry->DetectorName() << " not implemented "
617  << "in CreateSparseMap." << std::endl;
618 
619  std::vector<float> coordinates = { (float)globalWire, (float)globalTime, (float)wireid.TPC };
620 
621  if (usePixelTruth) {
622  // Get true information for this hit
623  std::vector<int> pdgs, tracks;
624  std::vector<float> energy;
625  std::vector<std::string> process;
626  GetHitTruth(clockData, cluster[iHit], pdgs, tracks, energy, process);
627  map.AddHit(globalPlane, coordinates, {cluster[iHit]->Integral()}, pdgs, tracks, energy, process);
628  } // if PixelTuth
629 
630  else {
631  map.AddHit(0, coordinates, {cluster[iHit]->Integral()});
632  }
633  } // for iHit
634 
635  return map;
636 
637  } // function PixelMapProducer::CreateSparseMap
638 
640  detinfo::DetectorClocksData const& clockData,
643 
644  // 3D coordinates (x,y,z) and a single 3D view
645  SparsePixelMap map(3, 1, true);
646 
649 
650  for (size_t iSP = 0; iSP < sp.size(); ++iSP) { // Loop over spacepoints
651 
652  std::vector<float> features(3, 0.); // charge on each plane
653  std::vector<float> coordinates(3);
654  const double *pos = sp[iSP]->XYZ();
655  for (size_t p = 0; p < 3; ++p) coordinates[p] = pos[p];
656 
657  std::vector<int> pdgs, tracks;
658  std::vector<float> energy;
659  std::vector<std::string> process;
660 
661  for (size_t iH = 0; iH < hit[iSP].size(); ++iH) { // Loop over this spacepoint's hits
662  features[hit[iSP][iH]->View()] += hit[iSP][iH]->Integral(); // Add hit integral to corresponding view's features
663  GetHitTruth(clockData, hit[iSP][iH], pdgs, tracks, energy, process);
664  map.AddHit(0, coordinates, features, pdgs, tracks, energy, process);
665  } // for hit iH
666  } // for spacepoint iSP
667 
668  return map;
669 
670  } // function PixelMapProducer::CreateSparseMap
671 
672 } // namespace cvn
double fTRes
Timing resolution for pixel map.
bool isVertical() const
Returns if this wire is vertical (theta_z ~ pi/2)
Definition: WireGeo.h:276
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:65
coordinates
Definition: geo_types.h:121
int PdgCode() const
Definition: MCParticle.h:212
void AddHit(unsigned int view, std::vector< float > coordinates, std::vector< float > features)
Default AddHit implementation, which just adds pixel value and coordinates.
std::vector< sim::TrackIDE > HitToTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
unsigned short fUnwrapped
Use unwrapped pixel maps?
int Mother() const
Definition: MCParticle.h:213
unsigned int NTdc() const
std::ostream & operator<<(std::ostream &os, const PixelMapProducer &p)
struct vector vector
double globalTime
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
Cluster finding and building.
std::string Process() const
Definition: MCParticle.h:215
Utility class for truth labels.
unsigned int fNTdc
Number of tdcs, width of pixel map.
double Width() const
Width is associated with x coordinate [cm].
Definition: TPCGeo.h:109
void GetDUNE10ktGlobalWireTDC(detinfo::DetectorPropertiesData const &detProp, unsigned int localWire, double localTDC, unsigned int plane, unsigned int tpc, unsigned int &globalWire, unsigned int &globalPlane, double &globalTDC) const
unsigned int Nwires(unsigned int p, unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wires in the specified plane.
const char features[]
Definition: feature_tests.c:2
art framework interface to geometry description
PixelMap CreateMapGivenBoundary(detinfo::DetectorPropertiesData const &detProp, const std::vector< const recob::Hit * > &cluster, const Boundary &bound)
void GetProtoDUNEGlobalWire(unsigned int localWire, unsigned int plane, unsigned int tpc, unsigned int &globalWire, unsigned int &globalPlane) const
def process(f, kind)
Definition: search.py:254
T abs(T value)
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -> PDG code
simb::MCParticle TrackIdToParticle(int const id) const
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
bt
Definition: tracks.py:83
double _getIntercept(geo::WireID wireid) const
std::string DetectorName() const
Returns a string with the name of the detector, as configured.
void Add(const unsigned int &wire, const double &tdc, const unsigned int &view, const double &pe)
Definition: PixelMap.cxx:46
bool fProtoDUNE
Do we want to use this for particle extraction from protoDUNE?
p
Definition: test.py:223
double DriftVelocity(double efield=0., double temperature=0.) const
cm/us
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
void GetProtoDUNEGlobalWireTDC(unsigned int localWire, double localTDC, unsigned int plane, unsigned int tpc, unsigned int &globalWire, double &globalTDC, unsigned int &globalPlane) const
Detector simulation of raw signals on wires.
std::vector< double > fVDPlane1
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
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
Definition: tracks.py:1
double TanThetaZ() const
Definition: WireGeo.h:266
unsigned int fNWire
Number of wires, length for pixel maps.
Declaration of signal hit object.
double DriftDistance() const
Definition: TPCGeo.h:155
void GetDUNEGlobalWireTDC(detinfo::DetectorPropertiesData const &detProp, unsigned int localWire, double localTDC, unsigned int plane, unsigned int tpc, unsigned int &globalWire, unsigned int &globalPlane, double &globalTDC) const
Contains all timing reference information for the detector.
unsigned int NWire() const
SparsePixelMap CreateSparseMap2D(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit > > &cluster, bool usePixelTruth=false)
double ActiveWidth() const
Width (associated with x coordinate) of active TPC volume [cm].
Definition: TPCGeo.h:97
PixelMap, basic input to CVN neural net.
Definition: PixelMap.h:22
def center(depos, point)
Definition: depos.py:117
PixelMapProducer for CVN.
void GetDUNEGlobalWire(unsigned int localWire, unsigned int plane, unsigned int tpc, unsigned int &globalWire, unsigned int &globalPlane) const
Function to convert to a global unwrapped wire number.
Access the description of detector geometry.
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
float pi
Definition: units.py:11
void GetDUNEVertDrift3ViewGlobalWire(unsigned int localWire, unsigned int plane, unsigned int tpc, unsigned int &globalWire, unsigned int &globalPlane) const
void GetCenter(double *xyz, double localz=0.0) const
Fills the world coordinate of a point on the wire.
Definition: WireGeo.cxx:73
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
geo::GeometryCore const * fGeometry
SparsePixelMap CreateSparseMap3D(detinfo::DetectorClocksData const &clockData, std::vector< art::Ptr< recob::SpacePoint > > &sp, std::vector< std::vector< art::Ptr< recob::Hit >>> &hit)
Producer algorithm for PixelMap, input to CVN neural net.
Boundary DefineBoundary(detinfo::DetectorPropertiesData const &detProp, const std::vector< const recob::Hit * > &cluster)
Get boundaries for pixel map representation of cluster.
Ionization energy from a Geant4 track.
Definition: SimChannel.h:25
void GetHitTruth(detinfo::DetectorClocksData const &clockData, art::Ptr< recob::Hit > &hit, std::vector< int > &pdgs, std::vector< int > &tracks, std::vector< float > &energies, std::vector< std::string > &processes)
Create sparse pixel map for SCN applications.
std::vector< double > fVDPlane0
QTextStream & endl(QTextStream &s)
PixelMap CreateMap(detinfo::DetectorPropertiesData const &detProp, const std::vector< art::Ptr< recob::Hit > > &slice)
WireGeo const * WirePtr(geo::WireID const &wireid) const
Returns the specified wire.