All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ICARUSParticleAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: ICARUSParticleAna
3 // Plugin Type: analyzer (art v3_02_06)
4 // File: ICARUSParticleAna_module.cc
5 //
6 // Generated at Fri Oct 11 12:08:44 2019 by Laura Domine using cetskelgen
7 // from cetlib version v3_07_02.
8 ////////////////////////////////////////////////////////////////////////
9 
10 #include "art/Framework/Core/EDAnalyzer.h"
11 #include "art/Framework/Core/ModuleMacros.h"
12 #include "art/Framework/Principal/Event.h"
13 #include "art/Framework/Principal/Handle.h"
14 #include "art/Framework/Principal/Run.h"
15 #include "art/Framework/Principal/SubRun.h"
16 #include "canvas/Utilities/InputTag.h"
17 #include "fhiclcpp/ParameterSet.h"
18 #include "messagefacility/MessageLogger/MessageLogger.h"
19 
20 #include "nusimdata/SimulationBase/MCTruth.h"
21 #include "nusimdata/SimulationBase/MCParticle.h"
22 //#include "lardataobj/Simulation/SimEnergyDeposit.h"
23 #include <TTree.h>
24 #include <TFile.h>
25 #include <TLorentzVector.h>
26 
27 class ICARUSParticleAna;
28 
29 
30 class ICARUSParticleAna : public art::EDAnalyzer {
31 public:
32  explicit ICARUSParticleAna(fhicl::ParameterSet const& p);
33  // The compiler-generated destructor is fine for non-base
34  // classes without bare pointers or other resource use.
35 
36  // Plugins should not be copied or assigned.
37  ICARUSParticleAna(ICARUSParticleAna const&) = delete;
38  ICARUSParticleAna(ICARUSParticleAna&&) = delete;
39  ICARUSParticleAna& operator=(ICARUSParticleAna const&) = delete;
40  ICARUSParticleAna& operator=(ICARUSParticleAna&&) = delete;
41 
42  // Required functions.
43  void analyze(art::Event const& e) override;
44  void beginJob() override;
45  void endJob() override;
46 
47 private:
48 
49  // Declare member data here.
50  TFile *_f;
51  std::string _output_fname;
52  std::string _particle_label;
53  std::string _trajectory_label;
54 
55  TTree* _particletree;
56  int _run, _event;
57  int _pdg_code;
58  double _x, _y, _z;
59  double _time;
60  double _energy;
61  double _momentum;
62  int _track_id;
63  std::vector<double> _x_v, _y_v, _z_v, _time_v, _energy_v;
64 };
65 
66 
67 ICARUSParticleAna::ICARUSParticleAna(fhicl::ParameterSet const& p)
68  : EDAnalyzer{p} // ,
69  // More initializers here.
70 {
71  // Call appropriate consumes<>() for any products to be retrieved by this module.
72  _output_fname = p.get<std::string>("OutputFileName");
73  _particle_label = p.get<std::string>("ParticleProducer");
74  _trajectory_label = p.get<std::string>("TrajectoryProducer");
75 }
76 
77 void ICARUSParticleAna::beginJob()
78 {
79  _f = TFile::Open(_output_fname.c_str(), "RECREATE");
80  std::string name = _particle_label + "_particletree";
81  _particletree = new TTree(name.c_str(), name.c_str());
82  _particletree->Branch("run", &_run, "run/I");
83  _particletree->Branch("event", &_event, "event/I");
84  _particletree->Branch("pdg_code", &_pdg_code, "pdg_code/I");
85  _particletree->Branch("x", &_x, "x/D");
86  _particletree->Branch("y", &_y, "y/D");
87  _particletree->Branch("z", &_z, "z/D");
88  _particletree->Branch("time", &_time, "time/D");
89  _particletree->Branch("energy", &_energy, "energy/D");
90  _particletree->Branch("momentum", &_momentum, "momentum/D");
91  _particletree->Branch("track_id", &_track_id, "track_id/I");
92  _particletree->Branch("x_v", &_x_v);
93  _particletree->Branch("y_v", &_y_v);
94  _particletree->Branch("z_v", &_z_v);
95  _particletree->Branch("time_v", &_time_v);
96  _particletree->Branch("energy_v", &_energy_v);
97 }
98 
99 void ICARUSParticleAna::endJob()
100 {
101  _f->cd();
102  _particletree->Write();
103  if (_f) _f->Close();
104 }
105 
106 void ICARUSParticleAna::analyze(art::Event const& e)
107 {
108  // Implementation of required member function here.
109  _event = e.id().event();
110  _run = e.id().run();
111 
112  art::Handle<std::vector< simb::MCTruth > > mctruth_h;
113  e.getByLabel(_particle_label, mctruth_h);
114  /*for (size_t idx = 0; idx < mctruth_h->size(); ++idx) {
115  auto const& mctruth = (*mctruth_h)[idx];
116  for (int part_idx = 0; part_idx < mctruth.NParticles(); ++part_idx) {
117  const auto& particle = mctruth.GetParticle(part_idx);
118  _pdg_code = particle.PdgCode();
119  _x = particle.Vx();
120  _y = particle.Vy();
121  _z = particle.Vz();
122  _time = particle.T();
123  _energy = particle.E();
124  _momentum = particle.P();
125  _track_id = particle.TrackId();
126 
127  const auto& trajectory = particle.Trajectory();
128  _x_v.resize(trajectory.size(), 0.);
129  _y_v.resize(trajectory.size(), 0.);
130  _z_v.resize(trajectory.size(), 0.);
131  for (size_t pt_idx = 0; pt_idx < trajectory.size(); ++pt_idx) {
132  _x_v[pt_idx] = trajectory.X(pt_idx);
133  _y_v[pt_idx] = trajectory.Y(pt_idx);
134  _z_v[pt_idx] = trajectory.Z(pt_idx);
135  _time_v[pt_idx] = trajectory.T(pt_idx);
136  _energy_v[pt_idx] = trajectory.E(pt_idx);
137  }
138  for (auto const & pair : trajectory) {
139  TLorentzVector const& pos = pair.first;
140  _x_v.push_back(pos.X());
141  _y_v.push_back(pos.Y());
142  _z_v.push_back(pos.Z());
143  _time_v.push_back(pos.T());
144  _energy_v.push_back(pos.E());
145  std::cout << pos << std::endl;
146  }
147  _particletree->Fill();
148  }
149  }*/
150 
151  art::Handle<std::vector< simb::MCParticle > > mcparticle_h;
152  e.getByLabel(_trajectory_label, mcparticle_h);
153  for (size_t idx = 0; idx < mcparticle_h->size(); ++idx) {
154  auto const& particle = (*mcparticle_h)[idx];
155  //for (auto const& particle : mcparticle_v) {
156  _pdg_code = particle.PdgCode();
157  _x = particle.Vx();
158  _y = particle.Vy();
159  _z = particle.Vz();
160  _time = particle.T();
161  _energy = particle.E();
162  _momentum = particle.P();
163  _track_id = particle.TrackId();
164 
165  const auto& trajectory = particle.Trajectory();
166  _x_v.resize(trajectory.size(), 0.);
167  _y_v.resize(trajectory.size(), 0.);
168  _z_v.resize(trajectory.size(), 0.);
169  _time_v.resize(trajectory.size(), 0.);
170  _energy_v.resize(trajectory.size(), 0.);
171  for (size_t pt_idx = 0; pt_idx < trajectory.size(); ++pt_idx) {
172  _x_v[pt_idx] = trajectory.X(pt_idx);
173  _y_v[pt_idx] = trajectory.Y(pt_idx);
174  _z_v[pt_idx] = trajectory.Z(pt_idx);
175  _time_v[pt_idx] = trajectory.T(pt_idx);
176  _energy_v[pt_idx] = trajectory.E(pt_idx);
177  }
178  _particletree->Fill();
179  //}
180  }
181 
182  // Also record sim::SimEnergyDeposit from largeant?
183  //art::Handle<std::vector< sim::SimEnergyDeposit > > simenergy_h;
184  //e.getByLabel(_simenergy_label, simenergy_h);
185  // Plus Kazu's algorithm to sparsify?
186 }
187 
188 DEFINE_ART_MODULE(ICARUSParticleAna)
ICARUSParticleAna::_time_v
std::vector< double > _time_v
Definition: ICARUSParticleAna_module.cc:63
p
pdgs p
Definition: selectors.fcl:22
ICARUSParticleAna::_z_v
std::vector< double > _z_v
Definition: ICARUSParticleAna_module.cc:63
ICARUSParticleAna::operator=
ICARUSParticleAna & operator=(ICARUSParticleAna const &)=delete
ICARUSParticleAna::_energy_v
std::vector< double > _energy_v
Definition: ICARUSParticleAna_module.cc:63
ICARUSParticleAna::_x_v
std::vector< double > _x_v
Definition: ICARUSParticleAna_module.cc:63
e
do i e
Definition: updateVersionList.sh:68
name
then echo fcl name
Definition: sbndpoms_genfclwithrunnumber_maker.sh:220
ICARUSParticleAna::_y_v
std::vector< double > _y_v
Definition: ICARUSParticleAna_module.cc:63
ICARUSParticleAna::analyze
void analyze(art::Event const &e) override
Definition: ICARUSParticleAna_module.cc:106
ICARUSParticleAna::ICARUSParticleAna
ICARUSParticleAna(fhicl::ParameterSet const &p)
Definition: ICARUSParticleAna_module.cc:67
Generated by   doxygen 1.8.5