HadCal_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: HadCal
3 // Module Type: analyzer
4 // File: HadCal_module.cc
5 //
6 // Generated at Wed Feb 8 09:38:58 2017 by Dorota Stefan using artmod
7 // from cetpkgsupport v1_11_00.
8 ////////////////////////////////////////////////////////////////////////
9 
10 #include "lardataobj/Simulation/SimChannel.h"
11 #include "larsim/Simulation/LArG4Parameters.h"
12 #include "larsim/Simulation/LArVoxelData.h"
13 #include "larsim/Simulation/LArVoxelList.h"
14 #include "larsim/Simulation/SimListUtils.h"
15 #include "larcore/Geometry/Geometry.h"
16 #include "larcorealg/Geometry/GeometryCore.h"
17 #include "lardataobj/RecoBase/Hit.h"
18 #include "lardataobj/RecoBase/Cluster.h"
19 #include "lardataobj/RecoBase/Track.h"
20 #include "lardataobj/RecoBase/PFParticle.h"
21 #include "lardataobj/AnalysisBase/Calorimetry.h"
22 #include "larreco/Calorimetry/CalorimetryAlg.h"
23 #include "lardataobj/RecoBase/TrackHitMeta.h"
24 #include "nusimdata/SimulationBase/MCParticle.h"
25 #include "nusimdata/SimulationBase/MCTruth.h"
26 #include "larcoreobj/SimpleTypesAndConstants/PhysicalConstants.h"
27 #include "lardata/Utilities/DatabaseUtil.h"
28 #include "lardata/DetectorInfoServices/DetectorClocksService.h"
29 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
30 
31 #include "art/Framework/Core/EDAnalyzer.h"
32 #include "art/Framework/Core/ModuleMacros.h"
33 #include "canvas/Persistency/Common/FindManyP.h"
34 #include "art/Framework/Principal/Event.h"
35 #include "art/Framework/Principal/Handle.h"
36 #include "art/Framework/Principal/Run.h"
37 #include "art/Framework/Principal/SubRun.h"
38 #include "art/Framework/Services/Registry/ServiceHandle.h"
39 #include "art_root_io/TFileService.h"
40 #include "canvas/Utilities/InputTag.h"
41 #include "fhiclcpp/ParameterSet.h"
42 #include "messagefacility/MessageLogger/MessageLogger.h"
43 
44 #include "lardata/ArtDataHelper/MVAReader.h"
45 
46 #include "TH1.h"
47 #include "TTree.h"
48 #include "TLorentzVector.h"
49 #include "TVector3.h"
50 
51 #include <fstream>
52 
53 #define MVA_LENGTH 4
54 
55 namespace proto
56 {
57  class HadCal;
58 }
59 
60 class proto::HadCal : public art::EDAnalyzer {
61 public:
62  explicit HadCal(fhicl::ParameterSet const & p);
63  // The destructor generated by the compiler is fine for classes
64  // without bare pointers or other resource use.
65 
66  // Plugins should not be copied or assigned.
67  HadCal(HadCal const &) = delete;
68  HadCal(HadCal &&) = delete;
69  HadCal & operator = (HadCal const &) = delete;
70  HadCal & operator = (HadCal &&) = delete;
71 
72  // Required functions.
73  void analyze(art::Event const & e) override;
74 
75  void beginJob() override;
76  void endJob() override;
77 
78  void reconfigure(fhicl::ParameterSet const& p) ;
79 
80 private:
81 
82  // Declare member data here.
83  double GetEdepHitsMeV(detinfo::DetectorClocksData const& clockData,
84  detinfo::DetectorPropertiesData const& detProp,
85  const std::vector< recob::Hit > & hits) const;
86  double GetEhitMeV(const recob::Hit & hit) const;
87  double GetEhitMeVcorrel(detinfo::DetectorClocksData const& clockData,
88  detinfo::DetectorPropertiesData const& detProp,
89  const recob::Hit & hit) const;
90  double GetEkinMeV(detinfo::DetectorClocksData const& clockData,
91  detinfo::DetectorPropertiesData const& detProp,
92  const std::vector < art::Ptr< recob::Hit > > & hits,
93  const std::vector < recob::TrackHitMeta const* > & data);
94  // EM
95  double GetEdepEM_MC(art::Event const & e) const;
96  double GetEdepEMAtt_MC(art::Event const & e) const;
97  double GetEdepEMhAtt_MC(art::Event const & e,
98  detinfo::DetectorClocksData const& clockData,
99  detinfo::DetectorPropertiesData const& detProp,
100  const std::vector< recob::Hit > & hits) const;
101 
102  double fEdepEM_MC;
103  double fEdepEMAtt_MC;
104  double fEdepEMhAtt_MC; // new
105 
106  double fEdepEM_reco;
107  double fEdepEMcorr_ellifetime_reco;
108 
109  // HAD
110  double GetEdepHAD_MC(art::Event const & e) const;
111  double GetEdepHADAtt_MC(art::Event const & e) const;
112  double GetEdepHADhAtt_MC(art::Event const & e,
113  detinfo::DetectorClocksData const& clockData,
114  detinfo::DetectorPropertiesData const& detProp,
115  const std::vector< recob::Hit > & hits) const;
116 
117  double fEdepHAD_MC; // new
118  double fEdepHADAtt_MC; // new
119  double fEdepHADhAtt_MC;
120 
121  double fEdepHAD_reco;
122  double fEdepHADcorr_ellifetime_reco;
123 
124  // MISSING ENERGY
125  double fMissEn_MC;
126 
127  void ResetVars();
128 
129  geo::GeometryCore const * fGeometry;
130 
131  TTree *fTree;
132  TTree *fTreeEntries;
133  int fRun;
134  int fEvent;
135  int fBestView;
136  int fNumberOfTracks;
137 
138 
139  double fdQdx;
140  double fdEdx;
141  double fdQ;
142  double fdx;
143  double fEdepSum;
144  double fEdepAllhits;
145  double fElectronsToGeV;
146  double fResRange;
147  double fEnGen;
148  double fEkGen;
149  double fT0;
150 
151  std::ofstream file;
152 
153  // Module labels to get data products
154  art::InputTag fNNetModuleLabel;
155  std::string fHitModuleLabel;
156  std::string fClusterModuleLabel;
157  std::string fTrackModuleLabel;
158  std::string fCalorimetryModuleLabel;
159  std::string fSimulationLabel;
160  calo::CalorimetryAlg fCalorimetryAlg;
161 
162  std::unordered_map< int, const simb::MCParticle* > fParticleMap;
163 };
164 
165 
166 proto::HadCal::HadCal(fhicl::ParameterSet const & p)
167  :
168  EDAnalyzer(p),
169  fCalorimetryAlg(p.get<fhicl::ParameterSet>("CalorimetryAlg"))
170  // More initializers here.
171 {
172  reconfigure(p);
173  // get a pointer to the geometry service provider
174  fGeometry = &*(art::ServiceHandle<geo::Geometry>());
175  file.open("dump.txt");
176 }
177 
178 void proto::HadCal::beginJob()
179 {
180  art::ServiceHandle<sim::LArG4Parameters> larParameters;
181  fElectronsToGeV = 1./larParameters->GeVToElectrons();
182 
183  // access art's TFileService, which will handle creating and writing hists
184  art::ServiceHandle<art::TFileService> tfs;
185 
186  fTree = tfs->make<TTree>("calibration","calibration tree");
187  fTree->Branch("fRun", &fRun, "fRun/I");
188  fTree->Branch("fEvent", &fEvent, "fEvent/I");
189  fTree->Branch("fEnGen", &fEnGen, "fEnGen/D");
190  fTree->Branch("fEkGen", &fEkGen, "fEkGeb/D");
191 
192  fTree->Branch("fEdepEM_MC", &fEdepEM_MC, "fEdepEM_MC/D");
193  fTree->Branch("fEdepEMAtt_MC", &fEdepEMAtt_MC, "fEdepEMAtt_MC/D");
194  fTree->Branch("fEdepEMhAtt_MC", &fEdepEMhAtt_MC, "fEdepEMhAtt_MC/D");
195 
196  fTree->Branch("fEdepEM_reco", &fEdepEM_reco, "fEdepEM_reco/D");
197  fTree->Branch("fEdepEMcorr_ellifetime_reco", &fEdepEMcorr_ellifetime_reco, "fEdepEMcorr_ellifetime_reco/D");
198 
199  fTree->Branch("fEdepHAD_MC", &fEdepHAD_MC, "fEdepHAD_MC/D");
200  fTree->Branch("fEdepHADAtt_MC", &fEdepHADAtt_MC, "fEdepHADAtt_MC/D");
201  fTree->Branch("fEdepHADhAtt_MC", &fEdepHADhAtt_MC, "fEdepHADhAtt_MC/D");
202 
203  fTree->Branch("fEdepHAD_reco", &fEdepHAD_reco, "fEdepHAD_reco/D");
204  fTree->Branch("fEdepHADcorr_ellifetime_reco", &fEdepHADcorr_ellifetime_reco, "fEdepHADcorr_ellifetime_reco/D");
205 
206  fTree->Branch("fMissEn_MC", &fMissEn_MC, "fMissEn_MC/D");
207 
208  fTree->Branch("fNumberOfTracks", &fNumberOfTracks, "fNumberOfTracks/I");
209 
210  fTree->Branch("fEdepSum", &fEdepSum, "fEdepSum/D");
211  fTree->Branch("fEdepAllhits", &fEdepAllhits, "fEdepAllhits/D");
212  fTree->Branch("fT0", &fT0, "fT0/D");
213 
214  fTreeEntries = tfs->make<TTree>("entries","entries tree");
215  fTreeEntries->Branch("fdQdx", &fdQdx, "fdQdx/D");
216  fTreeEntries->Branch("fdEdx", &fdEdx, "fdEdx/D");
217  fTreeEntries->Branch("fdQ", &fdQ, "fdQ/D");
218  fTreeEntries->Branch("fdx", &fdx, "fdx/D");
219 
220 }
221 
222 void proto::HadCal::endJob()
223 {
224  file.close();
225 }
226 
227 void proto::HadCal::reconfigure(fhicl::ParameterSet const & p)
228 {
229  fHitModuleLabel = p.get< std::string >("HitModuleLabel");
230  fClusterModuleLabel = p.get< std::string >("ClusterModuleLabel");
231  fTrackModuleLabel = p.get< std::string >("TrackModuleLabel");
232  fCalorimetryModuleLabel = p.get< std::string >("CalorimetryModuleLabel");
233  fSimulationLabel = p.get< std::string >("SimulationLabel");
234  fNNetModuleLabel = p.get< std::string >("NNetModuleLabel");
235 }
236 
237 void proto::HadCal::analyze(art::Event const & e)
238 {
239  // Implementation of required member function here.
240  ResetVars();
241 
242  fRun = e.run();
243  fEvent = e.id().event();
244 
245  // MC particle list
246  auto particleHandle = e.getValidHandle< std::vector<simb::MCParticle> >(fSimulationLabel);
247  bool flag = true;
248 
249  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService>()->DataFor(e);
250  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService>()->DataFor(e, clockData);
251 
252  for (auto const& p : *particleHandle)
253  {
254  fParticleMap[p.TrackId()] = &p;
255  if ((p.Process() == "primary") && flag)
256  {
257  fEnGen = p.P();
258  fEkGen = (std::sqrt(p.P()*p.P() + p.Mass()*p.Mass()) - p.Mass()) * 1000; // MeV
259  fT0 = p.T();
260  flag = false;
261  }
262  }
263 
264  fEdepEM_MC = GetEdepEM_MC(e);
265  fEdepEMAtt_MC = GetEdepEMAtt_MC(e);
266 
267  fEdepHAD_MC = GetEdepHAD_MC(e);
268  fEdepHADAtt_MC = GetEdepHADAtt_MC(e);
269 
270  // hits
271  const auto& hitListHandle = *e.getValidHandle< std::vector<recob::Hit> >(fHitModuleLabel);
272  fEdepAllhits = GetEdepHitsMeV(clockData, detProp, hitListHandle);
273 
274  // MC associated with a hit
275  fEdepEMhAtt_MC = GetEdepEMhAtt_MC(e, clockData, detProp, hitListHandle);
276  fEdepHADhAtt_MC = GetEdepHADhAtt_MC(e, clockData, detProp, hitListHandle);
277 
278  // MC missing energy
279  fMissEn_MC = fEkGen - (fEdepEM_MC + fEdepHAD_MC);
280 
281  // clusters
282  auto cluResults = anab::MVAReader< recob::Cluster, MVA_LENGTH>::create(e, fNNetModuleLabel);
283 
284  std::unordered_map<int, bool> hitIDE;
285  if (cluResults)
286  {
287  const art::FindManyP< recob::Hit > hitsFromCluster(cluResults->dataHandle(), e, cluResults->dataTag());
288 
289  for (size_t c = 0; c < cluResults->size(); ++c)
290  {
291  for (size_t h = 0; h < hitsFromCluster.at(c).size(); ++h)
292  {
293  if (hitsFromCluster.at(c)[h]->View() == fBestView)
294  {
295  hitIDE[hitsFromCluster.at(c)[h].key()] = false;
296  }
297  }
298  }
299  }
300 
301  // tracks
302  fNumberOfTracks = 0;
303  auto trkHandle = e.getValidHandle< std::vector<recob::Track> >(fTrackModuleLabel);
304  art::FindManyP< recob::PFParticle > pfpFromTrack(trkHandle, e, fTrackModuleLabel);
305 
306  const art::FindManyP< anab::Calorimetry > calFromTracks(trkHandle, e, fCalorimetryModuleLabel);
307  const art::FindManyP< recob::Hit > hitsFromTracks(trkHandle, e, fTrackModuleLabel);
308  const art::FindManyP<recob::Hit, recob::TrackHitMeta> fmthm(trkHandle, e, fTrackModuleLabel);
309 
310  if (fmthm.isValid())
311  {
312  // loop over tracks
313  for (size_t t = 0; t < trkHandle->size(); ++t)
314  {
315  auto vhit = fmthm.at(t);
316  auto vmeta = fmthm.data(t);
317 
318  // mark that hits have been used
319  for (size_t h = 0; h < hitsFromTracks.at(t).size(); ++h)
320  {
321  if (hitsFromTracks.at(t)[h]->View() == fBestView)
322  {
323  hitIDE[hitsFromTracks.at(t)[h].key()] = true;
324  }
325  }
326 
327  int nplanes = calFromTracks.at(t).size();
328 
329  for (size_t p = 0; p < pfpFromTrack.at(t).size(); ++p)
330  {
331  int pdg = pfpFromTrack.at(t)[p]->PdgCode();
332  // condition for tracks
333  if ((pdg != 11) && (pdg != -11))
334  {
335  fNumberOfTracks++;
336 
337  // for now it work for 3 planes and only collection view
338  if (nplanes == 3)
339  {
340  // fHadEnSum += GetEkinMeV(vhit, vmeta);
341  for (size_t h = 0; h < hitsFromTracks.at(t).size(); ++h)
342  {
343  if (hitsFromTracks.at(t)[h]->View() == fBestView)
344  {
345  fEdepHAD_reco += GetEhitMeV(*hitsFromTracks.at(t)[h]);
346  fEdepHADcorr_ellifetime_reco += GetEhitMeVcorrel(clockData, detProp, *hitsFromTracks.at(t)[h]);
347  }
348  }
349  }
350  }
351  else // ... and condition for em showers
352  {
353  for (size_t h = 0; h < hitsFromTracks.at(t).size(); ++h)
354  {
355  // kinetic energy without correction for recombination,
356  // only one view is considered
357  if (hitsFromTracks.at(t)[h]->View() == fBestView)
358  {
359  fEdepEM_reco += GetEhitMeV(*hitsFromTracks.at(t)[h]);
360  fEdepEMcorr_ellifetime_reco += GetEhitMeVcorrel(clockData, detProp, *hitsFromTracks.at(t)[h]);
361  }
362  }
363  }
364  }
365  }
366  }
367 
368  // output from cnn's
369  if (cluResults)
370  {
371  const art::FindManyP< recob::Hit > hitsFromCluster(cluResults->dataHandle(), e, cluResults->dataTag());
372 
373  // loop over clusters
374  for (size_t c = 0; c < cluResults->size(); ++c)
375  {
376  if (cluResults->item(c).View() == fBestView)
377  {
378  std::array< float, MVA_LENGTH > cnn_out = cluResults->getOutput(c);
379  double p_trk_or_sh = cnn_out[0] + cnn_out[1];
380  double pdg = 1;
381  if (p_trk_or_sh > 0) pdg = cnn_out[1] / p_trk_or_sh;
382 
383  for (size_t h = 0; h < hitsFromCluster.at(c).size(); ++h)
384  {
385  if (hitIDE[hitsFromCluster.at(c)[h].key()] == false)
386  {
387  hitIDE[hitsFromCluster.at(c)[h].key()] = true;
388 
389  if (pdg < 0.63)
390  {
391  fEdepHAD_reco += GetEhitMeV(*hitsFromCluster.at(c)[h]);
392  fEdepHADcorr_ellifetime_reco += GetEhitMeVcorrel(clockData, detProp, *hitsFromCluster.at(c)[h]);
393  }
394  else
395  {
396  fEdepEM_reco += GetEhitMeV(*hitsFromCluster.at(c)[h]);
397  fEdepEMcorr_ellifetime_reco += GetEhitMeVcorrel(clockData, detProp, *hitsFromCluster.at(c)[h]);
398  }
399  }
400 
401  fEdepSum += GetEhitMeV(*hitsFromCluster.at(c)[h]);
402  }
403  }
404  }
405  }
406 
407  /*file << fEvent << " " << fEkGen << " " << fEdepEM_MC << " " << fEdepHAD_MC << " " << fMissEn_MC
408  << " " << fEdepEMAtt_MC << " " << fEdepEMhAtt_MC
409  << " " << fEdepHADAtt_MC << " " << fEdepHADhAtt_MC
410  << " " << fEdepHADcorr_ellifetime_reco << " " << fEdepEMcorr_ellifetime_reco << std::endl;*/
411  fTree->Fill();
412 }
413 
414 // kinetic energy of hadronic part: correction for recombination applied according to Birks/Box model
415 
416 /******************************************************************************/
417 
418 double proto::HadCal::GetEkinMeV(detinfo::DetectorClocksData const& clockData,
419  detinfo::DetectorPropertiesData const& detProp,
420  const std::vector < art::Ptr< recob::Hit > > & hits,
421  const std::vector < recob::TrackHitMeta const* > & data)
422 {
423  double ekin = 0.0;
424 
425  if (!hits.size()) return ekin;
426 
427 
428  for (size_t h = 0; h < hits.size(); ++h)
429  {
430  double dqadc = hits[h]->Integral();
431  if (!std::isnormal(dqadc) || (dqadc < 0)) dqadc = 0.0;
432 
433  unsigned short plane = hits[h]->WireID().Plane;
434  unsigned short time = hits[h]->PeakTime();
435  double t0 = 0;
436 
437  fdQdx = 0.0;
438  fdQ = dqadc;
439  fdx = data[h]->Dx();
440  if ((fdx > 0) && (fdQ > 0))
441  {
442  fdQdx = fdQ/fdx;
443  fdEdx = fCalorimetryAlg.dEdx_AREA(clockData, detProp, fdQdx, time, plane, t0);
444  if (fdEdx > 35) fdEdx = 35;
445 
446  ekin += ( fdEdx * fdx );
447  }
448  else if ((fdx == 0) && (fdQ > 0))
449  {
450 
451  }
452 
453  fTreeEntries->Fill();
454  }
455 
456  return ekin;
457 }
458 
459 /******************************************************************************/
460 
461 double proto::HadCal::GetEdepHitsMeV(detinfo::DetectorClocksData const& clockData,
462  detinfo::DetectorPropertiesData const& detProp,
463  const std::vector< recob::Hit > & hits) const
464 {
465  if (!hits.size()) return 0.0;
466 
467  double dqsum = 0.0;
468  for (size_t h = 0; h < hits.size(); ++h)
469  {
470  double dqadc = hits[h].Integral();
471  if (!std::isnormal(dqadc) || (dqadc < 0)) continue;
472 
473  unsigned short plane = hits[h].WireID().Plane;
474  if (plane == fBestView)
475  {
476  double tdrift = hits[h].PeakTime();
477  double dqel = fCalorimetryAlg.ElectronsFromADCArea(dqadc, plane);
478 
479  double correllifetime = fCalorimetryAlg.LifetimeCorrection(clockData, detProp, tdrift, fT0);
480  double dq = dqel * correllifetime * fElectronsToGeV * 1000;
481 
482  if (!std::isnormal(dq) || (dq < 0)) continue;
483 
484  dqsum += dq;
485  }
486  }
487 
488  return dqsum;
489 }
490 
491 /******************************************************************************/
492 
493 double proto::HadCal::GetEhitMeVcorrel(detinfo::DetectorClocksData const& clockData,
494  detinfo::DetectorPropertiesData const& detProp,
495  const recob::Hit & hit) const
496 {
497  double dqadc = hit.Integral();
498  if (!std::isnormal(dqadc) || (dqadc < 0)) dqadc = 0.0;
499 
500  unsigned short plane = hit.WireID().Plane;
501  double tdrift = hit.PeakTime();
502  double dqel = fCalorimetryAlg.ElectronsFromADCArea(dqadc, plane);
503 
504  double correllifetime = fCalorimetryAlg.LifetimeCorrection(clockData, detProp, tdrift, fT0);
505  double dq = dqel * correllifetime * fElectronsToGeV * 1000;
506  if (!std::isnormal(dq) || (dq < 0)) dq = 0.0;
507 
508  return dq;
509 }
510 
511 /******************************************************************************/
512 
513 double proto::HadCal::GetEhitMeV(const recob::Hit & hit) const
514 {
515  double dqadc = hit.Integral();
516  if (!std::isnormal(dqadc) || (dqadc < 0)) dqadc = 0.0;
517 
518  unsigned short plane = hit.WireID().Plane;
519  // double tdrift = hit.PeakTime();
520  double dqel = fCalorimetryAlg.ElectronsFromADCArea(dqadc, plane);
521 
522  double correllifetime = 1; // fCalorimetryAlg.LifetimeCorrection(tdrift, fT0);
523  double dq = dqel * correllifetime * fElectronsToGeV * 1000;
524  if (!std::isnormal(dq) || (dq < 0)) dq = 0.0;
525 
526  return dq;
527 }
528 
529 /******************************************************************************/
530 
531 double proto::HadCal::GetEdepEM_MC(art::Event const & e) const
532 {
533  double enEM = 0.0;
534 
535  auto simchannelHandle = e.getHandle< std::vector<sim::SimChannel> >(fSimulationLabel);
536  if (simchannelHandle)
537  {
538  for ( auto const& channel : (*simchannelHandle) )
539  {
540  if (fGeometry->View(channel.Channel()) == fBestView)
541  {
542  // for every time slice in this channel:
543  auto const& timeSlices = channel.TDCIDEMap();
544  for ( auto const& timeSlice : timeSlices )
545  {
546  // loop over the energy deposits.
547  auto const& energyDeposits = timeSlice.second;
548 
549  for ( auto const& energyDeposit : energyDeposits )
550  {
551  double energy = energyDeposit.energy; // energy not attenuated
552  int trackID = energyDeposit.trackID;
553 
554  if (trackID < 0)
555  {
556  enEM += energy;
557  }
558  else if (trackID > 0)
559  {
560  auto search = fParticleMap.find(trackID);
561  bool found = true;
562  if (search == fParticleMap.end())
563  {
564  mf::LogWarning("TrainingDataAlg") << "PARTICLE NOT FOUND";
565  found = false;
566  }
567 
568  int pdg = 0;
569  if (found)
570  {
571  const simb::MCParticle& particle = *((*search).second);
572  if (!pdg) pdg = particle.PdgCode(); // not EM activity so read what PDG it is
573  }
574 
575  if ((pdg == 11) || (pdg == -11) || (pdg == 22)) enEM += energy;
576  }
577 
578  }
579  }
580  }
581  }
582  }
583 
584  return enEM;
585 }
586 
587 /******************************************************************************/
588 
589 double proto::HadCal::GetEdepEMAtt_MC(art::Event const & e) const
590 {
591  double enEM = 0.0;
592 
593  auto simchannelHandle = e.getHandle< std::vector<sim::SimChannel> >(fSimulationLabel);
594  if (simchannelHandle)
595  {
596  for ( auto const& channel : (*simchannelHandle) )
597  {
598  if (fGeometry->View(channel.Channel()) == fBestView)
599  {
600  // for every time slice in this channel:
601  auto const& timeSlices = channel.TDCIDEMap();
602  for ( auto const& timeSlice : timeSlices )
603  {
604  // loop over the energy deposits.
605  auto const& energyDeposits = timeSlice.second;
606 
607  for ( auto const& energyDeposit : energyDeposits )
608  {
609  double energy = energyDeposit.numElectrons * fElectronsToGeV * 1000; // attenuated
610  int trackID = energyDeposit.trackID;
611 
612  if (trackID < 0)
613  {
614  enEM += energy;
615  }
616  else if (trackID > 0)
617  {
618  auto search = fParticleMap.find(trackID);
619  bool found = true;
620  if (search == fParticleMap.end())
621  {
622  mf::LogWarning("TrainingDataAlg") << "PARTICLE NOT FOUND";
623  found = false;
624  }
625 
626  int pdg = 0;
627  if (found)
628  {
629  const simb::MCParticle& particle = *((*search).second);
630  if (!pdg) pdg = particle.PdgCode(); // not EM activity so read what PDG it is
631  }
632 
633  if ((pdg == 11) || (pdg == -11) || (pdg == 22)) enEM += energy;
634  }
635 
636  }
637  }
638  }
639  }
640  }
641  return enEM;
642 }
643 
644 
645 /******************************************************************************/
646 
647 double proto::HadCal::GetEdepEMhAtt_MC(art::Event const & e,
648  detinfo::DetectorClocksData const& clockData,
649  detinfo::DetectorPropertiesData const& detProp,
650  const std::vector< recob::Hit > & hits) const
651 {
652  double totemhit = 0.0;
653 
654  auto simChannelHandle = e.getValidHandle< std::vector<sim::SimChannel> >(fSimulationLabel);
655  const std::vector< sim::SimChannel > & channels = *simChannelHandle;
656 
657  for (auto const & hit: hits)
658  {
659  // the channel associate with this hit.
660  auto hitChannelNumber = hit.Channel();
661 
662  double hitEn = 0.0; double hitEnSh = 0;
663 
664  for (auto const & channel: channels)
665  {
666  if (channel.Channel() != hitChannelNumber) continue;
667  if (fGeometry->View(channel.Channel()) != fBestView) continue;
668 
669  // for every time slice in this channel:
670  auto const& timeSlices = channel.TDCIDEMap();
671  for ( auto const& timeSlice : timeSlices )
672  {
673  int time = timeSlice.first;
674  if (std::abs(hit.TimeDistanceAsRMS(time)) < 1.0)
675  {
676  // loop over the energy deposits.
677  auto const& energyDeposits = timeSlice.second;
678 
679  for ( auto const& energyDeposit : energyDeposits )
680  {
681  int trackID = energyDeposit.trackID;
682  double energy = energyDeposit.numElectrons * fElectronsToGeV * 1000; // attenuated
683  hitEn += energy;
684 
685  if (trackID < 0)
686  {
687  hitEnSh += energy;
688  }
689  else if (trackID > 0)
690  {
691  auto search = fParticleMap.find(trackID);
692  bool found = true;
693  if (search == fParticleMap.end())
694  {
695  mf::LogWarning("TrainingDataAlg") << "PARTICLE NOT FOUND";
696  found = false;
697  }
698 
699  int pdg = 0;
700  if (found)
701  {
702  const simb::MCParticle& particle = *((*search).second);
703  if (!pdg) pdg = particle.PdgCode(); // not EM activity so read what PDG it is
704  }
705 
706  if ((pdg == 11) || (pdg == -11) || (pdg == 22)) hitEnSh += energy;
707  }
708  }
709  }
710  }
711  }
712 
713  double ratio = 0.0;
714  if (hitEn > 0)
715  {
716  ratio = hitEnSh / hitEn;
717  }
718 
719  // use energy of hit corrected for electron lifetime
720  if (ratio > 0.5)
721  {
722  totemhit += GetEhitMeVcorrel(clockData, detProp, hit);
723  }
724  }
725 
726  return totemhit;
727 }
728 
729 /******************************************************************************/
730 
731 double proto::HadCal::GetEdepHAD_MC(art::Event const & e) const
732 {
733  double enHAD = 0.0;
734 
735  auto simchannelHandle = e.getHandle< std::vector<sim::SimChannel> >(fSimulationLabel);
736  if (simchannelHandle)
737  {
738 
739  for (auto const & channel : (*simchannelHandle) )
740  {
741  if (fGeometry->View(channel.Channel()) != fBestView) continue;
742 
743  // for every time slice in this channel:
744  auto const& timeSlices = channel.TDCIDEMap();
745  for (auto const& timeSlice : timeSlices)
746  {
747  // loop over the energy deposits.
748  auto const & energyDeposits = timeSlice.second;
749 
750  for (auto const & energyDeposit : energyDeposits)
751  {
752  double energy = energyDeposit.energy; // energy not attenuated
753  int trackID = energyDeposit.trackID;
754 
755  if (trackID < 0) { } // EM activity
756  else if (trackID > 0)
757  {
758  auto search = fParticleMap.find(trackID);
759 
760  if (search != fParticleMap.end())
761  {
762  const simb::MCParticle & particle = *((*search).second);
763  int pdg = particle.PdgCode(); // not EM activity so read what PDG it is
764 
765  if ((pdg == 11) || (pdg == -11) || (pdg == 22)) {}
766  else
767  {
768  enHAD += energy;
769  }
770  }
771  else { mf::LogWarning("TrainingDataAlg") << "PARTICLE NOT FOUND"; }
772  }
773  }
774  }
775  }
776  }
777 
778  return enHAD;
779 }
780 
781 /******************************************************************************/
782 
783 double proto::HadCal::GetEdepHADAtt_MC(art::Event const & e) const
784 {
785  double enHAD = 0.0;
786 
787  auto simchannelHandle = e.getHandle< std::vector<sim::SimChannel> >(fSimulationLabel);
788  if (simchannelHandle)
789  {
790 
791  for (auto const & channel : (*simchannelHandle) )
792  {
793  if (fGeometry->View(channel.Channel()) != fBestView) continue;
794 
795  // for every time slice in this channel:
796  auto const& timeSlices = channel.TDCIDEMap();
797  for (auto const& timeSlice : timeSlices)
798  {
799  // loop over the energy deposits.
800  auto const & energyDeposits = timeSlice.second;
801 
802  for (auto const & energyDeposit : energyDeposits)
803  {
804  double energy = energyDeposit.numElectrons * fElectronsToGeV * 1000; // attenuated
805  int trackID = energyDeposit.trackID;
806 
807  if (trackID < 0) { } // EM activity
808  else if (trackID > 0)
809  {
810  auto search = fParticleMap.find(trackID);
811 
812  if (search != fParticleMap.end())
813  {
814  const simb::MCParticle & particle = *((*search).second);
815  int pdg = particle.PdgCode(); // not EM activity so read what PDG it is
816 
817  if ((pdg == 11) || (pdg == -11) || (pdg == 22)) {}
818  else
819  {
820  enHAD += energy;
821  }
822  }
823  else { mf::LogWarning("TrainingDataAlg") << "PARTICLE NOT FOUND"; }
824  }
825  }
826  }
827  }
828  }
829 
830  return enHAD;
831 }
832 
833 /******************************************************************************/
834 
835 double proto::HadCal::GetEdepHADhAtt_MC(art::Event const & e,
836  detinfo::DetectorClocksData const& clockData,
837  detinfo::DetectorPropertiesData const& detProp,
838  const std::vector< recob::Hit > & hits) const
839 {
840  double tothadhit = 0.0;
841 
842  auto simChannelHandle = e.getValidHandle< std::vector<sim::SimChannel> >(fSimulationLabel);
843  const std::vector< sim::SimChannel > & channels = *simChannelHandle;
844 
845  for (auto const & hit: hits)
846  {
847  // the channel associated with this hit.
848  auto hitChannelNumber = hit.Channel();
849 
850  double hitEn = 0.0; double hitEnTrk = 0;
851 
852  for (auto const & channel : channels)
853  {
854  if (channel.Channel() != hitChannelNumber) continue;
855  if (fGeometry->View(channel.Channel()) != fBestView) continue;
856 
857  // for every time slice in this channel:
858  auto const& timeSlices = channel.TDCIDEMap();
859  for (auto const& timeSlice : timeSlices)
860  {
861  int time = timeSlice.first;
862  if (std::abs(hit.TimeDistanceAsRMS(time)) < 1.0)
863  {
864  // loop over the energy deposits.
865  auto const & energyDeposits = timeSlice.second;
866 
867  for (auto const & energyDeposit : energyDeposits)
868  {
869  int trackID = energyDeposit.trackID;
870 
871  double energy = energyDeposit.numElectrons * fElectronsToGeV * 1000; // attenuated
872  hitEn += energy;
873 
874  if (trackID < 0) { } // EM activity
875  else if (trackID > 0)
876  {
877  auto search = fParticleMap.find(trackID);
878  if (search != fParticleMap.end())
879  {
880  const simb::MCParticle & particle = *((*search).second);
881  int pdg = particle.PdgCode(); // not EM activity so read what PDG it is
882 
883  if ((pdg == 11) || (pdg == -11) || (pdg == 22)) {}
884  else {hitEnTrk += energy;}
885  }
886  else { mf::LogWarning("TrainingDataAlg") << "PARTICLE NOT FOUND"; }
887  }
888  }
889  }
890  }
891  }
892 
893  double ratio = 0.0;
894  if (hitEn > 0)
895  {
896  ratio = hitEnTrk / hitEn;
897  }
898  // use energy of hit corrected for electron lifetime
899  if (ratio > 0.5)
900  {
901  tothadhit += GetEhitMeVcorrel(clockData, detProp, hit);
902  }
903  }
904 
905  return tothadhit;
906 }
907 
908 /******************************************************************************/
909 
910 void proto::HadCal::ResetVars()
911 {
912  fRun = 0;
913  fEvent = 0;
914  fBestView = 2;
915  fNumberOfTracks = 0;
916 
917  fEdepEM_MC = 0.0;
918  fEdepEMAtt_MC = 0.0;
919  fEdepEMhAtt_MC = 0.0;
920 
921  fEdepHAD_MC = 0.0;
922  fEdepHADAtt_MC = 0.0;
923  fEdepHADhAtt_MC = 0.0;
924 
925  fMissEn_MC = 0.0;
926 
927  fEdepEM_reco = 0.0;
928  fEdepEMcorr_ellifetime_reco = 0.0;
929  fEdepHAD_reco = 0.0;
930  fEdepHADcorr_ellifetime_reco = 0.0;
931 
932  fEnGen = 0.0;
933  fEkGen = 0.0;
934 
935  fEdepSum = 0.0;
936  fEdepAllhits = 0.0;
937 
938  fdQdx = 0.0;
939  fdEdx = 0.0;
940  fdQ = 0.0;
941  fdx = 0.0;
942  fResRange = 0.0;
943  fT0 = 0.0;
944  fParticleMap.clear();
945 }
946 
947 DEFINE_ART_MODULE(proto::HadCal)
t0
code to link reconstructed objects back to the MC truth information
Definition: DirectHitParticleAssns_tool.cc:12
LArG4Parameters.h
Store parameters for running LArG4.
proto::HadCal::beginJob
void beginJob() override
Definition: HadCal_module.cc:178
simb::MCParticle::PdgCode
int PdgCode() const
Definition: MCParticle.h:212
proto::HadCal::endJob
void endJob() override
Definition: HadCal_module.cc:222
generate_CCQE_events.energy
string energy
Definition: generate_CCQE_events.py:15
proto::HadCal::fGeometry
geo::GeometryCore const * fGeometry
Definition: HadCal_module.cc:129
proto::HadCal::HadCal
HadCal(fhicl::ParameterSet const &p)
Definition: HadCal_module.cc:166
proto::HadCal::fHitModuleLabel
std::string fHitModuleLabel
Definition: HadCal_module.cc:155
LArVoxelList.h
Container of LAr voxel information.
art::DataViewImpl::getHandle
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
string
std::string string
Definition: nybbler.cc:12
proto::HadCal::file
std::ofstream file
Definition: HadCal_module.cc:151
recob::Hit::WireID
geo::WireID WireID() const
Definition: Hit.h:233
proto::HadCal::analyze
void analyze(art::Event const &e) override
Definition: HadCal_module.cc:237
proto::HadCal::GetEdepHitsMeV
double GetEdepHitsMeV(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const std::vector< recob::Hit > &hits) const
Definition: HadCal_module.cc:461
proto::HadCal::fTrackModuleLabel
std::string fTrackModuleLabel
Definition: HadCal_module.cc:157
vector
struct vector vector
TrackHitMeta.h
Class to keep data related to recob::Hit associated with recob::Track.
proto::HadCal::fEdepEMhAtt_MC
double fEdepEMhAtt_MC
Definition: HadCal_module.cc:104
proto::HadCal::fEdepHADhAtt_MC
double fEdepHADhAtt_MC
Definition: HadCal_module.cc:119
recob::Hit::Integral
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:224
channel
uint8_t channel
Definition: CRTFragment.hh:201
art::EDAnalyzer::EDAnalyzer
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
MCParticle.h
Particle class.
proto::HadCal::fSimulationLabel
std::string fSimulationLabel
Definition: HadCal_module.cc:159
proto::HadCal::reconfigure
void reconfigure(fhicl::ParameterSet const &p)
Definition: HadCal_module.cc:227
proto::HadCal::fEdepHAD_reco
double fEdepHAD_reco
Definition: HadCal_module.cc:121
Geometry.h
art framework interface to geometry description
proto::HadCal::fEdepHADcorr_ellifetime_reco
double fEdepHADcorr_ellifetime_reco
Definition: HadCal_module.cc:122
proto::HadCal::fCalorimetryModuleLabel
std::string fCalorimetryModuleLabel
Definition: HadCal_module.cc:158
proto::HadCal::fNNetModuleLabel
art::InputTag fNNetModuleLabel
Definition: HadCal_module.cc:154
calo::CalorimetryAlg::ElectronsFromADCArea
double ElectronsFromADCArea(double area, unsigned short plane) const
Definition: CalorimetryAlg.h:153
proto::HadCal::GetEhitMeVcorrel
double GetEhitMeVcorrel(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const recob::Hit &hit) const
Definition: HadCal_module.cc:493
reco_momentum_tuples.h
h
Definition: reco_momentum_tuples.py:66
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
proto::HadCal::fClusterModuleLabel
std::string fClusterModuleLabel
Definition: HadCal_module.cc:156
e
const double e
Definition: gUpMuFluxGen.cxx:165
proto::HadCal::operator=
HadCal & operator=(HadCal const &)=delete
LArVoxelData.h
Encapsulates the information we want store for a voxel.
DEFINE_ART_MODULE
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
reco_momentum_tuples.t
t
Definition: reco_momentum_tuples.py:25
fhicl::ParameterSet::get
T get(std::string const &key) const
Definition: ParameterSet.h:271
proto::HadCal::GetEdepHADAtt_MC
double GetEdepHADAtt_MC(art::Event const &e) const
Definition: HadCal_module.cc:783
art::DataViewImpl::getValidHandle
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:441
test.p
p
Definition: test.py:223
proto::HadCal::fEdepEMAtt_MC
double fEdepEMAtt_MC
Definition: HadCal_module.cc:103
proto
General LArSoft Utilities.
Definition: protoDUNEBeam.cxx:24
proto::HadCal::fElectronsToGeV
double fElectronsToGeV
Definition: HadCal_module.cc:145
search
Definition: search.py:1
art::DataViewImpl::run
RunNumber_t run() const
Definition: DataViewImpl.cc:71
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
geo::GeometryCore
Description of geometry of one entire detector.
Definition: GeometryCore.h:1468
geo::GeometryCore::View
View_t View(geo::PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
Definition: GeometryCore.cxx:601
proto::HadCal::fTreeEntries
TTree * fTreeEntries
Definition: HadCal_module.cc:132
hit
Detector simulation of raw signals on wires.
Definition: NumberOfHitsFilter_module.cc:25
make_THn_beam_input.data
data
Definition: make_THn_beam_input.py:80
proto::HadCal::GetEhitMeV
double GetEhitMeV(const recob::Hit &hit) const
Definition: HadCal_module.cc:513
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.
detinfo::DetectorClocksData
Contains all timing reference information for the detector.
Definition: DetectorClocksData.h:283
proto::HadCal::fCalorimetryAlg
calo::CalorimetryAlg fCalorimetryAlg
Definition: HadCal_module.cc:160
proto::HadCal::GetEdepEM_MC
double GetEdepEM_MC(art::Event const &e) const
Definition: HadCal_module.cc:531
proto::HadCal::fParticleMap
std::unordered_map< int, const simb::MCParticle * > fParticleMap
Definition: HadCal_module.cc:162
mf::LogWarning
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Definition: MessageLogger.h:210
Track.h
Provides recob::Track data product.
calo::CalorimetryAlg::dEdx_AREA
double dEdx_AREA(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
Definition: CalorimetryAlg.cxx:131
run_cnn_1class.pdg
pdg
Definition: run_cnn_1class.py:54
art::EventID::event
EventNumber_t event() const
Definition: EventID.h:116
proto::HadCal::GetEdepHADhAtt_MC
double GetEdepHADhAtt_MC(art::Event const &e, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const std::vector< recob::Hit > &hits) const
Definition: HadCal_module.cc:835
GeometryCore.h
Access the description of detector geometry.
proto::HadCal::fEdepHADAtt_MC
double fEdepHADAtt_MC
Definition: HadCal_module.cc:118
recob::Hit
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48
art::get
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
calo::CalorimetryAlg::LifetimeCorrection
double LifetimeCorrection(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, double time, double T0=0) const
Definition: CalorimetryAlg.cxx:246
PhysicalConstants.h
Collection of Physical constants used in LArSoft.
proto::HadCal::GetEdepEMhAtt_MC
double GetEdepEMhAtt_MC(art::Event const &e, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const std::vector< recob::Hit > &hits) const
Definition: HadCal_module.cc:647
proto::HadCal::GetEdepHAD_MC
double GetEdepHAD_MC(art::Event const &e) const
Definition: HadCal_module.cc:731
anab::MVAReader::create
static std::unique_ptr< MVAReader > create(const art::Event &evt, const art::InputTag &tag)
Definition: MVAReader.h:110
art::Event::id
EventID id() const
Definition: Event.cc:34
sim::LArG4Parameters::GeVToElectrons
double GeVToElectrons() const
Definition: LArG4Parameters.h:54
proto::HadCal::GetEdepEMAtt_MC
double GetEdepEMAtt_MC(art::Event const &e) const
Definition: HadCal_module.cc:589
proto::HadCal::GetEkinMeV
double GetEkinMeV(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const std::vector< art::Ptr< recob::Hit > > &hits, const std::vector< recob::TrackHitMeta const * > &data)
Definition: HadCal_module.cc:418
proto::HadCal::fEdepEMcorr_ellifetime_reco
double fEdepEMcorr_ellifetime_reco
Definition: HadCal_module.cc:107
Generated by   doxygen 1.8.11