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"
20 #include "art_root_io/TFileService.h"
21 #include "art_root_io/TFileDirectory.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  EDProducer(pset),
97  fWCHitFinderAlg(pset.get<fhicl::ParameterSet>("WCHitFinderAlg")),
98  fWCTrackAlg(pset.get<fhicl::ParameterSet>("WCTrackAlg")) {
99 
100  this->reconfigure(pset);
101  produces<std::vector<brb::WCTrack> >();
102  produces<art::Assns<brb::WCTrack, rawdata::RawBeamlineWC> >();
103  produces<std::vector<rb::Vertex> > ();
104  produces<art::Assns<rb::Vertex, brb::WCTrack> > ();
105 }
106 
107 // -----------------------------------------------------------------------
109  fRawWCDataLabel = pset.get<art::InputTag>("RawWCDataLabel");
110  fRawBeamlineConfigLabel = pset.get<std::string> ("RawBeamlineConfigLabel");
111  fMCenterDataLabel = pset.get<std::string> ("MCenterDataLabel");
112  fVerbose = pset.get<bool> ("Verbose", false);
113  fUseIFDB = pset.get<bool> ("UseIFDB");
114  fVertexTime = pset.get<double> ("VertexTime");
115 }
116 
117 // -----------------------------------------------------------------------
120 
121  fWCDist = tfs->make<TH1F>("WCCond","WC Conditions",7,0,7);
122  fReco_P = tfs->make<TH1F>("Reco_P","Reconstructed momentum", 180, 0, 1800);
123  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);
124  fX_Dist = tfs->make<TH1F>("X_Dist","X distance between US/DS tracks at midplane (mm)",1200,-60,1260); //Do not trust these values
125  fY_Dist = tfs->make<TH1F>("Y_Dist","Y distance between US/DS tracks at midplane (mm)",1200,-600,600); //Do not trust these values
126  fZ_Dist = tfs->make<TH1F>("Z_Dist","Z distance between US/DS tracks at midplane (mm)",1200,-60,1260); //Do not trust these values
127  fX_Mag_Int = tfs->make<TH1F>("X_Mag_Int","X intercept of WC track with magnet front face (cm)",100,-150,-100);
128  fX_Mag_Int->GetXaxis()->SetTitle("X intercept [cm]");
129  fX_Mag_Int->GetYaxis()->SetTitle("Tracks per 2 cm");
130  fY_Mag_Int = tfs->make<TH1F>("Y_Mag_Int","Y intercept of WC track with magnet front face (cm)",100,-10,10);
131  fY_Mag_Int->GetXaxis()->SetTitle("Y intercept [cm]");
132  fY_Mag_Int->GetYaxis()->SetTitle("Tracks per 0.22 cm");
133  fZ_Mag_Int = tfs->make<TH1F>("Z_Mag_Int","Z intercept of WC track with magnet front face (cm)",100,415,425);
134  fZ_Mag_Int->GetXaxis()->SetTitle("Z intercept [cm]");
135  fZ_Mag_Int->GetYaxis()->SetTitle("Tracks per 0.1 cm");
136  fDist_to_Mag_Axis = tfs->make<TH1F>("Dist_to_Mag_Axis","Transverse distance to central axis of magnet [cm]",80,-20,20);
137  fDist_to_Mag_Axis->GetXaxis()->SetTitle("Transverse distance from WC track to central axis of magnet [cm]");
138  fDist_to_Mag_Axis->GetXaxis()->SetTitle("Tracks per 0.5 cm");
139  fX_Face_Dist = tfs->make<TH1F>("X_Face","X Location of Track's Detector Entry (mm)",1600,-200,1400);
140  fY_Face_Dist = tfs->make<TH1F>("Y_Face","Y Location of Track's Detector Entry (mm)",800,-400,400);
141  fTheta_Dist = tfs->make<TH1F>("Theta","Track Theta (w.r.t. Detector Z axis), (radians),",400,-.4,0.4);
142  fPhi_Dist = tfs->make<TH1F>("Phi","Track Phi (w.r.t. Detector X axis), (radians)",2000,-6.28318,6.28318);
143  fReco_P->GetXaxis()->SetTitle("Reconstructed momentum (MeV/c)");
144  fReco_P->GetYaxis()->SetTitle("Tracks per 10 MeV/c");
145  fY_Kink->GetXaxis()->SetTitle("Reconstructed y_kink (radians)");
146  fY_Kink->GetYaxis()->SetTitle("Tracks per 0.000872 radians");
147  fX_Dist->GetXaxis()->SetTitle("X distance between US and DS track ends");
148  fX_Dist->GetYaxis()->SetTitle("Tracks per 1 mm - DO NOT TRUST THESE VALUES");
149  fY_Dist->GetXaxis()->SetTitle("Y distance between US and DS track ends");
150  fY_Dist->GetYaxis()->SetTitle("Tracks per 1 mm - DO NOT TRUST THESE VALUES");
151  fZ_Dist->GetXaxis()->SetTitle("Z distance between US and DS track ends");
152  fZ_Dist->GetYaxis()->SetTitle("Tracks per 1 mm - DO NOT TRUST THESE VALUES");
153  fX_Face_Dist->GetXaxis()->SetTitle("X (mm)");
154  fX_Face_Dist->GetYaxis()->SetTitle("Tracks per 1 mm");
155  fY_Face_Dist->GetXaxis()->SetTitle("Y (mm)");
156  fY_Face_Dist->GetYaxis()->SetTitle("Tracks per 1 mm");
157  fTheta_Dist->GetXaxis()->SetTitle("Theta (radians)");
158  fTheta_Dist->GetYaxis()->SetTitle("Tracks per .002 radians");
159  fPhi_Dist->GetXaxis()->SetTitle("Phi (radians)");
160  fPhi_Dist->GetYaxis()->SetTitle("Tracks per 0.00628 radians");
161  fTrack_Type = tfs->make<TH1F>("TrackType","WCTrack conditions: 1=missHit,2=uniqueHits,3=lonelyHit,4=socialHits",4,0,4);
162  fTrack_Type->GetYaxis()->SetTitle("# Events");
163  fTrack_Type->GetXaxis()->SetTitle("Track Conditions");
164 
165 }
166 
167 // -----------------------------------------------------------------------
169 
170  // data products to be made
171  std::unique_ptr<std::vector<brb::WCTrack> > wc_reco(new std::vector<brb::WCTrack>);
172  std::unique_ptr<art::Assns<brb::WCTrack, rawdata::RawBeamlineWC> >
174  std::unique_ptr<std::vector<rb::Vertex> > wc_vtx(new std::vector<rb::Vertex>);
175  std::unique_ptr<art::Assns<rb::Vertex, brb::WCTrack> >
176  vtx_to_wc_trk(new art::Assns<rb::Vertex, brb::WCTrack>);
177 
178  // get the raw WC digits
180  std::vector<art::Ptr<rawdata::RawBeamlineWC> > wcs;
181  if (evt.getByLabel(fRawWCDataLabel, wcHandle))
182  art::fill_ptr_vector(wcs, wcHandle);
183 
185  evt.getByLabel(fMCenterDataLabel, ifdb);
186  // get the magnetic field
187  float b_field = 0.;
188  // If spill is not found, magnet current will be set to -1e10.
189  // In this case, fallback on final else to use beamline run to find magnet current in fcl.
190  if ( fUseIFDB && ifdb.isValid() && ifdb->mc7magnet!=-1e10 )
191  b_field = fMagneticField->FieldFromCurrent(ifdb->mc7magnet);
192  else if ( fUseIFDB && !ifdb.isValid() ) {
193  mf::LogInfo("TBSpillInfo") << "MW: Change in reconstruction, March 27.\n"
194  << "The magnet information is now obtained from IFDB via the TBSpillInfo module, which need to be run first.\n"
195  << "Be sure to have run this module before WCTrackReco. (Alternatively, turn off the fUseIFDB flag (temporary))."
196  << std::endl;
197  abort();
198  }
199  else {
201  unsigned int beamline_run = 0;
202  //float b_field = -1.;
203  if (evt.getByLabel(fRawBeamlineConfigLabel, blConfigHandle)) {
204  //b_field = blConfigHandle->BField();
205  beamline_run = blConfigHandle->BeamlineRun();
206  }
207  try {
208  b_field = fMagneticField->MagneticField(beamline_run);
209  }
210  catch (cet::exception& e) {
211  e.what();
212  }
213  }
214  fWCTrackAlg.setMagneticField(b_field);
215 
216  // Find hits
217  std::vector<std::vector<WCHitList> > good_hits; // Two vectors: WC#, axis
219  good_hits,
220  fVerbose);
221 
222  // Find tracks
223  std::vector<double> reco_p_list;
224  std::vector<TVector3> mag_int_list;
225  std::vector<double> dist_to_mag_axis_list;
226  std::vector<TVector2> face_list;
227  std::vector<std::vector<TVector3> > wc_hit_pos_list;
228  std::vector<TVector3> dir_list;
229  std::vector<double> y_kink_list;
230  std::vector<TVector3> dist_list;
231  std::vector<WCHitList> final_tracks;
232  float residual;
233  int WCMissed;
234  fWCTrackAlg.reconstructTracks(reco_p_list,
235  mag_int_list,
236  dist_to_mag_axis_list,
237  face_list,
238  wc_hit_pos_list,
239  dir_list,
240  final_tracks,
241  good_hits,
242  y_kink_list,
243  dist_list,
244  WCMissed,
245  fWCDist,
246  residual);
247 
248  for (unsigned int iNewTrack = 0; iNewTrack < final_tracks.size(); ++iNewTrack) {
249 
250  WCHitList final_track = final_tracks[iNewTrack];
251 
252  // std::cout << "Making track with p " << reco_p_list[iNewTrack] << std::endl
253  // << "y kink " << y_kink_list[iNewTrack] << std::endl
254  // << "dist (" << dist_list[iNewTrack].X() << ", " << dist_list[iNewTrack].Y() << ", " << dist_list[iNewTrack].Z() << ")" << std::endl
255  // << "face (" << face_list[iNewTrack].X() << ", " << face_list[iNewTrack].Y() << ")" << std::endl
256  // << "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
257  // << "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
258  // << "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
259  // << "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
260  // << "dir " << dir_list[iNewTrack].X() << ", " << dir_list[iNewTrack].Y() << ", " << dir_list[iNewTrack].Z() << ")" << std::endl
261  // << "theta " << dir_list[iNewTrack].Theta() << std::endl
262  // << "phi " << dir_list[iNewTrack].Phi() << std::endl
263  // << "residual " << residual << std::endl;
264 
265  brb::WCTrack the_track(reco_p_list[iNewTrack],
266  y_kink_list[iNewTrack],
267  dist_list[iNewTrack],
268  mag_int_list[iNewTrack],
269  dist_to_mag_axis_list[iNewTrack],
270  face_list[iNewTrack],
271  wc_hit_pos_list[iNewTrack],
272  dir_list[iNewTrack],
273  dir_list[iNewTrack].Theta(),
274  dir_list[iNewTrack].Phi(),
275  residual);
276 
277  fReco_P->Fill(reco_p_list[iNewTrack]);
278  fY_Kink->Fill(y_kink_list[iNewTrack]);
279  fX_Dist->Fill(dist_list[iNewTrack].X());
280  fY_Dist->Fill(dist_list[iNewTrack].Y());
281  fZ_Dist->Fill(dist_list[iNewTrack].Z());
282  fX_Mag_Int->Fill(mag_int_list[iNewTrack].X());
283  fY_Mag_Int->Fill(mag_int_list[iNewTrack].Y());
284  fZ_Mag_Int->Fill(mag_int_list[iNewTrack].Z());
285  fDist_to_Mag_Axis->Fill(dist_to_mag_axis_list[iNewTrack]);
286  fX_Face_Dist->Fill(face_list[iNewTrack].X());
287  fY_Face_Dist->Fill(face_list[iNewTrack].Y());
288  fTheta_Dist->Fill(dir_list[iNewTrack].Theta());
289  fPhi_Dist->Fill(dir_list[iNewTrack].Phi());
290 
291  // create vertex for detector
292  rb::Vertex vtx(face_list[iNewTrack].X(),face_list[iNewTrack].Y(),0,fVertexTime);
293 
294  // add to output data product and make associations
295  (*wc_reco).push_back(the_track);
296  util::CreateAssn(evt, *(wc_reco.get()), wcs, *(wc_assn.get()));
297  (*wc_vtx).push_back(vtx);
298  util::CreateAssn(evt, *(vtx_to_wc_trk.get()), *(wc_vtx.get()), *(wc_reco.get()));
299 
300  }
301 
302  evt.put(std::move(wc_reco));
303  evt.put(std::move(wc_assn));
304  evt.put(std::move(wc_vtx));
305  evt.put(std::move(vtx_to_wc_trk));
306 
307  return;
308 
309 }
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
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
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
bool isValid() const
Definition: Handle.h:183
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
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
int evt
Encapsulation of reconstructed Wire Chamber track. Part of beamline reconstruction for NOvA test beam...
Vertex location in position and time.
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:291
Float_t e
Definition: plot.C:35
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:730
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
enum BeamMode string