VertexFinder_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief This is the version 2, works for multi-track event
3 /// \author ricken@fnal.gov
4 /// \date Sept. 2012
5 ////////////////////////////////////////////////////////////////////////
6 
7 #include <iostream>
8 #include <vector>
9 
10 #include "TVector3.h"
11 #include "TNtuple.h"
12 
19 
20 #include "RecoBase/Track.h"
21 #include "RecoBase/Vertex.h"
22 #include "VertexFinder/DOCAInfo.h"
23 
24 
25 namespace vf {
27  {
28  public:
29  explicit VertexFinder(fhicl::ParameterSet const &pset);
30  virtual ~VertexFinder();
31 
32  void produce(art::Event &evt);
33  void beginJob();
34  void reconfigure(fhicl::ParameterSet const &pset);
35 
36  static const int GAMMAID;
37  static const int PIZEROID;
38  static const int PIPLUSID;
39  static const int PIMINUSID;
40  static const int NUEID;
41  static const int NUMUID;
42  static const int NUTAUID;
43 
44 
45  private:
46 
47  //TNtuple* ntuple;
48 
49  int fMinTrack;
50  bool fPi0Reco;
51  double fXMax;
52  double fYMax;
53  double fZMin;
54  double fZMax;
55  double fDOCA_cut;
56  double fAngle_cut;
57  };
58 }
59 
60 
61 namespace vf
62 {
63  //................................................................
64  const int VertexFinder::GAMMAID = 22;
65  const int VertexFinder::PIZEROID = 111;
66  const int VertexFinder::PIPLUSID = 211;
67  const int VertexFinder::PIMINUSID = -211;
68  const int VertexFinder::NUEID = 12;
69  const int VertexFinder::NUMUID = 14;
70  const int VertexFinder::NUTAUID = 16;
71 
72  //................................................................
74  {
75  reconfigure(pset);
76  produces< std::vector<rb::Vertex> >();
77  }
78 
79 
80  //......................................................................
82  {
83  fPi0Reco =pset.get< bool >("Pi0Reco" );
84  fMinTrack =pset.get< int >("MinTrack" );
85  fXMax =pset.get< double >("XMax" );
86  fYMax =pset.get< double >("YMax" );
87  fZMin =pset.get< double >("ZMin" );
88  fZMax =pset.get< double >("ZMax" );
89  fDOCA_cut =pset.get< double >("DOCA_cut" );
90  fAngle_cut =pset.get< double >("Angle_cut" );
91  }
92 
93  //......................................................................
95  {
96  }
97 
98  //......................................................................
100  {
101  }
102 
103  //.........................................................................
105  {
106 
107  std::unique_ptr< std::vector<rb::Vertex> > vert(new std::vector<rb::Vertex> );
108 
110  try{ evt.getByLabel("trackmerge",trackhandle); }
111  catch(...){ std::cout<<"couldn't get the trackhandle "<<std::endl; return;}
112 
113  TVector3 mc_vertex,rc_vertex, delta;
114  double doca, angle_reco;
115  std::vector<int> trackID;
116 
117  int NbTrack = trackhandle->size();
118  if (NbTrack<fMinTrack) return;
119 
120  //Get vertex time: finding the earliest cell hit in the track
121 
122 
123  double vertime;
124  art::Ptr<rb::Track> track(trackhandle,0);
125  vertime = track->Cell(0)->TNS();
126  int cell_Nb = track->NCell();
127  for (int idx =1 ; idx!=cell_Nb; idx++){
128  if (vertime > track->Cell(idx)->TNS())
129  vertime = track->Cell(idx)->TNS();
130  }
131 
132  if (NbTrack==1) rc_vertex = track->Start();
133 
134  // Let's do this with the qualified Tracks:
135  int combNB=0;
136  for (int k=0; k<NbTrack; k++) trackID.push_back(k);
137 
138  // rc_vertex.SetXYZ(0.,0.,0.);
139 
140  for (int i=0; i<NbTrack-1;i++){
141  for(int j=i+1;j<NbTrack;j++ )
142  {
143 
144  art::Ptr<rb::Track> track1(trackhandle, trackID[i]);
145  art::Ptr<rb::Track> track2(trackhandle, trackID[j]);
146  DOCAInfo doca_info(*track1, *track2);
147 
148  doca = doca_info.Doca();
149  angle_reco = doca_info.AngleReco();
150 
151  //Choosing right combinations
152  if (doca > fDOCA_cut ||
153  fabs(angle_reco) > fAngle_cut ||
154  fabs(doca_info.VertexReco().x())>fXMax ||
155  fabs(doca_info.VertexReco().y())>fYMax ||
156  doca_info.VertexReco().z()<fZMin ||
157  doca_info.VertexReco().z()>fZMax ) continue;
158 
159  else {
160  rc_vertex = rc_vertex+doca_info.VertexReco();
161  combNB++;
162  }
163 
164  }
165  }
166 
167 
168  if (combNB<1) return;
169  rc_vertex = rc_vertex*(1.0/(combNB+0.0));
170 
171  vert->push_back(rb::Vertex(rc_vertex,vertime));
172 
173  evt.put(std::move(vert));
174  return;
175  }
176 
178 
179 
180 } // end namespace vf
181 //..........................................................................
This is the version 2, works for multi-track event.
float TNS() const
Definition: CellHit.h:46
A 3D position and time representing an interaction vertex.
Definition: Vertex.h:15
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
void reconfigure(fhicl::ParameterSet const &pset)
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
VertexFinder(fhicl::ParameterSet const &pset)
double delta
Definition: runWimpSim.h:98
static const int PIPLUSID
double AngleReco() const
Definition: DOCAInfo.cxx:65
static const int NUTAUID
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
static const int PIZEROID
static const int PIMINUSID
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
static const int GAMMAID
const double j
Definition: BetheBloch.cxx:29
Vertex location in position and time.
OStream cout
Definition: OStream.cxx:6
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
double Doca() const
Definition: DOCAInfo.cxx:59
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void produce(art::Event &evt)
Int_t trackID
Definition: plot.C:84
static const int NUMUID
TVector3 VertexReco() const
Definition: DOCAInfo.cxx:53
Float_t track
Definition: plot.C:35
static const int NUEID