CosmicRayVertex_module.cc
Go to the documentation of this file.
1 ///
2 /// \file CosmicRayVertex_module.cc
3 /// \brief Place a vertex on a cosmic ray track.
4 /// \author messier@indiana.edu
5 /// \version $Id:$
6 ///
13 
14 #include "RecoBase/FilterList.h"
15 #include "RecoBase/Cluster.h"
16 #include "RecoBase/Vertex.h"
17 #include "RecoBase/Track.h"
18 #include "Utilities/AssociationUtil.h"
20 #include "Geometry/Geometry.h"
21 #include "GeometryObjects/Geo.h"
22 
23 /// Place a vertex on the detector wall for a cosmic-ray track
24 namespace crvtx {
25  ///
26  /// Module for placing a vertex on a cosmic ray track
27  ///
29  private:
30  // folder labels
31  std::string fSliceLabel; ///< Where to get slices?
32  std::string fTrackLabel; ///< Where to look for tracks?
33  bool fObeyPreselection; ///< Check rb::IsFiltered?
34 
35  public:
36  //..................................................................
38  this->reconfigure(p);
39  this->produces< std::vector<rb::Vertex> > ();
40  this->produces< art::Assns<rb::Vertex, rb::Cluster> >();
41  }
42 
43  //..................................................................
45 # define SET(N) { f##N = p.get<__typeof__(f##N)>(#N); }
46  SET(ObeyPreselection);
47  SET(SliceLabel);
48  SET(TrackLabel);
49 # undef SET
50  }
51 
52  //..................................................................
53  virtual ~CosmicRayVertex() { }
54 
55  //..................................................................
56  void beginJob() { }
57 
58  //..................................................................
60  unsigned int i,j;
61  //
62  // We produce a vertex with an association to the slice clutser
63  // from which it came
64  //
65  std::unique_ptr< std::vector<rb::Vertex> >
66  vtx(new std::vector<rb::Vertex>);
67  std::unique_ptr< art::Assns<rb::Vertex, rb::Cluster> >
68  vtx_to_slice(new art::Assns<rb::Vertex, rb::Cluster>);
69 
70  //
71  // Find the slice we want to analyze and any associated tracks
72  //
74  evt.getByLabel(fSliceLabel,slice);
75 
76  art::FindMany<rb::Track> trackh(slice, evt, fTrackLabel);
77 
79 
80  //Vertex location
81  double enterp[3];
82 
83  for(i = 0; i < slice->size(); ++i) {
84  if(fObeyPreselection && rb::IsFiltered(evt, slice, i)) continue;
85  //
86  // Get a pointer to "this" slice. Skip it if its the noise slice
87  //
88  const rb::Cluster& s = (*slice)[i];
89  if (s.IsNoise()) continue;
90  //
91  // Find the tracks associated to this slice
92  //
93  std::vector<const rb::Track*> tracks = trackh.at(i);
94  //
95  // Determine the locations of the vertices for each track
96  //
97  for (j=0; j<tracks.size(); ++j) {
98  //(x,y,z) track start point
99  TVector3 st = tracks[j]->Start();
100  //track direction
101  TVector3 direct = tracks[j]->Dir();
102  //(x,y,z) track start point (double)
103  double start[3] = {st.X(),st.Y(),st.Z()};
104  //Backwards along track direction (towards edge of detector)
105  double dir[3] = {-direct.X(),-direct.Y(),-direct.Z()};
106  //
107  //Check if point in detector
108  //
109  bool is_inside = true;
110  is_inside &= (start[0]>-geom->DetHalfWidth())&&(start[0]<geom->DetHalfWidth());
111  is_inside &= (start[1]>-geom->DetHalfHeight())&&(start[1]<geom->DetHalfHeight());
112  is_inside &= (start[2]>0)&&(start[2]<geom->DetLength());
113  //If yes, vertex at extrapolated entrance point of cosmic ray
114  if (is_inside)
115  geom->DEdge(start,dir,enterp);
116  //If no, check to see if the track intercepts the detector in the opposite direction
117  else {
118  double oppDir[3] = {direct.X(), direct.Y(), direct.Z()};
119  bool intersects = geom->IntersectsDetector(start, oppDir);
120  //If it doesn't intersect, vertex at track start
121  if (!intersects) {
122  std::copy(std::begin(start), std::end(start), std::begin(enterp));
123  }
124  //If it does intersect, vertex at extrapolated entrance point of cosmic ray
125  //searched for from the outside
126  else {
127  TVector3 p1;
128  p1.SetX(st.X() + (1000*direct.X()));
129  p1.SetY(st.Y() + (1000*direct.Y()));
130  p1.SetZ(st.Z() + (1000*direct.Z()));
131  geo::ClampRayToDetectorVolume(&st, &p1, &(*geom));
132  double newStart[3] = {st.X(), st.Y(), st.Z()};
133  std::copy(std::begin(newStart), std::end(newStart), std::begin(enterp));
134  }
135  }
136  //Initial time of track
137  double t0 = tracks[j]->MinTNS();
138 
139  //
140  // Add to collection of vertices to store
141  //
142  rb::Vertex v(enterp[0],enterp[1],enterp[2],t0);
143  vtx->push_back(v);
144  //
145  // Create associations between slice and vertex
146  //
147  art::Ptr<rb::Cluster> sptr(slice, i);
148  util::CreateAssn(*this, evt, *vtx, sptr, *vtx_to_slice);
149  }
150  }
151 
152  evt.put(std::move(vtx));
153  evt.put(std::move(vtx_to_slice));
154  }
155 
156  //......................................................................
157  void endJob() { }
158  };
159 }
160 ////////////////////////////////////////////////////////////////////////
161 namespace crvtx
162 {
164 }
165 ////////////////////////////////////////////////////////////////////////
A 3D position and time representing an interaction vertex.
Definition: Vertex.h:15
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
bool fObeyPreselection
Check rb::IsFiltered?
const char * p
Definition: xmltok.h:285
CosmicRayVertex(fhicl::ParameterSet const &p)
A collection of associated CellHits.
Definition: Cluster.h:47
double DetLength() const
#define SET(N)
void reconfigure(const fhicl::ParameterSet &p)
bool ClampRayToDetectorVolume(TVector3 *p0, TVector3 *p1, const GeometryBase *geom)
If either endpoint is outside the detector move it to the edge.
Definition: Geo.cxx:640
DEFINE_ART_MODULE(TestTMapFile)
std::string fTrackLabel
Where to look for tracks?
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
const XML_Char * s
Definition: expat.h:262
double DEdge(const double *x0, const double *dx, double *exitp=0) const
Encapsulate the geometry of one entire detector (near, far, ndos)
int evt
Collect Geo headers and supply basic geometry functions.
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
Definition: FilterList.h:96
const double j
Definition: BetheBloch.cxx:29
double DetHalfHeight() const
Vertex location in position and time.
Place a vertex on the detector wall for a cosmic-ray track.
void produce(art::Event &evt)
std::string fSliceLabel
Where to get slices?
double DetHalfWidth() const
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void geom(int which=0)
Definition: geom.C:163
bool IntersectsDetector(double *xyz_cm, double *dxyz) const
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
Encapsulate the geometry of one entire detector (near, far, ndos)
enum BeamMode string