Validation_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: Validation
3 // Module Type: analyzer
4 // File: Validation_module.cc
5 //
6 // Generated at Mon Jun 8 16:52:56 2015 by John Matter using artmod
7 // from cetpkgsupport v1_08_05.
8 ////////////////////////////////////////////////////////////////////////
9 
17 #include "fhiclcpp/ParameterSet.h"
19 
21 #include <boost/algorithm/string.hpp>
22 #include "MCCheater/BackTracker.h"
23 #include "Simulation/Particle.h"
24 #include "RecoBase/CellHit.h"
25 #include <iostream>
26 #include <TTree.h>
27 
28 namespace mmc {
29  class Validation;
30 }
31 
33 public:
34  explicit Validation(fhicl::ParameterSet const & p);
35 
36  // Plugins should not be copied or assigned.
37  Validation(Validation const &) = delete;
38  Validation(Validation &&) = delete;
39  Validation & operator = (Validation const &) = delete;
40  Validation & operator = (Validation &&) = delete;
41 
42  // Required functions.
43  void analyze(art::Event const & e) override;
44  void beginJob() override;
45  void endJob() override;
46 
47 private:
51  std::map<std::string, double> et_, flsht_, cht_;
52 
53  template<class T>
54  void tset(T & tree, std::string const branch_name, double const value);
55 
56  double phi(simb::MCParticle const* particle) const;
57  double theta(simb::MCParticle const* particle) const;
58 
59  int GetMCTruth(const art::Event& evt,
60  std::vector<simb::MCTruth>& mctruth);
61 
62  void print_truth_info(art::Event const& e,
63  bool const& print_hit_details) const;
64 
66 };
67 
68 
70  EDAnalyzer(p),
71  branch_default_value_(-9999),
72  to_degrees_(180.0 / TMath::Pi())
73 {}
74 
76 {
77  //print_truth_info(e,true);
78 
79  // Initialize all branches to default value
80  for (auto & branch : et_)
82 
83  //-------------------------------------------------------
84  // MC info
85  if (bt_->HaveTruthInfo())
86  {
87  tset(et_, "is_mc", true);
88 
89  // First, get info from backtracker -------------------
91  tset(et_, "mc_bt_n_particles", pn.size());
92  if (pn.size() == 1)
93  {
94  sim::Particle const* bt_particle = pn.Primary(0);
95  tset(et_, "mc_bt_phi", phi(bt_particle));
96  tset(et_, "mc_bt_theta", theta(bt_particle));
97  tset(et_, "mc_bt_cos_theta", cos(theta(bt_particle)));
98  tset(et_, "mc_bt_beta", bt_particle->P(0) / bt_particle->Mass());
99  tset(et_, "mc_bt_n_daughters", bt_particle->NumberDaughters());
100 
101  // FLSHits ----------------------
102  double dx, dE, dEBirks, dEdx, dEBirksdx;
103  double Edep = 0;
104  double EBirksdep = 0;
105  double pathLength = 0;
106  const std::vector<sim::FLSHit>& flshits = bt_->ParticleToFLSHit(bt_particle->TrackId());
107  for (const auto flshit : flshits)
108  {
109  //Initialize branches
110  for (auto & branch : flsht_)
111  branch.second = branch_default_value_;
112 
113  // Get this hit's values
114  dx = flshit.GetTotalPathLength();
115  dE = flshit.GetEdep();
116  dEBirks = flshit.GetEdepBirks();
117  dEdx = dE / dx;
118  dEBirksdx = dEBirks / dx;
119 
120  tset(flsht_, "hit_dx", dx);
121  tset(flsht_, "hit_dE", dE);
122  tset(flsht_, "hit_dEBirks", dEBirks);
123  tset(flsht_, "hit_dEdx", dEdx);
124  tset(flsht_, "hit_dEBirksdx", dEBirksdx);
125 
126  tset(flsht_, "mc_bt_beta", bt_particle->P(0) / bt_particle->Mass());
127 
128  flshit_tree_->Fill();
129 
130  // Increment event-level totals
131  Edep += dE;
132  EBirksdep += dEBirks;
133  pathLength += dx;
134  } // loop over flshits
135  tset(et_, "mc_bt_dE", Edep);
136  tset(et_, "mc_bt_dEBirks", EBirksdep);
137  tset(et_, "mc_bt_pathlength", pathLength);
138 
139  // CellHits --------------------------------
141  e.getByLabel("calhit", hits);
142  std::vector<art::Ptr<rb::CellHit> > hitptrs;
143  for (unsigned i = 0; i < hits->size(); ++i)
144  hitptrs.emplace_back(hits, i);
145 
146  for (auto const& hitptr : hitptrs)
147  {
148  if (bt_->IsNoise(hitptr)) continue;
149 
150  //Initialize branches
151  for (auto & branch : cht_)
152  branch.second = branch_default_value_;
153 
154  tset(cht_, "view", hitptr->View());
155  tset(cht_, "plane", hitptr->Plane());
156  tset(cht_, "cell", hitptr->Cell());
157  tset(cht_, "tns", hitptr->TNS());
158  tset(cht_, "adc", hitptr->ADC()); //Rob: returns ADC(3)-ADC(0)
159 
160  tset(cht_, "mc_bt_beta", bt_particle->P(0) / bt_particle->Mass());
161 
162  cellhit_tree_->Fill();
163  }
164 
165  } else {
166  std::cerr << "\n\tParticle navigator contains "
167  << pn.size() << " particles. Check MCTruth info.\n"
168  << std::endl;
169  }
170 
171  // Then, get info from MCTruth a la Event Display -----
172  std::vector<simb::MCTruth> mctruth;
173  GetMCTruth(e, mctruth);
174  simb::MCParticle const mct_particle = mctruth[0].GetParticle(0); //unwise?
175  tset(et_, "mc_truth_phi", phi(&mct_particle));
176  tset(et_, "mc_truth_theta", theta(&mct_particle));
177  tset(et_, "mc_truth_cos_theta", cos(theta(&mct_particle)));
178  tset(et_, "mc_truth_beta", mct_particle.P(0) / mct_particle.Mass());
179 
180  } else {
181  tset(et_, "is_mc", false);
182  }
183  event_tree_->Fill();
184 }
185 
187 {
188  // FLSHit-Based Tree
189  std::vector<std::string> flshit_tree_names =
190  {"hit_dx",
191  "hit_dE",
192  "hit_dEBirks",
193  "hit_dEdx",
194  "hit_dEBirksdx",
195  "mc_bt_beta"};
196  flshit_tree_ = tfs_->make<TTree>("FLSHit", "Monopole FLSHit Tree");
197  for (auto const& name : flshit_tree_names)
198  flshit_tree_->Branch(name.c_str(), &flsht_[name], (name + "/D").c_str());
199 
200  // CellHit-Based Tree
201  std::vector<std::string> cellhit_tree_names =
202  {"view",
203  "plane",
204  "cell",
205  "tns",
206  "adc",
207  "mc_bt_beta"};
208  cellhit_tree_ = tfs_->make<TTree>("CellHit", "Monopole CellHit Tree");
209  for (auto const& name : cellhit_tree_names)
210  cellhit_tree_->Branch(name.c_str(), &cht_[name], (name + "/D").c_str());
211 
212  // Event-Based Tree
213  std::vector<std::string> event_tree_names =
214  {"is_mc",
215  "mc_bt_phi",
216  "mc_bt_theta",
217  "mc_bt_cos_theta",
218  "mc_bt_beta",
219  "mc_bt_n_particles",
220  "mc_bt_n_daughters",
221  "mc_bt_pathlength",
222  "mc_bt_dE",
223  "mc_bt_dEBirks",
224  "mc_truth_phi",
225  "mc_truth_theta",
226  "mc_truth_cos_theta",
227  "mc_truth_beta"};
228  event_tree_ = tfs_->make<TTree>("Event", "Monopole Event Tree");
229  for (auto const& name : event_tree_names)
230  event_tree_->Branch(name.c_str(), &et_[name], (name + "/D").c_str());
231 }
232 
234 {
235 }
236 
237 template<class T>
239 (T & tree, std::string const branch_name, double const value)
240 {
241  auto it = tree.find(branch_name);
242  if (it == tree.end())
243  {
244  std::cerr << "\n\t Branch " << branch_name << " does not exist! \n"
245  << std::endl;
246  assert(false);
247  }
248  it->second = value;
249 }
250 
251 double mmc::Validation::phi(simb::MCParticle const* particle) const
252 {
253  double Px(particle->Px()), Py(particle->Py());
254  return std::atan2(Py, Px);
255 }
256 
257 double mmc::Validation::theta(simb::MCParticle const* particle) const
258 {
259  double Px(particle->Px()), Py(particle->Py()), Pz(particle->Pz());
260  double Pt = std::sqrt(Px * Px + Py * Py);
261  return std::atan2(Pt, Pz);
262 }
263 
264 //GetMCTruth was appropriated from EventDisplay/SimulationDrawer
266  std::vector<simb::MCTruth>& mctruth)
267 {
268  mctruth.clear();
269 
270  // going to use getManyByType so that we don't have to
271  // worry about using a module label
272  std::vector< art::Handle< std::vector<simb::MCTruth> > > tcol;
273 
274  evt.getManyByType(tcol);
275  for(size_t mctc = 0; mctc < tcol.size(); ++mctc)
276  {
277  for(size_t mct = 0; mct < tcol[mctc]->size(); ++mct)
278  {
279  mctruth.push_back(tcol[mctc]->at(mct));
280  }
281  }
282 
283  return mctruth.size();
284 }
285 
286 
288 (art::Event const& e, bool const& print_hit_details) const
289 {
290  if (bt_->HaveTruthInfo())
291  {
293  unsigned n_particle = 0;
294 
295  for (auto p_handle : pn)
296  {
297  const sim::Particle* particle = p_handle.second;
298  std::cout << "MF: MC Particle " << ++n_particle
299  << "\tEvent " << e.event()
300  << "\nMF: Mass = " << particle->Mass()
301  << "\nMF: N Trajectory Points = "
302  << particle->NumberTrajectoryPoints()
303  << std::endl;
304 
305  for (unsigned trajectory = 0;
306  trajectory != particle->NumberTrajectoryPoints(); ++trajectory)
307  {
308  double to_degrees = 180 / TMath::Pi();
309  std::cout << "MF: Trajectory " << trajectory
310  << "\nMF:\t Energy = " << particle->E(trajectory)
311  << "\nMF:\t Momentum = " << particle->P(trajectory)
312  << "\nMF:\t Beta = "
313  << particle->P(trajectory) / particle->Mass()
314  << "\nMF:\t Time = " << particle->T(trajectory)
315  << "\nMF:\t Phi = " << to_degrees * phi(particle)
316  << "\nMF:\t Theta = " << to_degrees * theta(particle)
317  << "\nMF:\t [Vx, Vy, Vz] = ["
318  << particle->Vx(trajectory) << ", "
319  << particle->Vy(trajectory) << ", "
320  << particle->Vz(trajectory) << "]"
321  << "\nMF:\t [Px, Py, Pz] = ["
322  << particle->Px(trajectory) << ", "
323  << particle->Py(trajectory) << ", "
324  << particle->Pz(trajectory) << "]"
325  << std::endl;
326  }
327 
328  if (print_hit_details)
329  {
330  std::vector<sim::FLSHit> fls_hits =
331  bt_->ParticleToFLSHit(particle->TrackId());
332  std::cout << "MF: Hits:" << std::endl;
333  std::sort(fls_hits.begin(), fls_hits.end(),
334  [](sim::FLSHit const& lhs, sim::FLSHit const& rhs)
335  { return lhs.GetPlaneID() < rhs.GetPlaneID(); });
336 
337  for (auto const& fls_hit : fls_hits)
338  {
339  float dt = fls_hit.GetExitT() - fls_hit.GetEntryT();
340  float dx = fls_hit.GetTotalPathLength();
341  float v = dx / dt * 1.0e7;
342  float beta = v / 3.0e8;
343  float dEdx = fls_hit.GetEdep() / dx;
344 
345  std::cout << "MF:\t [plane, cell, beta, dEdx, path, Edep] = ["
346  << fls_hit.GetPlaneID() << ", "
347  << fls_hit.GetCellID() << ", "
348  << beta << ", "
349  << dEdx << ", "
350  << fls_hit.GetTotalPathLength() << ", "
351  << fls_hit.GetEdep() << "]" << std::endl;
352 
353  for (int point = 0; point != fls_hit.GetNPoints(); ++point)
354  {
355  std::cout << "MF:\t\t [x, y, z] = ["
356  << fls_hit.GetX(point) << ", "
357  << fls_hit.GetY(point) << ", "
358  << fls_hit.GetZ(point) << "]" << std::endl;
359  }
360  }
361  }
362  }
363  }
364 }
365 
double E(const int i=0) const
Definition: MCParticle.h:232
const XML_Char * name
Definition: expat.h:151
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:217
back track the reconstruction to the simulation
std::map< std::string, double > flsht_
Validation(fhicl::ParameterSet const &p)
double Py(const int i=0) const
Definition: MCParticle.h:230
set< int >::iterator it
int GetPlaneID() const
Plane ID.
Definition: FLSHit.h:37
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
size_type size() const
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
std::vector< sim::FLSHit > ParticleToFLSHit(const int &trackID) const
All the FLSHits that were created by the track id trackID, sorted from most to least light...
double Mass() const
Definition: MCParticle.h:238
double Px(const int i=0) const
Definition: MCParticle.h:229
Validation & operator=(Validation const &)=delete
OStream cerr
Definition: OStream.cxx:7
A single unit of energy deposition in the liquid scintillator.
Definition: FLSHit.h:19
Double_t beta
void analyze(art::Event const &e) override
DEFINE_ART_MODULE(TestTMapFile)
double dE
int NumberDaughters() const
Definition: MCParticle.h:216
bool IsNoise(const art::Ptr< rb::CellHit > &hit) const
Is this hit not associated with any particles?
int TrackId() const
Definition: MCParticle.h:209
void beginJob() override
Definition: Cand.cxx:23
void tset(T &tree, std::string const branch_name, double const value)
void hits()
Definition: readHits.C:15
double dx[NP][NC]
const sim::Particle * Primary(const int) const
const XML_Char int const XML_Char * value
Definition: expat.h:331
double P(const int i=0) const
Definition: MCParticle.h:233
double theta(simb::MCParticle const *particle) const
int evt
double T(const int i=0) const
Definition: MCParticle.h:223
void getManyByType(std::vector< Handle< PROD >> &results) const
Definition: DataViewImpl.h:446
void endJob() override
art::ServiceHandle< art::TFileService > tfs_
EventNumber_t event() const
Definition: Event.h:67
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
void print_truth_info(art::Event const &e, bool const &print_hit_details) const
const double branch_default_value_
OStream cout
Definition: OStream.cxx:6
double Vx(const int i=0) const
Definition: MCParticle.h:220
T * make(ARGS...args) const
const double to_degrees_
std::map< std::string, double > et_
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:133
double Pz(const int i=0) const
Definition: MCParticle.h:231
double Vz(const int i=0) const
Definition: MCParticle.h:222
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
double T
Definition: Xdiff_gwt.C:5
Float_t e
Definition: plot.C:35
std::map< std::string, double > cht_
double phi(simb::MCParticle const *particle) const
art::ServiceHandle< cheat::BackTracker > bt_
int GetMCTruth(const art::Event &evt, std::vector< simb::MCTruth > &mctruth)
double Vy(const int i=0) const
Definition: MCParticle.h:221
T atan2(T number)
Definition: d0nt_math.hpp:72
enum BeamMode string