WCTrackReco_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////
2 /// \file WCTrackReco_module.cc
3 /// \module producer
4 /// \brief Reconstruction for the wire chambers used in the beamline.
5 /// Part of the beamline reconstruction for NOvA test beam.
6 /// \author Mike Wallbank (University of Cincinnati) <wallbank@fnal.gov>
7 /// Fan Gao (University of Pittsburgh) <fag16@pitt.edu>
8 /// \date December 2018
9 ////////////////////////////////////////////////////////////////////////////
10 
11 // framework
15 #include "fhiclcpp/ParameterSet.h"
23 
24 // nova
30 #include "RawData/RawBeamline.h"
31 #include "RecoBase/Vertex.h"
33 #include "Utilities/AssociationUtil.h"
34 
35 // stl
36 #include <iostream>
37 
38 // root
39 
40 // -----------------------------------------------------------------------
41 namespace beamlinereco {
42 
43  class WCTrackReco : public art::EDProducer {
44 
45  public:
46 
47  explicit WCTrackReco(const fhicl::ParameterSet& pset);
48 
49  void reconfigure(const fhicl::ParameterSet& pset);
50  void beginJob();
51  void produce(art::Event& evt);
52 
53  private:
54 
55  // Data labels
59 
60  // Configuration
61  bool fVerbose;
63  bool fUseIFDB;
64  double fVertexTime;
65 
66  // Histograms
67  TH1F* fReco_P;
68  TH1F* fY_Kink;
69  TH1F* fX_Dist;
70  TH1F* fY_Dist;
71  TH1F* fZ_Dist;
72  TH1F* fX_Mag_Int;
73  TH1F* fY_Mag_Int;
74  TH1F* fZ_Mag_Int;
76  TH1F* fX_Face_Dist;
77  TH1F* fY_Face_Dist;
78  TH1F* fTheta_Dist;
79  TH1F* fPhi_Dist;
80  TH1F* fTrack_Type;
81  TH1F* fWCDist;
82 
83  // Algorithms and services
87 
88  };
89 
91 
92 }
93 
94 // -----------------------------------------------------------------------
96  fWCHitFinderAlg(pset.get<fhicl::ParameterSet>("WCHitFinderAlg")),
97  fWCTrackAlg(pset.get<fhicl::ParameterSet>("WCTrackAlg")) {
98 
99  this->reconfigure(pset);
100  produces<std::vector<brb::WCTrack> >();
101  produces<art::Assns<brb::WCTrack, rawdata::RawBeamlineWC> >();
102  produces<std::vector<rb::Vertex> > ();
103  produces<art::Assns<rb::Vertex, brb::WCTrack> > ();
104 }
105 
106 // -----------------------------------------------------------------------
108  fRawWCDataLabel = pset.get<art::InputTag>("RawWCDataLabel");
109  fRawBeamlineConfigLabel = pset.get<std::string> ("RawBeamlineConfigLabel");
110  fMCenterDataLabel = pset.get<std::string> ("MCenterDataLabel");
111  fVerbose = pset.get<bool> ("Verbose", false);
112  fUseIFDB = pset.get<bool> ("UseIFDB");
113  fVertexTime = pset.get<double> ("VertexTime");
114 }
115 
116 // -----------------------------------------------------------------------
119 
120  fWCDist = tfs->make<TH1F>("WCCond","WC Conditions",7,0,7);
121  fReco_P = tfs->make<TH1F>("Reco_P","Reconstructed momentum", 180, 0, 1800);
122  fY_Kink = tfs->make<TH1F>("Y_Kink","Angle between US/DS tracks in Y direction (degrees)",200,-10*3.1415926/180,10*3.141592654/180);
123  fX_Dist = tfs->make<TH1F>("X_Dist","X distance between US/DS tracks at midplane (mm)",1200,-60,1260); //Do not trust these values
124  fY_Dist = tfs->make<TH1F>("Y_Dist","Y distance between US/DS tracks at midplane (mm)",1200,-600,600); //Do not trust these values
125  fZ_Dist = tfs->make<TH1F>("Z_Dist","Z distance between US/DS tracks at midplane (mm)",1200,-60,1260); //Do not trust these values
126  fX_Mag_Int = tfs->make<TH1F>("X_Mag_Int","X intercept of WC track with magnet front face (cm)",100,-150,-100);
127  fX_Mag_Int->GetXaxis()->SetTitle("X intercept [cm]");
128  fX_Mag_Int->GetYaxis()->SetTitle("Tracks per 2 cm");
129  fY_Mag_Int = tfs->make<TH1F>("Y_Mag_Int","Y intercept of WC track with magnet front face (cm)",100,-10,10);
130  fY_Mag_Int->GetXaxis()->SetTitle("Y intercept [cm]");
131  fY_Mag_Int->GetYaxis()->SetTitle("Tracks per 0.22 cm");
132  fZ_Mag_Int = tfs->make<TH1F>("Z_Mag_Int","Z intercept of WC track with magnet front face (cm)",100,415,425);
133  fZ_Mag_Int->GetXaxis()->SetTitle("Z intercept [cm]");
134  fZ_Mag_Int->GetYaxis()->SetTitle("Tracks per 0.1 cm");
135  fDist_to_Mag_Axis = tfs->make<TH1F>("Dist_to_Mag_Axis","Transverse distance to central axis of magnet [cm]",80,-20,20);
136  fDist_to_Mag_Axis->GetXaxis()->SetTitle("Transverse distance from WC track to central axis of magnet [cm]");
137  fDist_to_Mag_Axis->GetXaxis()->SetTitle("Tracks per 0.5 cm");
138  fX_Face_Dist = tfs->make<TH1F>("X_Face","X Location of Track's Detector Entry (mm)",1600,-200,1400);
139  fY_Face_Dist = tfs->make<TH1F>("Y_Face","Y Location of Track's Detector Entry (mm)",800,-400,400);
140  fTheta_Dist = tfs->make<TH1F>("Theta","Track Theta (w.r.t. Detector Z axis), (radians),",400,-.4,0.4);
141  fPhi_Dist = tfs->make<TH1F>("Phi","Track Phi (w.r.t. Detector X axis), (radians)",2000,-6.28318,6.28318);
142  fReco_P->GetXaxis()->SetTitle("Reconstructed momentum (MeV/c)");
143  fReco_P->GetYaxis()->SetTitle("Tracks per 10 MeV/c");
144  fY_Kink->GetXaxis()->SetTitle("Reconstructed y_kink (radians)");
145  fY_Kink->GetYaxis()->SetTitle("Tracks per 0.000872 radians");
146  fX_Dist->GetXaxis()->SetTitle("X distance between US and DS track ends");
147  fX_Dist->GetYaxis()->SetTitle("Tracks per 1 mm - DO NOT TRUST THESE VALUES");
148  fY_Dist->GetXaxis()->SetTitle("Y distance between US and DS track ends");
149  fY_Dist->GetYaxis()->SetTitle("Tracks per 1 mm - DO NOT TRUST THESE VALUES");
150  fZ_Dist->GetXaxis()->SetTitle("Z distance between US and DS track ends");
151  fZ_Dist->GetYaxis()->SetTitle("Tracks per 1 mm - DO NOT TRUST THESE VALUES");
152  fX_Face_Dist->GetXaxis()->SetTitle("X (mm)");
153  fX_Face_Dist->GetYaxis()->SetTitle("Tracks per 1 mm");
154  fY_Face_Dist->GetXaxis()->SetTitle("Y (mm)");
155  fY_Face_Dist->GetYaxis()->SetTitle("Tracks per 1 mm");
156  fTheta_Dist->GetXaxis()->SetTitle("Theta (radians)");
157  fTheta_Dist->GetYaxis()->SetTitle("Tracks per .002 radians");
158  fPhi_Dist->GetXaxis()->SetTitle("Phi (radians)");
159  fPhi_Dist->GetYaxis()->SetTitle("Tracks per 0.00628 radians");
160  fTrack_Type = tfs->make<TH1F>("TrackType","WCTrack conditions: 1=missHit,2=uniqueHits,3=lonelyHit,4=socialHits",4,0,4);
161  fTrack_Type->GetYaxis()->SetTitle("# Events");
162  fTrack_Type->GetXaxis()->SetTitle("Track Conditions");
163 
164 }
165 
166 // -----------------------------------------------------------------------
168 
169  // data products to be made
170  std::unique_ptr<std::vector<brb::WCTrack> > wc_reco(new std::vector<brb::WCTrack>);
171  std::unique_ptr<art::Assns<brb::WCTrack, rawdata::RawBeamlineWC> >
173  std::unique_ptr<std::vector<rb::Vertex> > wc_vtx(new std::vector<rb::Vertex>);
174  std::unique_ptr<art::Assns<rb::Vertex, brb::WCTrack> >
175  vtx_to_wc_trk(new art::Assns<rb::Vertex, brb::WCTrack>);
176 
177  // get the raw WC digits
179  std::vector<art::Ptr<rawdata::RawBeamlineWC> > wcs;
180  if (evt.getByLabel(fRawWCDataLabel, wcHandle))
181  art::fill_ptr_vector(wcs, wcHandle);
182 
184  evt.getByLabel(fMCenterDataLabel, ifdb);
185  // get the magnetic field
186  float b_field = 0.;
187  // If spill is not found, magnet current will be set to -1e10.
188  // In this case, fallback on final else to use beamline run to find magnet current in fcl.
189  if ( fUseIFDB && ifdb.isValid() && ifdb->mc7magnet!=-1e10 )
190  b_field = fMagneticField->FieldFromCurrent(ifdb->mc7magnet);
191  else if ( fUseIFDB && !ifdb.isValid() ) {
192  mf::LogInfo("TBSpillInfo") << "MW: Change in reconstruction, March 27.\n"
193  << "The magnet information is now obtained from IFDB via the TBSpillInfo module, which need to be run first.\n"
194  << "Be sure to have run this module before WCTrackReco. (Alternatively, turn off the fUseIFDB flag (temporary))."
195  << std::endl;
196  abort();
197  }
198  else {
200  unsigned int beamline_run = 0;
201  //float b_field = -1.;
202  if (evt.getByLabel(fRawBeamlineConfigLabel, blConfigHandle)) {
203  //b_field = blConfigHandle->BField();
204  beamline_run = blConfigHandle->BeamlineRun();
205  }
206  try {
207  b_field = fMagneticField->MagneticField(beamline_run);
208  }
209  catch (cet::exception& e) {
210  e.what();
211  }
212  }
213  fWCTrackAlg.setMagneticField(b_field);
214 
215  // Find hits
216  std::vector<std::vector<WCHitList> > good_hits; // Two vectors: WC#, axis
218  good_hits,
219  fVerbose);
220 
221  // Find tracks
222  std::vector<double> reco_p_list;
223  std::vector<TVector3> mag_int_list;
224  std::vector<double> dist_to_mag_axis_list;
225  std::vector<TVector2> face_list;
226  std::vector<std::vector<TVector3> > wc_hit_pos_list;
227  std::vector<TVector3> dir_list;
228  std::vector<double> y_kink_list;
229  std::vector<TVector3> dist_list;
230  std::vector<WCHitList> final_tracks;
231  float residual;
232  int WCMissed;
233  fWCTrackAlg.reconstructTracks(reco_p_list,
234  mag_int_list,
235  dist_to_mag_axis_list,
236  face_list,
237  wc_hit_pos_list,
238  dir_list,
239  final_tracks,
240  good_hits,
241  y_kink_list,
242  dist_list,
243  WCMissed,
244  fWCDist,
245  residual);
246 
247  for (unsigned int iNewTrack = 0; iNewTrack < final_tracks.size(); ++iNewTrack) {
248 
249  WCHitList final_track = final_tracks[iNewTrack];
250 
251  // std::cout << "Making track with p " << reco_p_list[iNewTrack] << std::endl
252  // << "y kink " << y_kink_list[iNewTrack] << std::endl
253  // << "dist (" << dist_list[iNewTrack].X() << ", " << dist_list[iNewTrack].Y() << ", " << dist_list[iNewTrack].Z() << ")" << std::endl
254  // << "face (" << face_list[iNewTrack].X() << ", " << face_list[iNewTrack].Y() << ")" << std::endl
255  // << "wc 0 hit " << wc_hit_pos_list[iNewTrack][0].X() << ", " << wc_hit_pos_list[iNewTrack][0].Y() << ", " << wc_hit_pos_list[iNewTrack][0].Z() << ")" << std::endl
256  // << "wc 1 hit " << wc_hit_pos_list[iNewTrack][1].X() << ", " << wc_hit_pos_list[iNewTrack][1].Y() << ", " << wc_hit_pos_list[iNewTrack][1].Z() << ")" << std::endl
257  // << "wc 2 hit " << wc_hit_pos_list[iNewTrack][2].X() << ", " << wc_hit_pos_list[iNewTrack][2].Y() << ", " << wc_hit_pos_list[iNewTrack][2].Z() << ")" << std::endl
258  // << "wc 3 hit " << wc_hit_pos_list[iNewTrack][3].X() << ", " << wc_hit_pos_list[iNewTrack][3].Y() << ", " << wc_hit_pos_list[iNewTrack][3].Z() << ")" << std::endl
259  // << "dir " << dir_list[iNewTrack].X() << ", " << dir_list[iNewTrack].Y() << ", " << dir_list[iNewTrack].Z() << ")" << std::endl
260  // << "theta " << dir_list[iNewTrack].Theta() << std::endl
261  // << "phi " << dir_list[iNewTrack].Phi() << std::endl
262  // << "residual " << residual << std::endl;
263 
264  brb::WCTrack the_track(reco_p_list[iNewTrack],
265  y_kink_list[iNewTrack],
266  dist_list[iNewTrack],
267  mag_int_list[iNewTrack],
268  dist_to_mag_axis_list[iNewTrack],
269  face_list[iNewTrack],
270  wc_hit_pos_list[iNewTrack],
271  dir_list[iNewTrack],
272  dir_list[iNewTrack].Theta(),
273  dir_list[iNewTrack].Phi(),
274  residual);
275 
276  fReco_P->Fill(reco_p_list[iNewTrack]);
277  fY_Kink->Fill(y_kink_list[iNewTrack]);
278  fX_Dist->Fill(dist_list[iNewTrack].X());
279  fY_Dist->Fill(dist_list[iNewTrack].Y());
280  fZ_Dist->Fill(dist_list[iNewTrack].Z());
281  fX_Mag_Int->Fill(mag_int_list[iNewTrack].X());
282  fY_Mag_Int->Fill(mag_int_list[iNewTrack].Y());
283  fZ_Mag_Int->Fill(mag_int_list[iNewTrack].Z());
284  fDist_to_Mag_Axis->Fill(dist_to_mag_axis_list[iNewTrack]);
285  fX_Face_Dist->Fill(face_list[iNewTrack].X());
286  fY_Face_Dist->Fill(face_list[iNewTrack].Y());
287  fTheta_Dist->Fill(dir_list[iNewTrack].Theta());
288  fPhi_Dist->Fill(dir_list[iNewTrack].Phi());
289 
290  // create vertex for detector
291  rb::Vertex vtx(face_list[iNewTrack].X(),face_list[iNewTrack].Y(),0,fVertexTime);
292 
293  // add to output data product and make associations
294  (*wc_reco).push_back(the_track);
295  util::CreateAssn(*this, evt, *(wc_reco.get()), wcs, *(wc_assn.get()));
296  (*wc_vtx).push_back(vtx);
297  util::CreateAssn(*this, evt, *(vtx_to_wc_trk.get()), *(wc_vtx.get()), *(wc_reco.get()));
298 
299  }
300 
301  evt.put(std::move(wc_reco));
302  evt.put(std::move(wc_assn));
303  evt.put(std::move(wc_vtx));
304  evt.put(std::move(vtx_to_wc_trk));
305 
306  return;
307 
308 }
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.
WCTrackReco(const fhicl::ParameterSet &pset)
void createHits(const std::vector< art::Ptr< rawdata::RawBeamlineWC > > &wcs, std::vector< std::vector< WCHitList > > &good_hits, bool verbose)
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
void setMagneticField(float field)
Definition: WCTrackAlg.cxx:517
unsigned int BeamlineRun() const
Definition: RawBeamline.cxx:30
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void produce(art::Event &evt)
This module reads IFDB DB and then stores BNB spill info.
Definition of the beamline offline geometry. Also implementation of a service to obtain this informat...
DEFINE_ART_MODULE(TestTMapFile)
art::ServiceHandle< beamlineutil::BeamlineMagneticField > fMagneticField
Float_t Y
Definition: plot.C:38
Float_t Z
Definition: plot.C:38
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
bool isValid() const
Definition: Handle.h:189
Utility to provide the magnetic field in the beamline for a given run number/run conditions. NB/ Currently this will use a lookup defined in a fhicl file, it should and will be replaced by an interface to a database when the field is stored in such a way.
void reconfigure(const fhicl::ParameterSet &pset)
T get(std::string const &key) const
Definition: ParameterSet.h:231
Encapsulation of reconstructed Wire Chamber track. Part of beamline reconstruction for NOvA test beam...
fvar< T > Phi(const fvar< T > &x)
Definition: Phi.hpp:12
Vertex location in position and time.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Algorithm implementation for wire chamber reconstruction. Part of the beamline reconstruction for NOv...
float FieldFromCurrent(float magnet_current)
Convert magnet current to magnetic field using field mapping information.
double mc7magnet
Tertiary beamline magnet current [Amps].
Definition: MCenterData.h:29
Raw data definitions for beamline data used in NOvA test beam experiment.
float MagneticField(unsigned int beamline_run)
Return the magnetic field for a given beamline data run.
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
Float_t e
Definition: plot.C:35
Float_t X
Definition: plot.C:38
void reconstructTracks(std::vector< double > &reco_p_list, std::vector< TVector3 > &mag_int_list, std::vector< double > &dist_to_mag_axis_list, std::vector< TVector2 > &face_list, std::vector< std::vector< TVector3 > > &wc_hit_pos_list, std::vector< TVector3 > &dir_list, std::vector< WCHitList > &event_final_tracks, std::vector< std::vector< WCHitList > > &good_hits, std::vector< double > &y_kink_list, std::vector< TVector3 > &dist_list, int &WCMissed, TH1F *&WCdistribution, float &residual)
Main function called for each event.
Definition: WCTrackAlg.cxx:63
const float Theta() const
Definition: HoughPeak.h:88