ProtoDUNEBeamlineReco_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // ProtoDUNEBeamlineReco does tracking and reconstruction from the beamline
4 // monitors
5 //
6 // 2018 Jake Calcutt, calcuttj@msu.edu
7 //
8 ////////////////////////////////////////////////////////////////////////
9 
10 // Framework includes
15 #include "art_root_io/TFileService.h"
18 #include "fhiclcpp/ParameterSet.h"
20 #include "ProtoDUNETruthUtils.h"
21 
23 
24 #include "TTree.h"
25 
26 namespace protoana{
27  class ProtoDUNEBeamlineReco;
28 }
29 
31 public:
32 
33  explicit ProtoDUNEBeamlineReco(fhicl::ParameterSet const& pset);
34  virtual ~ProtoDUNEBeamlineReco();
35 
36  virtual void beginJob() override;
37  virtual void endJob() override;
38 
39  void analyze(art::Event const & evt) override;
40  void reconfigure(fhicl::ParameterSet const& pset);
41  void reset();
42 
43 private:
44 
48 
49  TTree * fOutTree;
50 
51  double tof;
52  int event, run;
53  int chan;
54  double momentum;
55  int c0, c1;
56  std::vector<int> glitches_p1, glitches_p2, glitches_p3;
58  int nMomenta;
59  std::vector<double> momenta;
60  std::vector<int> possible_pdg;
61  int nTracks;
62  std::vector<double> trackEndX, trackEndY, trackEndZ;
63  std::vector<double> dirX, dirY, dirZ;
64 
65 
66  std::vector<short> fibers_h_upstream, fibers_v_upstream,
69  fibers_p3;
70 
71  unsigned long long GenTrigTS;
72  bool perfectP;
73 
74  int true_PDG;
75  double true_P;
76  double true_X, true_Y, true_Z;
77 };
78 
79 //-----------------------------------------------------------------------
81  EDAnalyzer(pset),
82  fBeamlineUtils(pset.get<fhicl::ParameterSet>("BeamlineUtils")),
83  fDoGlitches( pset.get< bool >( "DoGlitches" ) ),
84  fReferenceMomentum( pset.get< int >( "ReferenceMomentum" ) )
85 {
86  this->reconfigure(pset);
87 }
88 
89 //-----------------------------------------------------------------------
91 
92 //-----------------------------------------------------------------------
94 
96 
97  fOutTree = tfs->make<TTree>("tree", "");
98  fOutTree->Branch("event", &event);
99  fOutTree->Branch("run", &run);
100  fOutTree->Branch("glitches_p1", &glitches_p1);
101  fOutTree->Branch("glitches_p2", &glitches_p2);
102  fOutTree->Branch("glitches_p3", &glitches_p3);
103  fOutTree->Branch("TOF", &tof);
104  fOutTree->Branch("perfectP", &perfectP);
105  fOutTree->Branch("Chan", &chan);
106  fOutTree->Branch("Momentum", &momentum);
107  fOutTree->Branch("C0",&c0);
108  fOutTree->Branch("C1",&c1);
109  fOutTree->Branch("nMomenta", &nMomenta);
110  fOutTree->Branch("Momenta", &momenta);
111  fOutTree->Branch("nTracks", &nTracks);
112  fOutTree->Branch("trackEndX", &trackEndX);
113  fOutTree->Branch("trackEndY", &trackEndY);
114  fOutTree->Branch("trackEndZ", &trackEndZ);
115  fOutTree->Branch("dirX", &dirX);
116  fOutTree->Branch("dirY", &dirY);
117  fOutTree->Branch("dirZ", &dirZ);
118  fOutTree->Branch("possible_pdg", &possible_pdg);
119 
120  //Tracking Monitors
121  fOutTree->Branch("fibers_h_upstream", &fibers_h_upstream);
122  fOutTree->Branch("fibers_v_upstream", &fibers_v_upstream);
123  fOutTree->Branch("fibers_h_downstream", &fibers_h_downstream);
124  fOutTree->Branch("fibers_v_downstream", &fibers_v_downstream);
125 
126  fOutTree->Branch("glitches_h_upstream", &glitches_h_upstream);
127  fOutTree->Branch("glitches_v_upstream", &glitches_v_upstream);
128  fOutTree->Branch("glitches_h_downstream", &glitches_h_downstream);
129  fOutTree->Branch("glitches_v_downstream", &glitches_v_downstream);
130 
131  //Momentum Monitors
132  fOutTree->Branch("fibers_p1", &fibers_p1);
133  fOutTree->Branch("fibers_p2", &fibers_p2);
134  fOutTree->Branch("fibers_p3", &fibers_p3);
135 
136  //Timestamp
137  fOutTree->Branch("GenTrigTS", &GenTrigTS);
138 
139  //True info
140  fOutTree->Branch("true_PDG", &true_PDG);
141  fOutTree->Branch("true_P", &true_P);
142  fOutTree->Branch("true_X", &true_X);
143  fOutTree->Branch("true_Y", &true_Y);
144  fOutTree->Branch("true_Z", &true_Z);
145 }
146 
147 //-----------------------------------------------------------------------
149 }
150 
151 //-----------------------------------------------------------------------
153 
154  reset();
155 
156  //std::cout << "Computed TOF " << fBeamlineUtils.ComputeTOF( kProton, 2.0 ) << std::endl;
157  //std::cout << "Computed Momentum " << fBeamlineUtils.ComputeMomentum( kProton, 130. ) << std::endl;
158 
159  c0 = -1;
160  c1 = -1;
161  momentum = -1.;
162  tof = -1.;
163  chan = -1;
164  GenTrigTS = 0;
165  true_PDG = -1;
166  true_P = -999.;
167  true_X = 0.;
168  true_Y = 0.;
169  true_Z = 0.;
170 
171  event = evt.id().event();
172  run = evt.run();
173 
174  //Access the Beam Event
175  /*
176  auto beamHandle = evt.getValidHandle<std::vector<beam::ProtoDUNEBeamEvent>>("beamevent");
177 
178  std::vector<art::Ptr<beam::ProtoDUNEBeamEvent>> beamVec;
179  if( beamHandle.isValid()){
180  art::fill_ptr_vector(beamVec, beamHandle);
181  }
182 
183  const beam::ProtoDUNEBeamEvent & beamEvent = *(beamVec.at(0)); //Should just have one
184  */
185  const beam::ProtoDUNEBeamEvent & beamEvent = fBeamlineUtils.GetBeamEvent(evt);
186  /////////////////////////////////////////////////////////////
187 
188 
189  //Check the quality of the event
190  std::cout << "Timing Trigger: " << beamEvent.GetTimingTrigger() << std::endl;
191  std::cout << "Is Matched: " << beamEvent.CheckIsMatched() << std::endl << std::endl;
192 
194  std::cout << "Failed quality check" << std::endl;
195  return;
196  }
197 
198  std::cout << "Passed quality check!" << std::endl << std::endl;
199  /////////////////////////////////////////////////////////////
200 
201 
202  //Access momentum
203  const std::vector< double > & the_momenta = beamEvent.GetRecoBeamMomenta();
204  std::cout << "Number of reconstructed momenta: " << the_momenta.size() << std::endl;
205 
206  if( the_momenta.size() > 0 )
207  std::cout << "Measured Momentum: " << the_momenta[0] << std::endl;
208 
209  if( the_momenta.size() == 1)
210  momentum = the_momenta[0];
211 
212  momenta.insert( momenta.end(), the_momenta.begin(), the_momenta.end() );
213  /*
214  for( size_t i = 0; i < the_momenta.size(); ++i ){
215  momenta.push_back( the_momenta[i] );
216  }
217  */
218  nMomenta = momenta.size();
219  /////////////////////////////////////////////////////////////
220 
221 
222  std::cout << "Current: " << beamEvent.GetMagnetCurrent() << std::endl;
223 
224  //Access time of flight
225  const std::vector< double > & the_tofs = beamEvent.GetTOFs();
226  const std::vector< int > & the_chans = beamEvent.GetTOFChans();
227 
228  std::cout << "Number of measured TOF: " << the_tofs.size() << std::endl;
229  std::cout << "First TOF: " << beamEvent.GetTOF() << std::endl;
230  std::cout << "First TOF Channel: " << beamEvent.GetTOFChan() << std::endl << std::endl;
231 
232  std::cout << "All (TOF, Channels): " << std::endl;
233  for( size_t i = 0; i < the_tofs.size(); ++i ){
234  std::cout << "\t(" << the_tofs[i] << ", " << the_chans[i] << ")" << std::endl;
235  }
236  std::cout << std::endl;
237 
238  if( the_tofs.size() > 0){
239  tof = the_tofs[0];
240  chan = the_chans[0];
241  }
242  /////////////////////////////////////////////////////////////
243 
244 
245  //Access Cerenkov info
246  std::cout << "Cerenkov status, pressure:" << std::endl;
247  std::cout << "C0: " << beamEvent.GetCKov0Status() << ", " << beamEvent.GetCKov0Pressure() << std::endl;
248  std::cout << "C1: " << beamEvent.GetCKov1Status() << ", " << beamEvent.GetCKov1Pressure() << std::endl << std::endl;
249  c0 = beamEvent.GetCKov0Status();
250  c1 = beamEvent.GetCKov1Status();
251  /////////////////////////////////////////////////////////////
252 
253 
254 
255  //Access PID
256  std::vector< int > pids = fBeamlineUtils.GetPID( beamEvent, fReferenceMomentum );
257 
258  std::cout << "Possible particles" << std::endl;
259 
260  for( size_t i = 0; i < pids.size(); ++i ){
261  std::cout << pids[i] << std::endl;
262  }
263  std::cout << std::endl;
264 
265  possible_pdg.insert(possible_pdg.end(), pids.begin(), pids.end());
266 
268  std::cout << std::left << std::setw(10) << "electron " << candidates.electron << std::endl;
269  std::cout << std::left << std::setw(10) << "muon " << candidates.muon << std::endl;
270  std::cout << std::left << std::setw(10) << "pion " << candidates.pion << std::endl;
271  std::cout << std::left << std::setw(10) << "kaon " << candidates.kaon << std::endl;
272  std::cout << std::left << std::setw(10) << "proton " << candidates.proton << std::endl << std::endl;
273 
274  std::string candidates_string = fBeamlineUtils.GetPIDCandidates( beamEvent, fReferenceMomentum );
275  std::cout << candidates_string << std::endl;
276  /////////////////////////////////////////////////////////////
277 
278  //Tracking info
279  nTracks = beamEvent.GetBeamTracks().size();
280  if( nTracks == 1 ){
281  std::cout << "X " << beamEvent.GetBeamTracks()[0].Trajectory().End().X() << std::endl;
282  std::cout << "Y " << beamEvent.GetBeamTracks()[0].Trajectory().End().Y() << std::endl;
283  std::cout << "Z " << beamEvent.GetBeamTracks()[0].Trajectory().End().Z() << std::endl;
284  }
285  for (size_t i = 0; i < beamEvent.GetBeamTracks().size(); ++i) {
286  trackEndX.push_back(beamEvent.GetBeamTracks()[i].Trajectory().End().X());
287  trackEndY.push_back(beamEvent.GetBeamTracks()[i].Trajectory().End().Y());
288  trackEndZ.push_back(beamEvent.GetBeamTracks()[i].Trajectory().End().Z());
289  dirX.push_back(beamEvent.GetBeamTracks()[i].Trajectory().EndDirection().X());
290  dirY.push_back(beamEvent.GetBeamTracks()[i].Trajectory().EndDirection().Y());
291  dirZ.push_back(beamEvent.GetBeamTracks()[i].Trajectory().EndDirection().Z());
292  }
293  /////////////////////////////////////////////////////////////
294 
295  //Fibers
296  std::cout << beamEvent.GetFiberTime( "XBPF022697" ) << std::endl;
297  std::cout.precision(20);
298  unsigned long test = (unsigned long)(beamEvent.GetT0().first*1e6) + (unsigned long)(beamEvent.GetT0().second/1.e3);
299  std::cout << beamEvent.GetT0().first << " " << beamEvent.GetT0().second << std::endl;
300  std::cout << test << std::endl;
301 
302  GenTrigTS = (unsigned long)(beamEvent.GetT0().first*1e6) + (unsigned long)(beamEvent.GetT0().second/1.e3);
303 
304  fibers_p1 = beamEvent.GetActiveFibers( "XBPF022697" );
305  fibers_p2 = beamEvent.GetActiveFibers( "XBPF022701" );
306  fibers_p3 = beamEvent.GetActiveFibers( "XBPF022702" );
307 
308 
309 
310 
311 
312  fibers_h_upstream = beamEvent.GetActiveFibers( "XBPF022707" );
313  fibers_v_upstream = beamEvent.GetActiveFibers( "XBPF022708" );
314  fibers_h_downstream = beamEvent.GetActiveFibers( "XBPF022716" );
315  fibers_v_downstream = beamEvent.GetActiveFibers( "XBPF022717" );
316 
317  if( fDoGlitches ){
318  std::cout << "Doing glitches" << std::endl;
319  std::array<short,192> glitches = beamEvent.GetFBM( "XBPF022697" ).glitch_mask;
320  for( size_t i = 0; i < 192; ++i ){
321  if( glitches[i] ) glitches_p1.push_back( i );
322  }
323 
324  glitches = beamEvent.GetFBM( "XBPF022701" ).glitch_mask;
325  for( size_t i = 0; i < 192; ++i ){
326  if( glitches[i] ) glitches_p2.push_back( i );
327  }
328 
329  glitches = beamEvent.GetFBM( "XBPF022702" ).glitch_mask;
330  for( size_t i = 0; i < 192; ++i ){
331  if( glitches[i] ) glitches_p3.push_back( i );
332  }
333 
334 
335  glitches = beamEvent.GetFBM( "XBPF022707" ).glitch_mask;
336  for( size_t i = 0; i < 192; ++i ){
337  if( glitches[i] ) glitches_h_upstream.push_back( i );
338  }
339 
340  glitches = beamEvent.GetFBM( "XBPF022708" ).glitch_mask;
341  for( size_t i = 0; i < 192; ++i ){
342  if( glitches[i] ) glitches_v_upstream.push_back( i );
343  }
344 
345  glitches = beamEvent.GetFBM( "XBPF022716" ).glitch_mask;
346  for( size_t i = 0; i < 192; ++i ){
347  if( glitches[i] ) glitches_h_downstream.push_back( i );
348  }
349 
350  glitches = beamEvent.GetFBM( "XBPF022717" ).glitch_mask;
351  for( size_t i = 0; i < 192; ++i ){
352  if( glitches[i] ) glitches_v_downstream.push_back( i );
353  }
354  }
355  /////////////////////////////////////////////////////////////
356 
357  std::cout << "Perfect Momentum?" << fBeamlineUtils.HasPerfectBeamMomentum( evt ) << std::endl;
359 
360  if (!evt.isRealData()) {
361  auto mcTruths = evt.getValidHandle<std::vector<simb::MCTruth>>("generator");
363  auto true_beam_particle = truthUtil.GetGeantGoodParticle((*mcTruths)[0],evt);
364  if (!true_beam_particle) {
365  std::cout << "No true beam particle" << std::endl;
366  }
367  else {
368  std::cout << "True PDG: " << true_beam_particle->PdgCode() << std::endl;
369  true_PDG = true_beam_particle->PdgCode();
370  true_P = true_beam_particle->P();
371  size_t nPoints = true_beam_particle->NumberTrajectoryPoints();
372  for (size_t i = 0; i < nPoints; ++i) {
373  if (true_beam_particle->Position(i).Z() > 0. &&
374  true_beam_particle->Position(i).Z() < 1.) {
375  true_X = true_beam_particle->Position(i).X();
376  true_Y = true_beam_particle->Position(i).Y();
377  true_Z = true_beam_particle->Position(i).Z();
378  break;
379  }
380  }
381  }
382  }
383 
384  fOutTree->Fill();
385 }
386 
388 
390  momenta.clear();
391 
392  fibers_p1.clear();
393  fibers_p2.clear();
394  fibers_p3.clear();
395 
396  fibers_h_upstream.clear();
397  fibers_v_upstream.clear();
398  fibers_h_downstream.clear();
399  fibers_v_downstream.clear();
400  glitches_p1.clear();
401  glitches_p2.clear();
402  glitches_p3.clear();
403 
404  glitches_h_upstream.clear();
405  glitches_v_upstream.clear();
406  glitches_h_downstream.clear();
407  glitches_v_downstream.clear();
408 
409  perfectP = false;
410  possible_pdg.clear();
411 
412  trackEndX.clear();
413  trackEndY.clear();
414  trackEndZ.clear();
415 
416  dirX.clear();
417  dirY.clear();
418  dirZ.clear();
419 }
420 
const FBM & GetFBM(std::string) const
PossibleParticleCands GetPIDCandidates(beam::ProtoDUNEBeamEvent const &beamevt, double nominal_momentum)
void reconfigure(fhicl::ParameterSet const &pset)
const int & GetTimingTrigger() const
const std::vector< double > & GetTOFs() const
std::string string
Definition: nybbler.cc:12
const double & GetCKov0Pressure() const
const short & GetCKov1Status() const
bool HasPerfectBeamMomentum(art::Event const &evt) const
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
protoana::ProtoDUNEBeamlineUtils fBeamlineUtils
bool isRealData() const
int precision() const
Definition: qtextstream.h:259
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
const std::vector< short > & GetActiveFibers(std::string) const
const std::vector< double > & GetRecoBeamMomenta() const
ProtoDUNEBeamlineReco(fhicl::ParameterSet const &pset)
const double & GetTOF() const
const double & GetFiberTime(std::string) const
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:441
bool IsGoodBeamlineTrigger(art::Event const &evt) const
void analyze(art::Event const &evt) override
RunNumber_t run() const
Definition: DataViewImpl.cc:71
const std::pair< double, double > & GetT0() const
Q_EXPORT QTSManip setw(int w)
Definition: qtextstream.h:331
std::vector< int > GetPID(beam::ProtoDUNEBeamEvent const &beamevt, double nominal_momentum)
const std::vector< int > & GetTOFChans() const
Provides recob::Track data product.
std::array< short, 192 > glitch_mask
EventNumber_t event() const
Definition: EventID.h:116
const std::vector< recob::Track > & GetBeamTracks() const
const short & GetCKov0Status() const
TCEvent evt
Definition: DataStructs.cxx:7
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
const double & GetMagnetCurrent() const
const bool & CheckIsMatched() const
int bool
Definition: qglobal.h:345
const int & GetTOFChan() const
EventID id() const
Definition: Event.cc:34
const beam::ProtoDUNEBeamEvent GetBeamEvent(art::Event const &evt) const
QTextStream & endl(QTextStream &s)
Event finding and building.
const double & GetCKov1Pressure() const