protonbeamana_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: protonbeamana
3 // Plugin Type: analyzer (art v2_11_02)
4 // File: protonbeamana_module.cc
5 // Protonbeamana module: Contact Heng-Ye Liao (liao@phys.ksu.edu)
6 // from cetlib version v3_03_01.
7 ////////////////////////////////////////////////////////////////////////
8 
9 // Framework includes
10 #include "art/Framework/Core/EDAnalyzer.h"
11 #include "art/Framework/Principal/Event.h"
12 #include "art/Framework/Principal/Handle.h"
13 #include "art/Framework/Services/Registry/ServiceHandle.h"
14 #include "art_root_io/TFileService.h"
15 #include "art_root_io/TFileDirectory.h"
16 #include "art/Framework/Core/ModuleMacros.h"
17 #include "messagefacility/MessageLogger/MessageLogger.h"
18 #include "fhiclcpp/ParameterSet.h"
19 
20 #include "art/Framework/Principal/Run.h"
21 #include "art/Framework/Principal/SubRun.h"
22 
23 #include "canvas/Utilities/InputTag.h"
24 #include "canvas/Persistency/Common/FindManyP.h"
25 
26 #include "larsim/MCCheater/BackTrackerService.h"
27 #include "larsim/MCCheater/ParticleInventoryService.h"
28 
29 #include "lardataobj/RecoBase/Hit.h"
30 #include "lardataobj/RecoBase/Track.h"
31 #include "lardataobj/RecoBase/Shower.h"
32 #include "lardataobj/RecoBase/PFParticle.h"
33 #include "nusimdata/SimulationBase/MCParticle.h"
34 #include "nusimdata/SimulationBase/MCTruth.h"
35 
36 #include "lardataobj/AnalysisBase/CosmicTag.h"
37 #include "lardataobj/AnalysisBase/T0.h"
38 #include "lardataobj/AnalysisBase/Calorimetry.h"
39 #include "lardataobj/RecoBase/SpacePoint.h"
40 #include "lardataobj/RecoBase/PointCharge.h"
41 #include "lardataobj/RecoBase/Track.h"
42 #include "lardataobj/RawData/RDTimeStamp.h"
43 
44 #include "dune/DuneObj/ProtoDUNEBeamEvent.h"
45 //#include "dune/Protodune/Analysis/ProtoDUNETrackUtils.h"
46 //#include "dune/Protodune/Analysis/ProtoDUNEShowerUtils.h"
47 //#include "dune/Protodune/Analysis/ProtoDUNETruthUtils.h"
48 //#include "dune/Protodune/Analysis/ProtoDUNEPFParticleUtils.h"
49 //#include "dune/Protodune/Analysis/ProtoDUNEDataUtils.h"
50 //#include "dune/Protodune/Analysis/ProtoDUNEBeamlineUtils.h"
51 
52 //#include "dune/Protodune/singlephase/DataUtils/ProtoDUNETrackUtils.h"
53 //#include "dune/Protodune/singlephase/DataUtils/ProtoDUNEShowerUtils.h"
54 //#include "dune/Protodune/singlephase/DataUtils/ProtoDUNETruthUtils.h"
55 //#include "dune/Protodune/singlephase/DataUtils/ProtoDUNEPFParticleUtils.h"
56 //#include "dune/Protodune/singlephase/DataUtils/ProtoDUNEDataUtils.h"
57 //#include "dune/Protodune/singlephase/DataUtils/ProtoDUNEBeamlineUtils.h"
58 
59 #include "dune/Protodune/singlephase/DataUtils/ProtoDUNEDataUtils.h"
60 #include "protoduneana/Utilities/ProtoDUNETrackUtils.h"
61 #include "protoduneana/Utilities/ProtoDUNEShowerUtils.h"
62 #include "protoduneana/Utilities/ProtoDUNETruthUtils.h"
63 #include "protoduneana/Utilities/ProtoDUNEPFParticleUtils.h"
64 //#include "protoduneana/Utilities/ProtoDUNEDataUtils.h"
65 #include "protoduneana/Utilities/ProtoDUNEBeamlineUtils.h"
66 #include "protoduneana/Utilities/ProtoDUNEBeamCuts.h"
67 #include "protoduneana/Utilities/ProtoDUNEEmptyEventFinder.h"
68 
69 #include "protoduneana/Utilities/ProtoDUNECalibration.h"
70 
71 
72 #include "larevt/SpaceChargeServices/SpaceChargeService.h"
73 
74 // ROOT includes
75 #include "TTree.h"
76 #include "TTimeStamp.h"
77 #include "TVector3.h"
78 #include "TFile.h"
79 #include "TString.h"
80 #include "TH1.h"
81 #include "TLorentzVector.h"
82 #include "TProfile.h"
83 
84 //const int kMaxTracks = 1000;
85 //const int kMaxHits = 10000;
86 
87 using namespace std;
88 
89 namespace protoana {
90  class protonbeamana;
91 }
92 
93 
94 
95 
96 class protoana::protonbeamana : public art::EDAnalyzer {
97  public:
98  explicit protonbeamana(fhicl::ParameterSet const & p);
99  // The compiler-generated destructor is fine for non-base
100  // classes without bare pointers or other resource use.
101 
102  // Plugins should not be copied or assigned.
103  protonbeamana(protonbeamana const &) = delete;
104  protonbeamana(protonbeamana &&) = delete;
105  protonbeamana & operator = (protonbeamana const &) = delete;
106  protonbeamana & operator = (protonbeamana &&) = delete;
107 
108  // Required functions.
109  void analyze(art::Event const & e) override;
110 
111  // Selected optional functions.
112  void beginJob() override;
113  void endJob() override;
114 
115  void reset();
116 
117  private:
118 
119  //const art::InputTag fSpacePointModuleLabel;
120  const art::InputTag fTrackModuleLabel;
121  //const art::InputTag fTimeDecoderModuleLabel;
122 
123  const art::InputTag fBeamModuleLabel;
124  //std::string fBeamModuleLabel;
125 
126 
127  TTree *fTree;
128  // Run information
129  int run;
130  int subrun;
131  int event;
132  //int trigger;
133  double evttime;
134  int fNactivefembs[6];
135 
136  // beam information
137  bool beam_inst_valid;
138  int beam_inst_trigger;
139  int beam_inst_nTracks, beam_inst_nMomenta;
140  std::vector<double> beam_inst_TOF;
141  std::vector< int > beam_inst_PDG_candidates, beam_inst_TOF_Chan;
142 
143  std::vector<double> beamPosx;
144  std::vector<double> beamPosy;
145  std::vector<double> beamPosz;
146 
147  std::vector<double> beamDirx;
148  std::vector<double> beamDiry;
149  std::vector<double> beamDirz;
150 
151  std::vector<double> beamMomentum;
152 
153  double tof;
154  std::vector<double> tofs;
155  std::vector<int> ch_tofs;
156  short low_pressure_status;
157  short high_pressure_status;
158  double low_pressure;
159  double high_pressure;
160 
161  //Beamline utils
162  protoana::ProtoDUNEBeamlineUtils fBeamlineUtils;
163  protoana::ProtoDUNEEmptyEventFinder fEmptyEventFinder;
164  protoana::ProtoDUNEDataUtils dataUtil;
165  bool reco_reconstructable_beam_event;
166 
167  // fcl parameters for PFP particles
168  std::string fCalorimetryTag;
169  std::string fTrackerTag;
170  std::string fHitTag;
171  std::string fShowerTag;
172  std::string fPFParticleTag;
173  //std::string fGeneratorTag;
174 
175  //Beam Momentum
176  fhicl::ParameterSet fBeamPars;
177  //bool fUseCERNCalibSelection;
178  bool fVerbose;
179 
180  //geometry
181  geo::GeometryCore const * fGeometry;
182 
183  // define parameters for primary tracks
184  std::vector<double> primtrk_startx;
185  std::vector<double> primtrk_starty;
186  std::vector<double> primtrk_startz;
187 
188  std::vector<double> primtrk_endx;
189  std::vector<double> primtrk_endy;
190  std::vector<double> primtrk_endz;
191 
192  std::vector<double> primtrk_Dirx;
193  std::vector<double> primtrk_Diry;
194  std::vector<double> primtrk_Dirz;
195  std::vector<double> primtrklen;
196  std::vector<double> primtrkID;
197  std::vector<int> primtrk_trktag;
198  //int fprimaryID;
199 
200  //Pandora slice
201  int fisprimarytrack;
202  int fisprimaryshower;
203 
204  //carlo info
205  std::vector<double> primtrk_dqdx;
206  std::vector<double> primtrk_resrange;
207  //std::vector<double> primtrk_deadwireresrange; //segfault, do NOT use to access this par
208  std::vector<double> primtrk_dedx;
209  std::vector<double> primtrk_range;
210  std::vector<double> primtrk_ke;
211  std::vector<double> primtrk_hitx;
212  std::vector<double> primtrk_hity;
213  std::vector<double> primtrk_hitz;
214  std::vector<double> primtrk_pitch;
215  std::vector<int> primtrk_wid;
216  std::vector<double> primtrk_pt;
217  std::vector<size_t> primtrk_calo_hit_index;
218  std::vector<int> primtrk_ch;
219 
220  double cosine_beam_primtrk;
221 
222  std::vector<int> pdg_code;
223  std::vector<int> n_daughter;
224  std::vector<int> n_beamparticle;
225  std::vector<int> isPrimary;
226  std::vector<int> pfp_self;
227  //std::vector<int> pfp_parent;
228  std::vector<int> pfp_daughter;
229 
230  //label overlapping tracks
231  std::vector<double> Zintersection;
232  std::vector<double> Yintersection;
233  std::vector<double> Xintersection;
234  std::vector<double> timeintersection;
235 
236 
237 };
238 
239 
240 protoana::protonbeamana::protonbeamana(fhicl::ParameterSet const & p)
241  :
242  EDAnalyzer(p),
243  //fSpacePointModuleLabel(p.get< art::InputTag >("SpacePointModuleLabel")),
244  fTrackModuleLabel(p.get< art::InputTag >("TrackModuleLabel")),
245  fBeamModuleLabel(p.get< art::InputTag >("BeamModuleLabel")),
246  fBeamlineUtils(p.get<fhicl::ParameterSet>("BeamlineUtils")),
247  fEmptyEventFinder(p.get< fhicl::ParameterSet >("EmptyEventFinder")),
248 
249  //fBeamModuleLabel(p.get<std::string>("BeamModuleLabel")),
250 
251  //fTimeDecoderModuleLabel(p.get< art::InputTag >("TimeDecoderModuleLabel")),
252 
253  dataUtil(p.get<fhicl::ParameterSet>("DataUtils")),
254  fCalorimetryTag(p.get<std::string>("CalorimetryTag")),
255  fTrackerTag(p.get<std::string>("TrackerTag")),
256  fHitTag(p.get<std::string>("HitTag")),
257  fShowerTag(p.get<std::string>("ShowerTag")),
258  fPFParticleTag(p.get<std::string>("PFParticleTag")),
259  //fGeneratorTag(p.get<std::string>("GeneratorTag")),
260  fBeamPars(p.get<fhicl::ParameterSet>("BeamPars")),
261  //fUseCERNCalibSelection(p.get<bool>("UseCERNCalibSelection")),
262  fVerbose(p.get<bool>("Verbose"))
263 {
264  //if (fSaveTrackInfo == false) fSaveCaloInfo = false;
265 
266 }
267 
268 void protoana::protonbeamana::analyze(art::Event const & evt)
269 {
270  //reset containers
271  protoana::protonbeamana::reset();
272 
273 
274  //Access the Beam Event ======================================================================//
275  auto beamHandle = evt.getValidHandle<std::vector<beam::ProtoDUNEBeamEvent>>("beamevent");
276  std::vector<art::Ptr<beam::ProtoDUNEBeamEvent>> beamVec;
277  if( beamHandle.isValid()){
278  art::fill_ptr_vector(beamVec, beamHandle);
279  }
280  const beam::ProtoDUNEBeamEvent & beamEvent = *(beamVec.at(0)); //Should just have one
281  beam_inst_trigger = beamEvent.GetTimingTrigger();
282  /////////////////////////////////////////////////////////////
283 
284  //Check the quality of the event
285  std::cout << "Timing Trigger: " << beamEvent.GetTimingTrigger() << std::endl;
286  std::cout << "Is Matched: " << beamEvent.CheckIsMatched() << std::endl;
287 
288  if( !fBeamlineUtils.IsGoodBeamlineTrigger( evt ) ){
289  std::cout << "Failed beam quality check!\n" << std::endl;
290  beam_inst_valid = false;
291  return; //returns, does nothing when !=IsGoodBeamlineTrigger
292  }
293 
294  std::cout << "Passed quality check!" << std::endl << std::endl;
295 
296  int nTracks = beamEvent.GetBeamTracks().size();
297  /////////////////////////////////////////////////////////////
298 
299  if (evt.isRealData()) {
300  std::vector< int > pdg_cands = fBeamlineUtils.GetPID( beamEvent, 1. );
301  beam_inst_PDG_candidates.insert( beam_inst_PDG_candidates.end(), pdg_cands.begin(), pdg_cands.end() ); //save info of incoming particle candidates
302  }
303 
304 
305  //Access momentum
306  const std::vector< double > & momenta = beamEvent.GetRecoBeamMomenta();
307  int nMomenta = momenta.size();
308  std::cout << "Number of reconstructed momenta: " << nMomenta << std::endl;
309 
310  beam_inst_nTracks = nTracks;
311  beam_inst_nMomenta = nMomenta;
312 
313  //if( momenta.size() > 0 )
314  //std::cout << "Measured Momentum: " << momenta.at(0) << std::endl;
315  /////////////////////////////////////////////////////////////
316 
317  //Access time of flight
318  const std::vector< double > & the_tofs = beamEvent.GetTOFs();
319  const std::vector< int > & the_chans = beamEvent.GetTOFChans();
320 
321  std::cout<<"run/subrun/event:"<<evt.run()<<"/"<<evt.subRun()<<"/"<<evt.id().event()<<std::endl;
322  std::cout << "Number of measured TOF: " << the_tofs.size() << std::endl;
323  std::cout << "First TOF: " << beamEvent.GetTOF() << std::endl;
324  std::cout << "First TOF Channel: " << beamEvent.GetTOFChan() << std::endl << std::endl;
325 
326  std::cout << "All (TOF, Channels): " << std::endl;
327  for( size_t i = 0; i < the_tofs.size(); ++i ){
328  std::cout << "\t(" << the_tofs[i] << ", " << the_chans[i] << ")" << std::endl;
329  beam_inst_TOF.push_back(the_tofs[i]);
330  beam_inst_TOF_Chan.push_back(the_chans[i]);
331  }
332  std::cout << std::endl;
333  /////////////////////////////////////////////////////////////
334 
335  //Access Cerenkov info
336  std::cout << "Cerenkov status, pressure:" << std::endl;
337  std::cout << "C0: " << beamEvent.GetCKov0Status() << ", " << beamEvent.GetCKov0Pressure() << std::endl;
338  std::cout << "C1: " << beamEvent.GetCKov1Status() << ", " << beamEvent.GetCKov1Pressure() << std::endl << std::endl;
339 
340  //if (beamEvent.GetBITrigger() == 1) {
341  //ckov0status = beamEvent.GetCKov0Status(); //low_pressure_status
342  //ckov1status = beamEvent.GetCKov1Status(); //high_pressure_status
343 
344  //ckov0pressure=beamEvent.GetCKov0Pressure();
345  //ckov1pressure=beamEvent.GetCKov1Pressure();
346  //}
347 
348  if(beamEvent.GetCKov0Pressure() < beamEvent.GetCKov1Pressure() ){
349  high_pressure_status = beamEvent.GetCKov1Status();
350  low_pressure_status = beamEvent.GetCKov0Status();
351 
352  high_pressure=beamEvent.GetCKov1Pressure();
353  low_pressure=beamEvent.GetCKov0Pressure();
354  }
355  else{
356  high_pressure_status = beamEvent.GetCKov0Status();
357  low_pressure_status = beamEvent.GetCKov1Status();
358 
359  high_pressure=beamEvent.GetCKov0Pressure();
360  low_pressure=beamEvent.GetCKov1Pressure();
361  }
362 
363  /////////////////////////////////////////////////////////////
364 
365  double nom_beam_mon=fBeamPars.get<double>("Momentum");
366  std::cout<<"Selected nominal beam momentum is: "<<nom_beam_mon<<" GeV/c"<<std::endl;
367 
368  //Access PID
369  //std::vector< int > pids = fBeamlineUtils.GetPID( beamEvent, nom_beam_mon);
370  //for( size_t i = 0; i < pids.size(); ++i ){
371  //std::cout << "pid["<<i<<"]: "<< pids[i] << std::endl;
372  //}
373 
374  PossibleParticleCands candidates = fBeamlineUtils.GetPIDCandidates( beamEvent, nom_beam_mon);
375  std::cout << "electron " << candidates.electron << std::endl;
376  std::cout << "muon " << candidates.muon << std::endl;
377  std::cout << "pion " << candidates.pion << std::endl;
378  std::cout << "kaon " << candidates.kaon << std::endl;
379  std::cout << "proton " << candidates.proton << std::endl;
380 
381  std::cout << std::endl;
382  //if (candidates.proton==0) {
383  //std::cout<<"no proton!!!!!!"<<std::endl;
384  //}
385  /////////////////////////////////////////////////////////////
386 
387  //Manual Event Selection -----------------------------------------------------------------------------------------------------------------------------------------------//
388  /*bool IsProton=false;
389  //old Tof: (!fUseCERNCalibSelection&&tof>170.)
390  if (nom_beam_mon==1.) {
391  if ((beamEvent.GetTOF()>110.&&beamEvent.GetTOF()<160.)&&low_pressure_status==0) {
392  IsProton=true;
393  std::cout<<"-->> This is a proton candidate..."<<std::endl;
394  }
395  }
396  */
397 
398  //HY::For Carlo info
399  std::vector<art::Ptr<recob::Track> > tracklist;
400  auto trackListHandle = evt.getHandle< std::vector<recob::Track> >(fTrackModuleLabel);
401  if (trackListHandle) art::fill_ptr_vector(tracklist, trackListHandle);
402  else return;
403  art::FindManyP<anab::Calorimetry> fmcal(trackListHandle, evt, fCalorimetryTag);
404  art::FindManyP<recob::PFParticle> pfp_trk_assn(trackListHandle, evt, "pandoraTrack");
405  auto allHits = evt.getValidHandle<std::vector<recob::Hit> >(fHitTag);
406 
407  //HY::Cluster info
408  //art::Handle< std::vector<recob::PFParticle> > PFPListHandle;
409  //std::vector<art::Ptr<recob::PFParticle> > pfplist;
410  //if(evt.getByLabel("pandoraTrack",trackListHandle)) art::fill_ptr_vector(tracklist, trackListHandle);
411  //if(evt.getByLabel("pandora",PFPListHandle)) art::fill_ptr_vector(pfplist, PFPListHandle);
412 
413  auto trackListHandle2 = evt.getHandle< std::vector<recob::Track> >("pandoraTrack");
414  if (trackListHandle2) art::fill_ptr_vector(tracklist, trackListHandle2);
415 
416  std::vector<art::Ptr<recob::PFParticle> > pfplist;
417  auto PFPListHandle = evt.getHandle< std::vector<recob::PFParticle> > ("pandora");
418  if (PFPListHandle) art::fill_ptr_vector(pfplist, PFPListHandle);
419 
420  // to get information about the hits
421  //art::Handle< std::vector<recob::Cluster> > clusterListHandle; // to get information about the hits
422  //std::vector<art::Ptr<recob::Cluster>> clusterlist;
423  //if(evt.getByLabel("pandora", clusterListHandle))
424  //art::fill_ptr_vector(clusterlist, clusterListHandle);
425 
426  std::vector<art::Ptr<recob::Cluster>> clusterlist;
427  auto clusterListHandle = evt.getHandle< std::vector<recob::Cluster> > ("pandora");
428  if (clusterListHandle) art::fill_ptr_vector(clusterlist, clusterListHandle);
429 
430  art::FindManyP<recob::Cluster> fmcp(PFPListHandle,evt,"pandora");
431  art::FindManyP<recob::Track> pftrack(PFPListHandle,evt,"pandoraTrack");
432  std::cout<<"number of pfp_particles "<<pfplist.size()<<std::endl;
433  //std::cout<<" size of pfParticles "<<pfParticles.size()<<std::endl;
434 
435  //HY::call geometry service
436  fGeometry = &*(art::ServiceHandle<geo::Geometry>());
437 
438 
439 
440  // Implementation of required member function here.
441  run = evt.run();
442  subrun = evt.subRun();
443  event = evt.id().event();
444  art::Timestamp ts = evt.time();
445 
446  // Is this a reconstructable beam event?
447  reco_reconstructable_beam_event = !fEmptyEventFinder.IsEmptyEvent(evt);
448  std::cout<<"reco_reconstructable_beam_event:"<<reco_reconstructable_beam_event<<std::endl;
449 
450  //std::cout<<ts.timeHigh()<<" "<<ts.timeLow()<<std::endl;
451  if (ts.timeHigh() == 0){
452  //TTimeStamp tts(ts.timeLow());
453  //evttime = tts.AsDouble();
454  TTimeStamp ts2(ts.timeLow());
455  evttime = ts2.AsDouble();
456  }
457  else{
458  //TTimeStamp tts(ts.timeHigh(), ts.timeLow());
459  //evttime = tts.AsDouble();
460  TTimeStamp ts2(ts.timeHigh(), ts.timeLow());
461  evttime = ts2.AsDouble();
462  }
463 
464  // Get number of active fembs
465  if(!evt.isRealData()){
466  for(int ik=0; ik<6; ik++)
467  fNactivefembs[ik]=20;
468  }
469  else{
470  for(int ik=0; ik<6; ik++)
471  fNactivefembs[ik]=dataUtil.GetNActiveFembsForAPA(evt,ik);
472  }
473 
474 
475  if (beamVec.size()&&candidates.proton==1) { //if beam proton
476  if (beamEvent.GetTimingTrigger()==12) { //get beam timing trigger
477  if (beamEvent.CheckIsMatched()){ //if CheckIsMatched
478  //Get TOF info
479  if (beamEvent.GetTOFChan() != -1) {//if TOFChan == -1, then there was not a successful match, if it's 0, 1, 2, or 3, then there was a good match corresponding to the different pair-wise combinations of the upstream and downstream channels
480  tof = beamEvent.GetTOF();
481  for( size_t ii = 0; ii < the_tofs.size(); ++ii ){
482  tofs.push_back(the_tofs[ii]);
483  ch_tofs.push_back(the_chans[ii]);
484  }
485  }
486  //Get beam particle trajectory info
487  //auto & tracks = beaminfo[0]->GetBeamTracks();
488  auto & tracks = beamEvent.GetBeamTracks();
489  std::cout<<"ToF:"<<tof<<" [ns]"<<std::endl;
490  std::cout<<"beam trk size:"<<tracks.size()<<std::endl;
491  for (size_t i = 0; i<tracks.size(); ++i){
492  beamPosx.push_back(tracks[i].End().X());
493  beamPosy.push_back(tracks[i].End().Y());
494  beamPosz.push_back(tracks[i].End().Z());
495  beamDirx.push_back(tracks[i].StartDirection().X());
496  beamDiry.push_back(tracks[i].StartDirection().Y());
497  beamDirz.push_back(tracks[i].StartDirection().Z());
498 
499  std::cout<<"run/subrun/evt:"<<run<<"/"<<subrun<<"/"<<event<<std::endl;
500  std::cout<<"beamPosx/beamPosy/beamPosz:"<<tracks[i].End().X()<<"/"<<tracks[i].End().Y()<<"/"<<tracks[i].End().Z()<<std::endl;
501  std::cout<<"beamDirx/beamDiry/beamDirz:"<<tracks[i].StartDirection().X()<<"/"<<tracks[i].StartDirection().Y()<<"/"<<tracks[i].StartDirection().Z()<<std::endl;
502  //std::cout<<"beamDirx^2+beamDiry^2+beamDirz^2:"<<tracks[i].StartDirection().X()*tracks[i].StartDirection().X()+tracks[i].StartDirection().Y()*tracks[i].StartDirection().Y()+tracks[i].StartDirection().Z()*tracks[i].StartDirection().Z()<<std::endl;
503  }
504  //Get reconstructed beam momentum info
505  //auto & beammom = beaminfo[0]->GetRecoBeamMomenta();
506  //auto & beammom = beamEvent.GetRecoBeamMomenta();
507  //auto & beammom = beamEvent.GetRecoBeamMomenta();
508  //std::cout<<"==============================================================="<<std::endl;
509  std::cout<<"beam mom size:"<<momenta.size()<<std::endl;
510  for (size_t i = 0; i<momenta.size(); ++i){
511  beamMomentum.push_back(momenta[i]);
512  std::cout<<"beam mom["<<i<<"]:"<<momenta[i]<<" [GeV]"<<std::endl;
513  }
514  //std::cout<<"==============================================================="<<std::endl;
515 
516  //put the track parameters in the cryo here ----------------------------------------------------------------------------//
517  /*
518  // Now we want to access the output from Pandora. This comes in the form of particle flow objects (recob::PFParticle).
519  // The primary PFParticles are those we want to consider and these PFParticles then have a hierarchy of daughters that
520  // describe the whole interaction of a given primary particle
521  //
522  // / daughter track
523  // /
524  // primary track /
525  // ---------------- ---- daughter track
526  // \
527  // /\-
528  // /\\-- daughter shower
529  //
530  // The above primary PFParticle will have links to three daughter particles, two track-like and one shower-like
531  */
532 
533  std::cout<<"\n*******************************************************"<<std::endl;
534  std::cout<<"Moving on to the PFParticle section..."<<std::endl;
535 
536  // Get the PFParticle utility
537  protoana::ProtoDUNEPFParticleUtils pfpUtil;
538 
539  // Get the track utility
540  protoana::ProtoDUNETrackUtils trackUtil;
541 
542  // Get all of the PFParticles, by default from the "pandora" product
543  auto recoParticles = evt.getValidHandle<std::vector<recob::PFParticle>>(fPFParticleTag);
544 
545  // We'd like to find the beam particle. Pandora tries to do this for us, so let's use the PFParticle utility
546  // to look for it. Pandora reconstructs slices containing one (or sometimes more) primary PFParticles. These
547  // are tagged as either beam or cosmic for ProtoDUNE. This function automatically considers only those
548  // PFParticles considered as primary
549 
550  /// Use the pandora metadata to tell us if this is a beam particle or not
551  //bool isBeamParticle=dataUtil.IsBeamParticle(fPFParticleTag, evt, );
552  //bool IsBeamParticle(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel) const;
553 
554 
555  //std::vector<const recob::PFParticle*> beamParticles = pfpUtil.GetPFParticlesFromBeamSlice(evt,fPFParticleTag);
556  auto beamParticles = pfpUtil.GetPFParticlesFromBeamSlice(evt,fPFParticleTag);
557  //HY:: From the new version, prepend recob::PFParticle* with "const"
558  if(beamParticles.size() == 0){
559  std::cout << "We found no beam particles for this event... moving on" << std::endl;
560  //return;
561  }
562 
563  //unsigned int nPrimPFPartices=pfpUtil.GetNumberPrimaryPFParticle(evt,fPFParticleTag); //count all the particles
564  //std::cout<<"--> Number of all prim. PFPartices:"<<nPrimPFPartices<<std::endl;
565 
566  //int tmp_counter=0;
567  n_beamparticle.push_back(beamParticles.size());
568  std::cout<<"we have "<<beamParticles.size()<<" beam particle(s)"<<std::endl;
569  for(const recob::PFParticle* particle : beamParticles){
570  int nTrack=0;
571  int nShower=0;
572 
573  // "particle" is the pointer to our beam particle. The recob::Track or recob::Shower object
574  // of this particle might be more helpful. These return null pointers if not track-like / shower-like
575  const recob::Track* thisTrack = pfpUtil.GetPFParticleTrack(*particle,evt,fPFParticleTag,fTrackerTag);
576  const recob::Shower* thisShower = pfpUtil.GetPFParticleShower(*particle,evt,fPFParticleTag,fShowerTag);
577  if(thisTrack != 0x0) { //thisTrack
578  std::cout << "Beam particle is track-like" << std::endl;
579 
580  fisprimarytrack = 1;
581  fisprimaryshower = 0;
582 
583  nTrack++;
584  primtrk_trktag.push_back(1);
585 
586  std::vector<anab::Calorimetry> calovector = trackUtil.GetRecoTrackCalorimetry(*thisTrack, evt, fTrackerTag, fCalorimetryTag);
587  for (auto & calo : calovector){
588  if (calo.PlaneID().Plane == 2){ //only collection plane
589  //std::cout<<"ck0"<<std::endl;
590  primtrk_range.push_back(calo.Range());
591  //std::cout<<"ck1"<<std::endl;
592  primtrk_ke.push_back(calo.KineticEnergy());
593  //std::cout<<"ck2"<<std::endl;
594  //std::cout<<"calo.Range():"<<calo.Range()<<std::endl;
595  //std::cout<<"calo.KineticEnergy():"<<calo.KineticEnergy()<<std::endl;
596  //std::cout<<"calo.dQdx().size():"<<calo.dQdx().size()<<std::endl;
597  for (size_t ihit = 0; ihit < calo.dQdx().size(); ++ihit){ //loop over hits
598  //std::cout<<"ck3--0"<<std::endl;
599  //std::cout<<"calo.dQdx["<<ihit<<"]:"<<calo.dQdx()[ihit]<<std::endl;
600  //std::cout<<"calo.dEdx["<<ihit<<"]:"<<calo.dEdx()[ihit]<<std::endl;
601  //std::cout<<"calo.DeadWireResRC["<<ihit<<"]:"<<calo.DeadWireResRC()[ihit]<<std::endl;
602  primtrk_dqdx.push_back(calo.dQdx()[ihit]);
603  primtrk_resrange.push_back(calo.ResidualRange()[ihit]);
604  //primtrk_deadwireresrange.push_back(calo.DeadWireResRC()[ihit]);
605  primtrk_dedx.push_back(calo.dEdx()[ihit]);
606  primtrk_pitch.push_back(calo.TrkPitchVec()[ihit]);
607  //std::cout<<"ck3--1"<<std::endl;
608 
609  const auto &primtrk_pos=(calo.XYZ())[ihit];
610  primtrk_hitx.push_back(primtrk_pos.X());
611  primtrk_hity.push_back(primtrk_pos.Y());
612  primtrk_hitz.push_back(primtrk_pos.Z());
613  //std::cout<<"ck3--2"<<std::endl;
614 
615  primtrk_calo_hit_index.push_back(calo.TpIndices()[ihit]);
616  const recob::Hit & theHit = (*allHits)[ calo.TpIndices()[ihit] ];
617  primtrk_wid.push_back(theHit.WireID().Wire);
618  primtrk_pt.push_back(theHit.PeakTime());
619  primtrk_ch.push_back(theHit.Channel());
620  //std::cout<<"ck3--3"<<std::endl;
621  //std::cout<<theHit.WireID().Wire<<std::endl;
622 
623  //double dedx_range=17.*pow(calo.ResidualRange()[ihit],-0.42);
624  //if (dedx_range>=17.79&&dedx_range<=18.7&&thisTrack->ID()==38) {
625  //if (dedx_range>=17.79&&dedx_range<=18.7) {
626  //std::cout<<"dqdx="<<calo.dQdx()[ihit]<<"; resrange="<<calo.ResidualRange()[ihit]<<"; pitch="<<calo.TrkPitchVec()[ihit]<<std::endl;
627  //std::cout<<"(X,Y,Z)="<<"("<<calo.XYZ()[ihit].X()<<","<<calo.XYZ()[ihit].Y()<<","<<calo.XYZ()[ihit].Z()<<")"<<std::endl;
628  //std::cout<<"wid/pt/ch/triID:"<<theHit.WireID().Wire<<"/"<<theHit.PeakTime()<<"/"<<theHit.Channel()<<"/"<<thisTrack->ID()<<std::endl;
629  //std::cout<<"wid/pt/ch/tpc:"<<theHit.WireID().Wire<<"/"<<theHit.PeakTime()<<"/"<<theHit.Channel()<<"/"<<theHit.WireID().TPC<<std::endl;
630  //}
631  //std::cout<<"(X,Y,Z)="<<"("<<primtrk_pos.X()<<","<<primtrk_pos.Y()<<","<<primtrk_pos.Z()<<")"<<std::endl;
632  //std::cout<<"(X,Y,Z)="<<"("<<tmp_primtrk_hitx[ihit]<<","<<tmp_primtrk_hity[ihit]<<","<<tmp_primtrk_hitz[ihit]<<")"<<std::endl;
633  } //loop over hits
634  //std::cout<<"ck3"<<std::endl;
635  } //only collection plane
636  }
637 
638  //HY::Here comes the start/end position of primary track
639  // Find the particle vertex. We need the tracker tag here because we need to do a bit of
640  // additional work if the PFParticle is track-like to find the vertex.
641  const TVector3 vtx = pfpUtil.GetPFParticleVertex(*particle,evt,fPFParticleTag,fTrackerTag);
642  const TVector3 vtx_end(thisTrack->Trajectory().End().X(), thisTrack->Trajectory().End().Y(), thisTrack->Trajectory().End().Z());
643  primtrk_startx.push_back(vtx.X());
644  primtrk_starty.push_back(vtx.Y());
645  primtrk_startz.push_back(vtx.Z());
646  std::cout<<"vtx_X:"<<vtx.X()<<" ; vtx_Y:"<<vtx.Y()<<" ; vtx_Z:"<<vtx.Z()<<std::endl;
647 
648  primtrk_endx.push_back(vtx_end.X());
649  primtrk_endy.push_back(vtx_end.Y());
650  primtrk_endz.push_back(vtx_end.Z());
651  std::cout<<"vtx_end_X:"<<vtx_end.X()<<" ; vtx_end_Y:"<<vtx_end.Y()<<" ; vtx_end_Z:"<<vtx_end.Z()<<std::endl;
652 
653  if (vtx.Z()==vtx_end.Z()) { //warning message if Z_start=Z_end
654  std::cout<<"WARNING!! StartZ and EndZ are the same!!"<<std::endl;
655  }
656 
657  //Use Ajib's intersection calculation function
658  //[1]identify the beam track and tag other tracks
659  std::vector<float> Stw, Endw, Stt, Endt, Stwires, Endwires, Stticks, Endticks, TPCb, TPCcl;
660  Stw.clear(); Endw.clear(); Stt.clear(); Endt.clear(); Stwires.clear(); Endwires.clear(); Stticks.clear(); Endticks.clear(); TPCb.clear(); TPCcl.clear();
661  float den;
662  float numw, numt,wire_no,ticks_no;
663 
664  for(size_t p1=0;p1<pfplist.size();p1++){
665  std::vector<art::Ptr<recob::Track>> trk=pftrack.at(p1);
666  std::vector<art::Ptr<recob::Cluster>> allClusters=fmcp.at(p1);
667  //std::cout<<"fprimaryID:"<<fprimaryID<<std::endl;
668  for(size_t c1=0;c1<allClusters.size();c1++){
669  if(allClusters[c1]->Plane().Plane!=2) continue;
670  if(trk.size() && int(trk[0].key())==(int)(thisTrack->ID())){
671  Stw.push_back(allClusters[c1]->StartWire());
672  Endw.push_back(allClusters[c1]->EndWire());
673  Stt.push_back(allClusters[c1]->StartTick());
674  Endt.push_back(allClusters[c1]->EndTick());
675  TPCb.push_back(allClusters[c1]->Plane().TPC);
676  //std::cout<<"\nst/end wire:"<<allClusters[c1]->StartWire()<<"/"<<allClusters[c1]->EndWire()<<std::endl;
677  }
678  else{
679  Stwires.push_back(allClusters[c1]->StartWire());
680  Endwires.push_back(allClusters[c1]->EndWire());
681  Stticks.push_back(allClusters[c1]->StartTick());
682  Endticks.push_back(allClusters[c1]->EndTick());
683  TPCcl.push_back(allClusters[c1]->Plane().TPC);
684  }
685  }
686  }
687  //[2]find interaction points if any (assuming all tracks are straight, find the interaction points)
688  for(size_t clt=0;clt<Stw.size();clt++){
689  for(size_t cl1=0;cl1<Stwires.size();cl1++){
690  if(TPCcl[cl1]!=TPCb[clt]) continue;
691  //std::cout<<"tpc are equal "<<std::endl;
692  den=(Stw[clt]-Endw[clt])*(Stticks[cl1]-Endticks[cl1])-(Stt[clt]-Endt[clt])*(Stwires[cl1]-Endwires[cl1]);
693  if(den==0) continue;
694  numw=(Stw[clt]*Endt[clt]-Stt[clt]*Endw[clt])*(Stwires[cl1]-Endwires[cl1])-(Stw[clt]-Endw[clt])*(Stwires[cl1]*Endticks[cl1]-Stticks[cl1]*Endwires[cl1]);
695  numt=(Stw[clt]*Endt[clt]-Stt[clt]*Endw[clt])*(Stticks[cl1]-Endticks[cl1])-(Stt[clt]-Endt[clt])*(Stwires[cl1]*Endticks[cl1]-Stticks[cl1]*Endwires[cl1]);
696  wire_no=numw/den;
697  // std::cout<<"wireno and ticks not solution "<<wire_no<<" "<<ticks_no<<std::endl;
698  ticks_no=numt/den;
699  if(((Stw[clt]<wire_no && Endw[clt]>wire_no)||(Stw[clt]>wire_no && Endw[clt]<wire_no))&&((Stt[clt]<ticks_no && Endt[clt]>ticks_no)||(Stt[clt]>ticks_no && Endt[clt]<ticks_no)) && ((Stwires[cl1]<wire_no && Endwires[cl1]>wire_no)||(Stwires[cl1]>wire_no && Endwires[cl1]<wire_no)) && ((Stticks[cl1]<ticks_no && Endticks[cl1]>ticks_no)||(Stticks[cl1]>ticks_no && Endticks[cl1]<ticks_no))) {
700  std::cout<<"intersection wire and ticks are "<<std::round(wire_no)<<" "<<ticks_no<<" Stw Endw StT EndT "<<Stwires[cl1]<<" "<<Endwires[cl1]<<" "<<Stticks[cl1]<<" "<<Endticks[cl1]<<std::endl;
701  double xyzStart[3];
702  double xyzEnd[3];
703  unsigned int wireno=std::round(wire_no);
704  geo::WireID wireid(0,TPCb[clt],2,wireno);
705  fGeometry->WireEndPoints(0,TPCb[clt],2,wireno, xyzStart, xyzEnd);
706  //fGeometry->WireEndPoints(wireid, xyzStart, xyzEnd);
707  std::cout<<"Z position of intersection = "<<xyzStart[2]<<" "<<xyzEnd[2]<<" "<<wireno<<std::endl;
708  Zintersection.push_back(xyzStart[2]);
709  Yintersection.push_back(xyzStart[1]);
710  Xintersection.push_back(xyzStart[0]);
711  timeintersection.push_back(ticks_no);
712  }
713  }
714  }
715  } //thisTrack
716  if(thisShower != 0x0) { //thisShower
717  std::cout << "Beam particle is shower-like" << std::endl;
718  fisprimarytrack = 0;
719  fisprimaryshower = 1;
720 
721  nShower++;
722  primtrk_trktag.push_back(-1);
723  } //thisShower
724 
725  //HY Add
726  pdg_code.push_back(particle->PdgCode());
727  n_daughter.push_back(particle->NumDaughters());
728  isPrimary.push_back(particle->IsPrimary());
729  pfp_self.push_back(particle->Self());
730  //pfp_parent.push_back(particle->Parent());
731 
732  std::cout<<"pdg code:"<<particle->PdgCode()<<std::endl;
733  std::cout<<"IsPrimary:"<<particle->IsPrimary()<<std::endl;
734  std::cout<<"NumDaughters:"<<particle->NumDaughters()<<std::endl;
735  std::cout<<"Self:"<<particle->Self()<<std::endl;
736  std::cout<<"Parent:"<<particle->Parent()<<std::endl;
737 
738  if ((particle->NumDaughters())>0) {
739  for (int ii=0; ii<(particle->NumDaughters());++ii) {
740  std::cout<<"Daughter["<<ii<<"]:"<<particle->Daughter(ii)<<std::endl;
741  pfp_daughter.push_back(particle->Daughter(ii));
742  }
743  }
744  else {
745  pfp_daughter.push_back(-99);
746  }
747 
748  //int pdg_code=pfpUtil.PdgCode(*particle,evt,fPFParticleTag,fShowerTag);
749  //std::cout<<"IsDaughters:"<<particle->Daughters()<<std::endl;
750  //HY::Add parameters for the primary tracks here ---------------//
751 
752  // Get track direction
753  if (thisTrack) {
754  //pdg_code=thisTrack->PdgCode();
755  //std::cout<<"pdg code:"<<pdg_code<<std::endl;
756 
757  auto trackdir = thisTrack->StartDirection();
758  std::cout<<"run/subrun/event:"<<run<<"/"<<subrun<<"/"<<event<<std::endl;
759  std::cout<<"trkDirx/trkDiry/trkDirz:"<<trackdir.X()<<"/"<<trackdir.Y()<<"/"<<trackdir.Z()<<std::endl;
760  primtrk_Dirx.push_back(trackdir.X());
761  primtrk_Diry.push_back(trackdir.Y());
762  primtrk_Dirz.push_back(trackdir.Z());
763 
764  primtrklen.push_back(thisTrack->Length()); //track length
765  std::cout<<"trk length: "<<thisTrack->Length()<<" [cm]"<<std::endl;
766  primtrkID.push_back(thisTrack->ID());
767  std::cout<<"trk ID: "<<thisTrack->ID()<<""<<std::endl; //HY::Fix me::trk ID seems wrong
768  //std::cout<<"TrackId: "<<thisTrack->TrackId()<<std::endl;
769 
770  //std::cout<<"trkDirx^2+trkDiry^2+trkDirz^2:"<<trackdir.X()*trackdir.X()+trackdir.Y()*trackdir.Y()+trackdir.Z()*trackdir.Z()<<std::endl;
771  //int nn=tracks.size()-1;
772  //cosine_beam_primtrk=tracks[nn].StartDirection().X()*trackdir.X()+tracks[nn].StartDirection().Y()*trackdir.Y()+tracks[nn].StartDirection().Z()*trackdir.Z();
773  if (tracks.size()){
774  cosine_beam_primtrk=tracks[0].StartDirection().X()*trackdir.X()+tracks[0].StartDirection().Y()*trackdir.Y()+tracks[0].StartDirection().Z()*trackdir.Z();
775  }
776  // fill a histogram of trackdir.X()*beamdir.X() + .....
777  // try to get calorimetry info of this track
778 
779  }
780  // Now we can look for the interaction point of the particle if one exists, i.e where the particle
781  // scatters off an argon nucleus. Shower-like objects won't have an interaction point, so we can
782  // check this by making sure we get a sensible position
783  const TVector3 interactionVtx = pfpUtil.GetPFParticleSecondaryVertex(*particle,evt,fPFParticleTag,fTrackerTag);
784 
785 
786  // Let's get the daughter PFParticles... we can do this simply without the utility
787  for(const int daughterID : particle->Daughters()){
788  // Daughter ID is the element of the original recoParticle vector
789  const recob::PFParticle *daughterParticle = &(recoParticles->at(daughterID));
790  std::cout << "Daughter " << daughterID << " has " << daughterParticle->NumDaughters() << " daughters" << std::endl;
791  }
792 
793  // For actually studying the objects, it is easier to have the daughters in their track and shower forms.
794  // We can use the utility to get a vector of track-like and a vector of shower-like daughters
795  const std::vector<const recob::Track*> trackDaughters = pfpUtil.GetPFParticleDaughterTracks(*particle,evt,fPFParticleTag,fTrackerTag);
796  const std::vector<const recob::Shower*> showerDaughters = pfpUtil.GetPFParticleDaughterShowers(*particle,evt,fPFParticleTag,fShowerTag);
797  std::cout << "Beam particle has " << trackDaughters.size() << " track-like daughters and " << showerDaughters.size() << " shower-like daughters." << std::endl;
798 
799  std::cout<<"# total Tracks:"<<nTrack<<std::endl;
800  std::cout<<"# total Showers:"<<nShower<<std::endl;
801  //tmp_counter++;
802  }
803  //std::cout<<"tmp_counter:"<<tmp_counter<<"\n\n\n"<<std::endl;
804  //'tmp_counter' should be the same as 'nBeamP'
805  std::cout<<"*******************************************************"<<std::endl;
806 
807  //put the track parameters in the cryo here ----------------------------------------------------------------------------//
808 
809 
810  } //if CheckIsMatched
811  } //get beam timing trigger
812 
813  } //if beam proton
814 
815  fTree->Fill();
816  }
817 
818  void protoana::protonbeamana::beginJob()
819  {
820  art::ServiceHandle<art::TFileService> tfs;
821  fTree = tfs->make<TTree>("beamana","beam analysis tree");
822  fTree->Branch("run",&run,"run/I");
823  fTree->Branch("subrun",&subrun,"subrun/I");
824  fTree->Branch("event",&event,"event/I");
825  //fTree->Branch("trigger",&trigger,"trigger/I");
826  fTree->Branch("evttime",&evttime,"evttime/D");
827  fTree->Branch("reco_reconstructable_beam_event",&reco_reconstructable_beam_event);
828 
829  fTree->Branch("Nactivefembs",&fNactivefembs,"Nactivefembs[5]/I");
830 
831  fTree->Branch("beam_inst_valid", &beam_inst_valid);
832  fTree->Branch("beam_inst_trigger", &beam_inst_trigger);
833  fTree->Branch("beam_inst_nTracks", &beam_inst_nTracks);
834  fTree->Branch("beam_inst_nMomenta", &beam_inst_nMomenta);
835 
836  fTree->Branch("beam_inst_PDG_candidates", &beam_inst_PDG_candidates);
837 
838  fTree->Branch("beam_inst_TOF", &beam_inst_TOF);
839  fTree->Branch("beam_inst_TOF_Chan", &beam_inst_TOF_Chan);
840 
841  fTree->Branch("isprimarytrack", &fisprimarytrack, "isprimarytrack/I");
842  fTree->Branch("isprimaryshower", &fisprimaryshower, "isprimaryshower/I");
843 
844  fTree->Branch("beamPosx",&beamPosx);
845  fTree->Branch("beamPosy",&beamPosy);
846  fTree->Branch("beamPosz",&beamPosz);
847  fTree->Branch("beamDirx",&beamDirx);
848  fTree->Branch("beamDiry",&beamDiry);
849  fTree->Branch("beamDirz",&beamDirz);
850  fTree->Branch("beamMomentum",&beamMomentum);
851  fTree->Branch("tof", &tof, "tof/D");
852  fTree->Branch("tofs",&tofs);
853  fTree->Branch("ch_tofs",&ch_tofs);
854  fTree->Branch("low_pressure_status", &low_pressure_status, "low_pressure_status/S");
855  fTree->Branch("high_pressure_status", &high_pressure_status, "high_pressure_status/S");
856  fTree->Branch("low_pressure", &low_pressure, "low_pressure/D");
857  fTree->Branch("high_pressure", &high_pressure, "high_pressure/D");
858 
859  fTree->Branch("cosine_beam_primtrk", &cosine_beam_primtrk, "cosine_beam_primtrk/D");
860  fTree->Branch("primtrk_startx",&primtrk_startx);
861  fTree->Branch("primtrk_starty",&primtrk_starty);
862  fTree->Branch("primtrk_startz",&primtrk_startz);
863 
864  fTree->Branch("primtrk_endx",&primtrk_endx);
865  fTree->Branch("primtrk_endy",&primtrk_endy);
866  fTree->Branch("primtrk_endz",&primtrk_endz);
867 
868  fTree->Branch("primtrk_Dirx",&primtrk_Dirx);
869  fTree->Branch("primtrk_Diry",&primtrk_Diry);
870  fTree->Branch("primtrk_Dirz",&primtrk_Dirz);
871  fTree->Branch("primtrklen",&primtrklen);
872  fTree->Branch("primtrkID",&primtrkID);
873 
874  fTree->Branch("primtrk_dqdx",&primtrk_dqdx);
875  fTree->Branch("primtrk_dedx",&primtrk_dedx);
876  fTree->Branch("primtrk_resrange",&primtrk_resrange);
877  //fTree->Branch("primtrk_deadwireresrange",&primtrk_deadwireresrange);
878  fTree->Branch("primtrk_range",&primtrk_range);
879  fTree->Branch("primtrk_hitx",&primtrk_hitx);
880  fTree->Branch("primtrk_hity",&primtrk_hity);
881  fTree->Branch("primtrk_hitz",&primtrk_hitz);
882  fTree->Branch("primtrk_ke",&primtrk_ke);
883  fTree->Branch("primtrk_pitch",&primtrk_pitch);
884  fTree->Branch("primtrk_wid",&primtrk_wid);
885  fTree->Branch("primtrk_pt",&primtrk_pt);
886  fTree->Branch("primtrk_calo_hit_index",&primtrk_calo_hit_index);
887  fTree->Branch("primtrk_ch",&primtrk_ch);
888 
889  fTree->Branch("pdg_code", &pdg_code);
890  fTree->Branch("n_beamparticle", &n_beamparticle);
891  fTree->Branch("n_daughter", &n_daughter);
892  fTree->Branch("isPrimary", &isPrimary);
893  fTree->Branch("pfp_self", &pfp_self);
894  //fTree->Branch("pfp_parent", &pfp_parent);
895  fTree->Branch("pfp_daughter", &pfp_daughter);
896  fTree->Branch("primtrk_trktag", &primtrk_trktag);
897 
898  fTree->Branch("Zintersection",&Zintersection);
899  fTree->Branch("Yintersection",&Yintersection);
900  fTree->Branch("Xintersection",&Xintersection);
901  fTree->Branch("timeintersection",&timeintersection);
902 
903  }
904 
905  void protoana::protonbeamana::endJob()
906  {
907 
908  }
909 
910  void protoana::protonbeamana::reset()
911  {
912  for(int k=0; k < 5; k++)
913  fNactivefembs[k] = -999;
914 
915 
916  beam_inst_valid = true;
917  beam_inst_trigger = -999;
918  beam_inst_nTracks = -999;
919  beam_inst_nMomenta = -999;
920  reco_reconstructable_beam_event = false;
921 
922  beam_inst_PDG_candidates.clear();
923  beam_inst_TOF.clear();
924  beam_inst_TOF_Chan.clear();
925 
926  tof = -1;
927  tofs.clear();
928  ch_tofs.clear();
929 
930  low_pressure_status = -1;
931  high_pressure_status = -1;
932  low_pressure = -99;
933  high_pressure = -99;
934 
935  beamPosx.clear();
936  beamPosy.clear();
937  beamPosz.clear();
938  beamDirx.clear();
939  beamDiry.clear();
940  beamDirz.clear();
941  beamMomentum.clear();
942 
943  fisprimarytrack = 0;
944  fisprimaryshower = 0;
945 
946  primtrk_startx.clear();
947  primtrk_starty.clear();
948  primtrk_startz.clear();
949 
950  primtrk_endx.clear();
951  primtrk_endy.clear();
952  primtrk_endz.clear();
953 
954  primtrk_Dirx.clear();
955  primtrk_Diry.clear();
956  primtrk_Dirz.clear();
957  primtrklen.clear();
958  primtrkID.clear();
959  primtrk_trktag.clear();
960 
961  //primtrk_dqdx->clear();
962  //primtrk_resrange->clear();
963  //primtrk_dedx->clear();
964  primtrk_dqdx.clear();
965  primtrk_resrange.clear();
966  //primtrk_deadwireresrange.clear();
967  primtrk_dedx.clear();
968  primtrk_range.clear();
969  primtrk_hitx.clear();
970  primtrk_hity.clear();
971  primtrk_hitz.clear();
972  primtrk_ke.clear();
973  primtrk_pitch.clear();
974  primtrk_wid.clear();
975  primtrk_pt.clear();
976  primtrk_calo_hit_index.clear();
977  primtrk_ch.clear();
978 
979  pdg_code.clear();
980  n_beamparticle.clear();
981  n_daughter.clear();
982 
983  isPrimary.clear();
984  pfp_self.clear();
985  //pfp_parent.clear();
986  pfp_daughter.clear();
987 
988  cosine_beam_primtrk=-99;
989 
990  Zintersection.clear();
991  Yintersection.clear();
992  Xintersection.clear();
993  timeintersection.clear();
994 
995 
996  }
997 
998  DEFINE_ART_MODULE(protoana::protonbeamana)
rootstat.analyze
def analyze(root, level, gtrees, gbranches, doprint)
Definition: rootstat.py:69
protoana::protonbeamana::primtrk_hitx
std::vector< double > primtrk_hitx
Definition: protonbeamana_module.cc:211
protoana::protonbeamana::beamMomentum
std::vector< double > beamMomentum
Definition: protonbeamana_module.cc:151
protoana::ProtoDUNEBeamlineUtils::GetPIDCandidates
PossibleParticleCands GetPIDCandidates(beam::ProtoDUNEBeamEvent const &beamevt, double nominal_momentum)
Definition: ProtoDUNEBeamlineUtils.cxx:445
protoana::protonbeamana::protonbeamana
protonbeamana(fhicl::ParameterSet const &p)
Definition: protonbeamana_module.cc:240
protoana::protonbeamana::fBeamPars
fhicl::ParameterSet fBeamPars
Definition: protonbeamana_module.cc:176
art::Timestamp::timeLow
constexpr std::uint32_t timeLow() const
Definition: Timestamp.h:29
beam::ProtoDUNEBeamEvent::GetTimingTrigger
const int & GetTimingTrigger() const
Definition: ProtoDUNEBeamEvent.h:100
recob::PFParticle::NumDaughters
int NumDaughters() const
Returns the number of daughter particles flowing from this one.
Definition: PFParticle.h:89
beam::ProtoDUNEBeamEvent::GetTOFs
const std::vector< double > & GetTOFs() const
Definition: ProtoDUNEBeamEvent.h:69
protoana::protonbeamana::beam_inst_TOF_Chan
std::vector< int > beam_inst_TOF_Chan
Definition: protonbeamana_module.cc:141
protoana::protonbeamana::primtrk_Dirz
std::vector< double > primtrk_Dirz
Definition: protonbeamana_module.cc:194
protoana::ProtoDUNEPFParticleUtils::GetPFParticleShower
const recob::Shower * GetPFParticleShower(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel, const std::string showerLabel) const
Get the shower associated to this particle. Returns a null pointer if not found.
Definition: ProtoDUNEPFParticleUtils.cxx:496
art::DataViewImpl::getHandle
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
string
std::string string
Definition: nybbler.cc:12
recob::Hit::WireID
geo::WireID WireID() const
Definition: Hit.h:233
beam::ProtoDUNEBeamEvent::GetCKov0Pressure
const double & GetCKov0Pressure() const
Definition: ProtoDUNEBeamEvent.h:41
recob::Track::Trajectory
const recob::TrackTrajectory & Trajectory() const
Access to the stored recob::TrackTrajectory.
Definition: Track.h:98
protoana::ProtoDUNETrackUtils::GetRecoTrackCalorimetry
std::vector< anab::Calorimetry > GetRecoTrackCalorimetry(const recob::Track &track, art::Event const &evt, const std::string trackModule, const std::string caloModule) const
Get the Calorimetry(s) from a given reco track.
Definition: ProtoDUNETrackUtils.cxx:72
art::Timestamp::timeHigh
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:34
beam::ProtoDUNEBeamEvent::GetCKov1Status
const short & GetCKov1Status() const
Definition: ProtoDUNEBeamEvent.h:38
protoana::ProtoDUNEPFParticleUtils::GetPFParticleVertex
const TVector3 GetPFParticleVertex(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel, const std::string trackLabel) const
Function to find the interaction vertex of a primary PFParticle.
Definition: ProtoDUNEPFParticleUtils.cxx:352
protoana::protonbeamana::n_beamparticle
std::vector< int > n_beamparticle
Definition: protonbeamana_module.cc:224
recob::tracking::Plane::Plane
Plane(const Point_t &planePos, const Vector_t &planeDir)
Constructor from reference position on the plane and direction orthogonal to the plane.
Definition: TrackingPlane.h:61
protoana::protonbeamana::primtrk_wid
std::vector< int > primtrk_wid
Definition: protonbeamana_module.cc:215
std
STL namespace.
PointCharge.h
Information about charge reconstructed in the active volume.
protoana::protonbeamana::primtrk_trktag
std::vector< int > primtrk_trktag
Definition: protonbeamana_module.cc:197
geo::WireID::Wire
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
protoana::protonbeamana::primtrk_hitz
std::vector< double > primtrk_hitz
Definition: protonbeamana_module.cc:213
protoana::protonbeamana::primtrk_endy
std::vector< double > primtrk_endy
Definition: protonbeamana_module.cc:189
MCParticle.h
Particle class.
protoana::ProtoDUNEEmptyEventFinder::IsEmptyEvent
bool IsEmptyEvent(const art::Event &evt) const
Definition: ProtoDUNEEmptyEventFinder.cxx:22
protoana::protonbeamana::primtrk_startx
std::vector< double > primtrk_startx
Definition: protonbeamana_module.cc:184
protoana::protonbeamana::beamPosx
std::vector< double > beamPosx
Definition: protonbeamana_module.cc:143
protoana::protonbeamana::Xintersection
std::vector< double > Xintersection
Definition: protonbeamana_module.cc:233
protoana::protonbeamana::beamDiry
std::vector< double > beamDiry
Definition: protonbeamana_module.cc:148
protoana::protonbeamana::primtrk_endx
std::vector< double > primtrk_endx
Definition: protonbeamana_module.cc:188
art::DataViewImpl::isRealData
bool isRealData() const
Definition: DataViewImpl.cc:113
reset
QTextStream & reset(QTextStream &s)
Definition: qtextstream.cpp:2048
recob::Track::StartDirection
Vector_t StartDirection() const
Access to track direction at different points.
Definition: Track.h:131
protoana::protonbeamana::fBeamlineUtils
protoana::ProtoDUNEBeamlineUtils fBeamlineUtils
Definition: protonbeamana_module.cc:162
e
const double e
Definition: gUpMuFluxGen.cxx:165
recob::Track::Length
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:167
protoana::protonbeamana::primtrk_endz
std::vector< double > primtrk_endz
Definition: protonbeamana_module.cc:190
protoana::protonbeamana::primtrk_calo_hit_index
std::vector< size_t > primtrk_calo_hit_index
Definition: protonbeamana_module.cc:217
art::DataViewImpl::time
Timestamp time() const
Definition: DataViewImpl.cc:106
DEFINE_ART_MODULE
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
breakpoints::beginJob
void beginJob()
Definition: Breakpoints.cc:14
protoana::protonbeamana::primtrk_hity
std::vector< double > primtrk_hity
Definition: protonbeamana_module.cc:212
protoana::protonbeamana::beam_inst_PDG_candidates
std::vector< int > beam_inst_PDG_candidates
Definition: protonbeamana_module.cc:141
wirecell.dfp.graph.key
def key(type, name=None)
Definition: graph.py:13
protoana::protonbeamana::primtrk_dedx
std::vector< double > primtrk_dedx
Definition: protonbeamana_module.cc:208
beam::ProtoDUNEBeamEvent::GetRecoBeamMomenta
const std::vector< double > & GetRecoBeamMomenta() const
Definition: ProtoDUNEBeamEvent.h:91
protoana::protonbeamana::primtrk_resrange
std::vector< double > primtrk_resrange
Definition: protonbeamana_module.cc:206
protoana::protonbeamana::Zintersection
std::vector< double > Zintersection
Definition: protonbeamana_module.cc:231
beam::ProtoDUNEBeamEvent::GetTOF
const double & GetTOF() const
Definition: ProtoDUNEBeamEvent.h:66
fhicl::ParameterSet::get
T get(std::string const &key) const
Definition: ParameterSet.h:271
protoana::ProtoDUNEPFParticleUtils::GetPFParticleTrack
const recob::Track * GetPFParticleTrack(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel, const std::string trackLabel) const
Get the track associated to this particle. Returns a null pointer if not found.
Definition: ProtoDUNEPFParticleUtils.cxx:476
protoana::protonbeamana::primtrk_Diry
std::vector< double > primtrk_Diry
Definition: protonbeamana_module.cc:193
art::DataViewImpl::getValidHandle
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:441
test.p
p
Definition: test.py:223
art::DataViewImpl::subRun
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:78
protoana::protonbeamana::timeintersection
std::vector< double > timeintersection
Definition: protonbeamana_module.cc:234
protoana::ProtoDUNEPFParticleUtils::GetPFParticlesFromBeamSlice
const std::vector< const recob::PFParticle * > GetPFParticlesFromBeamSlice(art::Event const &evt, const std::string particleLabel) const
Return the pointers for the PFParticles in the beam slice. Returns an empty vector is no beam slice w...
Definition: ProtoDUNEPFParticleUtils.cxx:175
protoana::protonbeamana::fEmptyEventFinder
protoana::ProtoDUNEEmptyEventFinder fEmptyEventFinder
Definition: protonbeamana_module.cc:163
protoana::protonbeamana::primtrklen
std::vector< double > primtrklen
Definition: protonbeamana_module.cc:195
protoana::ProtoDUNEBeamlineUtils::IsGoodBeamlineTrigger
bool IsGoodBeamlineTrigger(art::Event const &evt) const
Definition: ProtoDUNEBeamlineUtils.cxx:649
art::DataViewImpl::run
RunNumber_t run() const
Definition: DataViewImpl.cc:71
protoana::protonbeamana::primtrk_startz
std::vector< double > primtrk_startz
Definition: protonbeamana_module.cc:186
geo::GeometryCore
Description of geometry of one entire detector.
Definition: GeometryCore.h:1468
protoana::protonbeamana::Yintersection
std::vector< double > Yintersection
Definition: protonbeamana_module.cc:232
protoana::protonbeamana::fTrackModuleLabel
const art::InputTag fTrackModuleLabel
Definition: protonbeamana_module.cc:120
protoana::protonbeamana::beamPosz
std::vector< double > beamPosz
Definition: protonbeamana_module.cc:145
keras_to_tensorflow.int
int
Definition: keras_to_tensorflow.py:69
protoana::protonbeamana::primtrkID
std::vector< double > primtrkID
Definition: protonbeamana_module.cc:196
tracks
Definition: tracks.py:1
protoana::protonbeamana::beamDirx
std::vector< double > beamDirx
Definition: protonbeamana_module.cc:147
recob::Track::ID
int ID() const
Definition: Track.h:198
recob::Hit::PeakTime
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:218
Hit.h
Declaration of signal hit object.
protoana::ProtoDUNEBeamlineUtils::GetPID
std::vector< int > GetPID(beam::ProtoDUNEBeamEvent const &beamevt, double nominal_momentum)
Definition: ProtoDUNEBeamlineUtils.cxx:459
protoana::protonbeamana::analyze
void analyze(art::Event const &e) override
Definition: protonbeamana_module.cc:268
recob::PFParticle
Hierarchical representation of particle flow.
Definition: PFParticle.h:44
beam::ProtoDUNEBeamEvent::GetTOFChans
const std::vector< int > & GetTOFChans() const
Definition: ProtoDUNEBeamEvent.h:70
protoana::protonbeamana::pfp_daughter
std::vector< int > pfp_daughter
Definition: protonbeamana_module.cc:228
protoana::protonbeamana::primtrk_starty
std::vector< double > primtrk_starty
Definition: protonbeamana_module.cc:185
End
void End(void)
Definition: gXSecComp.cxx:210
protoana::ProtoDUNEDataUtils::GetNActiveFembsForAPA
int GetNActiveFembsForAPA(art::Event const &evt, int apa) const
Get number of active fembs in an APA.
Definition: ProtoDUNEDataUtils.cxx:49
protoana::ProtoDUNEPFParticleUtils::GetPFParticleDaughterShowers
const std::vector< const recob::Shower * > GetPFParticleDaughterShowers(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel, const std::string showerLabel) const
Get the daughter showers from the PFParticle.
Definition: ProtoDUNEPFParticleUtils.cxx:716
protoana::protonbeamana::fGeometry
geo::GeometryCore const * fGeometry
Definition: protonbeamana_module.cc:181
protoana::protonbeamana::primtrk_pt
std::vector< double > primtrk_pt
Definition: protonbeamana_module.cc:216
recob::TrackTrajectory::End
Point_t const & End() const
Returns the position of the last valid point of the trajectory [cm].
Definition: TrackTrajectory.h:257
Track.h
Provides recob::Track data product.
protoana::protonbeamana::primtrk_ke
std::vector< double > primtrk_ke
Definition: protonbeamana_module.cc:210
geo::GeometryCore::WireEndPoints
void WireEndPoints(geo::WireID const &wireid, double *xyzStart, double *xyzEnd) const
Fills two arrays with the coordinates of the wire end points.
Definition: GeometryCore.cxx:1243
protoana::protonbeamana::primtrk_range
std::vector< double > primtrk_range
Definition: protonbeamana_module.cc:209
protoana::ProtoDUNEPFParticleUtils::GetPFParticleSecondaryVertex
const TVector3 GetPFParticleSecondaryVertex(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel, const std::string trackLabel) const
Function to find the secondary interaction vertex of a primary PFParticle.
Definition: ProtoDUNEPFParticleUtils.cxx:407
art::EventID::event
EventNumber_t event() const
Definition: EventID.h:116
protoana::protonbeamana::primtrk_pitch
std::vector< double > primtrk_pitch
Definition: protonbeamana_module.cc:214
beam::ProtoDUNEBeamEvent::GetBeamTracks
const std::vector< recob::Track > & GetBeamTracks() const
Definition: ProtoDUNEBeamEvent.h:192
beam::ProtoDUNEBeamEvent::GetCKov0Status
const short & GetCKov0Status() const
Definition: ProtoDUNEBeamEvent.h:37
protoana::protonbeamana::beam_inst_TOF
std::vector< double > beam_inst_TOF
Definition: protonbeamana_module.cc:140
recob::Hit
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48
tca::evt
TCEvent evt
Definition: DataStructs.cxx:7
protoana::protonbeamana::fBeamModuleLabel
const art::InputTag fBeamModuleLabel
Definition: protonbeamana_module.cc:123
art::get
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
art::fill_ptr_vector
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
protoana::protonbeamana::tofs
std::vector< double > tofs
Definition: protonbeamana_module.cc:154
trkf::Plane
recob::tracking::Plane Plane
Definition: TrackState.h:17
beam::ProtoDUNEBeamEvent::CheckIsMatched
const bool & CheckIsMatched() const
Definition: ProtoDUNEBeamEvent.h:86
protoana::protonbeamana::beamPosy
std::vector< double > beamPosy
Definition: protonbeamana_module.cc:144
bool
int bool
Definition: qglobal.h:345
recob::Hit::Channel
raw::ChannelID_t Channel() const
ID of the readout channel the hit was extracted from.
Definition: Hit.h:230
beam::ProtoDUNEBeamEvent::GetTOFChan
const int & GetTOFChan() const
Definition: ProtoDUNEBeamEvent.h:67
protoana::ProtoDUNEPFParticleUtils::GetPFParticleDaughterTracks
const std::vector< const recob::Track * > GetPFParticleDaughterTracks(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel, const std::string trackLabel) const
Get the daughter tracks from the PFParticle.
Definition: ProtoDUNEPFParticleUtils.cxx:695
protoana::protonbeamana::primtrk_dqdx
std::vector< double > primtrk_dqdx
Definition: protonbeamana_module.cc:205
art::Event::id
EventID id() const
Definition: Event.cc:34
recob::Track
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
Definition: Track.h:49
calo
calorimetry
Definition: CaloChecker_module.cc:43
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
protoana::protonbeamana::dataUtil
protoana::ProtoDUNEDataUtils dataUtil
Definition: protonbeamana_module.cc:164
event
Event finding and building.
Definition: EventCheater_module.cc:32
protoana::protonbeamana::primtrk_Dirx
std::vector< double > primtrk_Dirx
Definition: protonbeamana_module.cc:192
protoana::protonbeamana::beamDirz
std::vector< double > beamDirz
Definition: protonbeamana_module.cc:149
beam::ProtoDUNEBeamEvent::GetCKov1Pressure
const double & GetCKov1Pressure() const
Definition: ProtoDUNEBeamEvent.h:42
Generated by   doxygen 1.8.11