ProtoDUNETruthBeamParticle_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Example how to read and save the beam truth particle.
4 // Georgios Christodoulou
5 // Georgios.Christodoulou@cern.ch
6 //
7 ////////////////////////////////////////////////////////////////////////
8 
9 // Framework includes
14 #include "art_root_io/TFileService.h"
17 #include "fhiclcpp/ParameterSet.h"
20 
21 // ROOT includes
22 #include "TTree.h"
23 #include "TFile.h"
24 #include "TString.h"
25 
26 // C++ Includes
27 #include <fstream>
28 #include <string>
29 #include <sstream>
30 #include <cmath>
31 #include <algorithm>
32 #include <iostream>
33 #include <vector>
34 
35 const int NMAXBEAMPARTICLES = 500;
36 
37 namespace protoana{
38  class ProtoDUNETruthBeamParticle;
39 }
40 
42 public:
43 
44  explicit ProtoDUNETruthBeamParticle(fhicl::ParameterSet const& pset);
46 
47  void beginJob();
48  void analyze(const art::Event& evt);
49  void reconfigure(fhicl::ParameterSet const& pset);
50  void endJob();
51 
52 private:
53 
57  int event[NMAXBEAMPARTICLES];
63  float p[NMAXBEAMPARTICLES][3];
64 
65  TTree *truth;
66 
67  // Parameters from datacard
68  std::vector<int> p_pdg;
70  bool selectall;
71 
72 };
73 
74 //-----------------------------------------------------------------------
76 
77  this->reconfigure(pset);
78 
79 }
80 
81 //-----------------------------------------------------------------------
83 
84 //-----------------------------------------------------------------------
86 
88 
89  // Define output tree
90  truth = tfs->make<TTree>("truth", "Beam Particle Truth Tree");
91  truth->Branch("nbeamparticles", &nbeamparticles, "nbeamparticles/I");
92  truth->Branch("run", run, "run[nbeamparticles]/I");
93  truth->Branch("subrun", subrun, "subrun[nbeamparticles]/I");
94  truth->Branch("event", event, "event[nbeamparticles]/I");
95  truth->Branch("pdg", pdg, "pdg[nbeamparticles]/I");
96  truth->Branch("goodparticle", goodparticle, "goodparticle[nbeamparticles]/I");
97  truth->Branch("mom", mom, "mom[nbeamparticles]/F");
98  truth->Branch("ene", ene, "ene[nbeamparticles]/F");
99  truth->Branch("startpos", startpos, "startpos[nbeamparticles][4]/F");
100  truth->Branch("p", p, "p[nbeamparticles][3]/F");
101 
102 }
103 
104 //-----------------------------------------------------------------------
106 
107  p_savegoodparticle = pset.get<bool>("SaveOnlyGoodParticle");
108  p_pdg = pset.get< std::vector<int> >("Pdg");
109  if(p_pdg[0] == 0)
110  selectall = true;
111  else
112  selectall = false;
113 
114 }
115 
116 //-----------------------------------------------------------------------
118 
119  // Access truth info
120  std::vector< art::Ptr<simb::MCTruth> > mclist;
121  auto mctruthhandle = evt.getHandle< std::vector<simb::MCTruth> >("generator");
122  if (mctruthhandle){
123  art::fill_ptr_vector(mclist,mctruthhandle);
124  }
125  else{
126  mf::LogError("protoana::ProtoDUNETruthBeamParticle::analyze") << "Requested protoDUNE beam generator information does not exist!";
127  return;
128  }
129 
130  art::Ptr<simb::MCTruth> mctruth;
131  mctruth = mclist[0];
132 
133  nbeamparticles = 0;
134  for(int i = 0; i < mctruth->NParticles(); ++i){
135  const simb::MCParticle& part(mctruth->GetParticle(i));
136 
137  // Select partciles with the chosen pdg
138  if(!selectall){
139  bool found = false;
140  for(unsigned int j = 0; j < p_pdg.size(); j++){
141  if(part.PdgCode() == p_pdg[j]){
142  found = true;
143  break;
144  }
145  }
146  if(!found) continue;
147  }
148 
149  // Option to select only good beam particles
150  bool isgoodparticle = (part.Process() == "primary");
151  if(p_savegoodparticle && !isgoodparticle) continue;
152 
153  nbeamparticles++;
154 
155  run[i] = evt.run();
156  subrun[i] = evt.subRun();
157  event[i] = evt.id().event();
158  goodparticle[i] = (part.Process() == "primary");
159  startpos[i][0] = part.Vx();
160  startpos[i][1] = part.Vy();
161  startpos[i][2] = part.Vz();
162  startpos[i][3] = part.T();
163  p[i][0] = part.Px();
164  p[i][1] = part.Py();
165  p[i][2] = part.Pz();
166  mom[i] = part.P();
167  ene[i] = part.E();
168  pdg[i] = part.PdgCode();
169  }
170 
171  if(nbeamparticles > 0)
172  truth->Fill();
173 
174 }
175 
177 
const int NMAXBEAMPARTICLES
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
ProtoDUNETruthBeamParticle(fhicl::ParameterSet const &pset)
int NParticles() const
Definition: MCTruth.h:75
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
Particle class.
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
T get(std::string const &key) const
Definition: ParameterSet.h:271
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:78
RunNumber_t run() const
Definition: DataViewImpl.cc:71
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:76
EventNumber_t event() const
Definition: EventID.h:116
TCEvent evt
Definition: DataStructs.cxx:7
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
void reconfigure(fhicl::ParameterSet const &pset)
EventID id() const
Definition: Event.cc:34
Event finding and building.