EdepCal_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: EdepCal
3 // Module Type: analyzer
4 // File: EdepCal_module.cc
5 //
6 //
7 //If you have some questions you can write to dorota.stefan@cern.ch
8 //or robert.sulej@cern.ch
9 //
10 ////////////////////////////////////////////////////////////////////////
11 
28 
31 #include "canvas/Persistency/Common/FindManyP.h"
37 #include "art_root_io/TFileService.h"
39 #include "fhiclcpp/ParameterSet.h"
41 
42 #include "TH1.h"
43 #include "TTree.h"
44 #include "TLorentzVector.h"
45 #include "TVector3.h"
46 
47 #include <cmath>
48 
49 namespace proto
50 {
51  struct bHitInfo;
52  class EdepCal;
53 }
54 
55 struct proto::bHitInfo
56 {
57  bHitInfo(size_t i, double x, double e, int w) :
58  Index(i), dE(e), dx(x), wire(w)
59  { }
60  size_t Index;
61  double dE, dx;
62  int wire;
63 };
64 
66 public:
67  explicit EdepCal(fhicl::ParameterSet const & p);
68  // The destructor generated by the compiler is fine for classes
69  // without bare pointers or other resource use.
70 
71  // Plugins should not be copied or assigned.
72  EdepCal(EdepCal const &) = delete;
73  EdepCal(EdepCal &&) = delete;
74  EdepCal & operator = (EdepCal const &) = delete;
75  EdepCal & operator = (EdepCal &&) = delete;
76 
77  // Required functions.
78  void analyze(art::Event const & e) override;
79 
80  void beginJob() override;
81 
82  void beginRun(const art::Run& run) override;
83 
84  void reconfigure(fhicl::ParameterSet const& p) ;
85 private:
86 
87  // Declare member data here.
88  void ResetVars();
89 
92 
93  bool Has(std::vector<int> v, int i) const;
94 
95  double GetEdepMC(art::Event const & e) const;
96 
97  double GetEdepAttenuatedMC(art::Event const & e) const;
98 
99  double GetEdepEM_MC(art::Event const & e) const;
100 
101  double GetEdepEMAttenuated_MC(art::Event const & e) const;
102 
103 // double GetEdepTotVox(art::Event const & e) const;
104 
105  double GetEdepHits(detinfo::DetectorClocksData const& clockData,
106  detinfo::DetectorPropertiesData const& detProp,
107  const std::vector< recob::Hit > & hits) const;
108 
109  double GetEdepHits(detinfo::DetectorClocksData const& clockData,
110  detinfo::DetectorPropertiesData const& detProp,
111  const std::vector< art::Ptr<recob::Hit> > & hits) const;
112 
113  double GetEdepHitsMeV(const std::vector< recob::Hit > & hits) const;
114 
116 
117  double fT0;
118 
119  //////
120  TTree *fTree;
121  int fRun;
122  int fEvent;
123  double fEnGen;
124  double fEkGen;
125  double fEdep;
126  double fEdepMeV;
127  double fEdepCl;
128  double fEdepMC;
129  double fEdepAttMC;
130 // double fEdepMCTotV;
131  double fEdepEMMC;
132  double fEdepEMAttMC;
133  double fRatioTot;
134  double fRatioEM;
135  double fRatioHad;
136  //////
137 
138  //std::vector< art::Ptr<simb::MCParticle> > fSimlist;
139  //std::vector< art::Ptr<recob::Hit> > fHitlist;
140  //std::vector< art::Ptr<recob::Cluster> > fClusterlist;
141 
142  // Module labels to get data products
146 
148 
149  std::unordered_map< int, const simb::MCParticle* > fParticleMap;
150 };
151 
153  :
154  EDAnalyzer(p),
155  fCalorimetryAlg(p.get<fhicl::ParameterSet>("CalorimetryAlg"))
156  // More initializers here.
157 {
158  // get a pointer to the geometry service provider
160 
161  reconfigure(p);
162 }
163 
165 {
167  fElectronsToGeV = 1./larParameters->GeVToElectrons();
168 }
169 
171 {
172  // access art's TFileService, which will handle creating and writing hists
174 
175  fTree = tfs->make<TTree>("calibration","calibration tree");
176  fTree->Branch("fRun", &fRun, "fRun/I");
177  fTree->Branch("fEvent", &fEvent, "fEvent/I");
178  fTree->Branch("fEnGen", &fEnGen, "fEnGen/D");
179  fTree->Branch("fEkGen", &fEkGen, "fEkGen/D");
180  fTree->Branch("fEdep", &fEdep, "fEdep/D");
181  fTree->Branch("fEdepMeV", &fEdepMeV, "fEdepMeV/D");
182  fTree->Branch("fEdepCl", &fEdepCl, "fEdepCl/D");
183  fTree->Branch("fEdepMC", &fEdepMC, "fEdepMC/D");
184  fTree->Branch("fEdepAttMC", &fEdepAttMC, "fEdepAttMC/D");
185 // fTree->Branch("fEdepMCTotV", &fEdepMCTotV, "fEdepMCTotV/D");
186  fTree->Branch("fEdepEMMC", &fEdepEMMC, "fEdepEMMC/D");
187  fTree->Branch("fEdepEMAttMC", &fEdepEMAttMC, "fEdepEMAttMC/D");
188  fTree->Branch("fRatioTot", &fRatioTot, "fRatioTot/D");
189  fTree->Branch("fRatioEM", &fRatioEM, "fRatioEM/D");
190  fTree->Branch("fRatioHad", &fRatioHad, "fRatioHad/D");
191 }
192 
194 {
195  fSimulationLabel = p.get< std::string >("SimulationLabel");
196  fHitsModuleLabel = p.get< std::string >("HitsModuleLabel");
197  fClusterModuleLabel = p.get< std::string >("ClusterModuleLabel");
198 
199  fBestview = p.get<int>("Bestview");
200 }
201 
203 {
204  // Implementation of required member function here.
205  ResetVars();
206 
207  fRun = e.run();
208  fEvent = e.id().event();
209 
210  // MC particle list
211  auto particleHandle = e.getValidHandle< std::vector<simb::MCParticle> >(fSimulationLabel);
212  bool flag = true;
213 
214  for (auto const& p : *particleHandle)
215  {
216  fParticleMap[p.TrackId()] = &p;
217  if ((p.Process() == "primary") && flag)
218  {
219  fEnGen = p.P();
220  fEkGen = (std::sqrt(p.P()*p.P() + p.Mass()*p.Mass()) - p.Mass()) * 1000; // MeV
221  fT0 = p.T();
222  flag = false;
223  }
224  }
225 
226  // MC
227  fEdepMC = GetEdepMC(e);
229 // fEdepMCTotV = GetEdepTotVox(e);
230  fEdepEMMC = GetEdepEM_MC(e);
232 
233  //sim::IDE
234  //auto simchannel = e.getValidHandle< std::vector<sim::SimChannel> >(fSimulationLabel);
235 
236  // hits
237  fEdep = 0.0;
238  const auto& hitListHandle = *e.getValidHandle< std::vector<recob::Hit> >(fHitsModuleLabel);
239 
240  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService>()->DataFor(e);
241  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService>()->DataFor(e, clockData);
242  fEdep = GetEdepHits(clockData, detProp, hitListHandle);
243  fEdepMeV = GetEdepHitsMeV(hitListHandle);
244 
245  // clusters
246  auto clListHandle = e.getValidHandle< std::vector<recob::Cluster> >(fClusterModuleLabel);
247 
248  // look for all hits associated to selected cluster
249  art::FindManyP< recob::Hit > hitsFromClusters(clListHandle, e, fClusterModuleLabel);
250 
251  fEdepCl = 0.0;
252  for (size_t c = 0; c < clListHandle->size(); ++c)
253  {
254  fEdepCl += GetEdepHits(clockData, detProp, hitsFromClusters.at(c));
255  }
256 
257 
258  if (fEdepMC > 0.0)
259  {
260  fRatioTot = fEdep / fEdepMC;
261  }
262  if (fEdepEMMC > 0.0)
263  {
265  }
266 
267  double edephad = fEdepMC - fEdepEMMC;
268  if (edephad > 0)
269  {
270  fRatioHad = (fEdep - fEdepCl) / edephad;
271  }
272 
273  fTree->Fill();
274 }
275 
277 {
278  double energy = 0.0;
279 
280  auto simchannelHandle = e.getHandle< std::vector<sim::SimChannel> >(fSimulationLabel);
281  if (simchannelHandle)
282  {
283  for ( auto const& channel : (*simchannelHandle) )
284  {
285  if (fGeometry->View(channel.Channel()) == fBestview)
286  {
287  // for every time slice in this channel:
288  auto const& timeSlices = channel.TDCIDEMap();
289  for ( auto const& timeSlice : timeSlices )
290  {
291  // loop over the energy deposits.
292  auto const& energyDeposits = timeSlice.second;
293 
294  for ( auto const& energyDeposit : energyDeposits )
295  {
296  energy += energyDeposit.energy;
297  }
298  }
299  }
300  }
301  }
302 
303  return energy;
304 }
305 
307 {
308  double energy = 0.0;
309 
310  auto simchannelHandle = e.getHandle< std::vector<sim::SimChannel> >(fSimulationLabel);
311  if (simchannelHandle)
312  {
313  for ( auto const& channel : (*simchannelHandle) )
314  {
315  if (fGeometry->View(channel.Channel()) == fBestview)
316  {
317  // for every time slice in this channel:
318  auto const& timeSlices = channel.TDCIDEMap();
319  for ( auto const& timeSlice : timeSlices )
320  {
321  // loop over the energy deposits.
322  auto const& energyDeposits = timeSlice.second;
323 
324  for ( auto const& energyDeposit : energyDeposits )
325  {
326  energy += energyDeposit.numElectrons * fElectronsToGeV * 1000;
327  }
328  }
329  }
330  }
331  }
332 
333  return energy;
334 
335 }
336 
338 {
339  double enEM = 0.0;
340 
341  auto simchannelHandle = e.getHandle< std::vector<sim::SimChannel> >(fSimulationLabel);
342  if (simchannelHandle)
343  {
344  for ( auto const& channel : (*simchannelHandle) )
345  {
346  if (fGeometry->View(channel.Channel()) == fBestview)
347  {
348  // for every time slice in this channel:
349  auto const& timeSlices = channel.TDCIDEMap();
350  for ( auto const& timeSlice : timeSlices )
351  {
352  // loop over the energy deposits.
353  auto const& energyDeposits = timeSlice.second;
354 
355  for ( auto const& energyDeposit : energyDeposits )
356  {
357  double energy = energyDeposit.energy;
358  int trackID = energyDeposit.trackID;
359 
360  if (trackID < 0)
361  {
362  enEM += energy;
363  }
364  else if (trackID > 0)
365  {
366  auto search = fParticleMap.find(trackID);
367  bool found = true;
368  if (search == fParticleMap.end())
369  {
370  mf::LogWarning("TrainingDataAlg") << "PARTICLE NOT FOUND";
371  found = false;
372  }
373 
374  int pdg = 0;
375  if (found)
376  {
377  const simb::MCParticle& particle = *((*search).second);
378  if (!pdg) pdg = particle.PdgCode(); // not EM activity so read what PDG it is
379  }
380 
381  if ((pdg == 11) || (pdg == -11) || (pdg == 22)) enEM += energy;
382  }
383 
384  }
385  }
386  }
387  }
388  }
389  return enEM;
390 }
391 
393 {
394  double enEM = 0.0;
395 
396 
397  auto simchannelHandle = e.getHandle< std::vector<sim::SimChannel> >(fSimulationLabel);
398  if (simchannelHandle)
399  {
400  for ( auto const& channel : (*simchannelHandle) )
401  {
402  if (fGeometry->View(channel.Channel()) == fBestview)
403  {
404  // for every time slice in this channel:
405  auto const& timeSlices = channel.TDCIDEMap();
406  for ( auto const& timeSlice : timeSlices )
407  {
408  // loop over the energy deposits.
409  auto const& energyDeposits = timeSlice.second;
410 
411  for ( auto const& energyDeposit : energyDeposits )
412  {
413  double energy = energyDeposit.numElectrons * fElectronsToGeV * 1000;
414  int trackID = energyDeposit.trackID;
415 
416  if (trackID < 0)
417  {
418  enEM += energy;
419  }
420  else if (trackID > 0)
421  {
422  auto search = fParticleMap.find(trackID);
423  bool found = true;
424  if (search == fParticleMap.end())
425  {
426  mf::LogWarning("TrainingDataAlg") << "PARTICLE NOT FOUND";
427  found = false;
428  }
429 
430  int pdg = 0;
431  if (found)
432  {
433  const simb::MCParticle& particle = *((*search).second);
434  if (!pdg) pdg = particle.PdgCode(); // not EM activity so read what PDG it is
435  }
436 
437  if ((pdg == 11) || (pdg == -11) || (pdg == 22)) enEM += energy;
438  }
439 
440  }
441  }
442  }
443  }
444  }
445  return enEM;
446 }
447 // before recombination takes place and other attenuations
448 /*double proto::EdepCal::GetEdepTotVox(art::Event const & e) const
449 {
450  double en = 0.0;
451 
452  const sim::LArVoxelList voxels = sim::SimListUtils::GetLArVoxelList(e, fSimulationLabel);
453 
454  sim::LArVoxelList::const_iterator vxitr;
455  for (vxitr = voxels.begin(); vxitr != voxels.end(); vxitr++)
456  {
457  const sim::LArVoxelData &vxd = (*vxitr).second;
458  for (size_t partIdx = 0; partIdx < vxd.NumberParticles(); partIdx++)
459  {
460  en += vxd.Energy(partIdx);
461  }
462  }
463 
464  return en;
465 }*/
466 
468  detinfo::DetectorPropertiesData const& detProp,
469  const std::vector< recob::Hit > & hits) const
470 {
471  if (!hits.size()) return 0.0;
472 
473  double dqsum = 0.0;
474  for (size_t h = 0; h < hits.size(); ++h)
475  {
476  unsigned short plane = hits[h].WireID().Plane;
477  if (plane != fBestview) continue;
478 
479  double dqadc = hits[h].Integral();
480  if (!std::isnormal(dqadc) || (dqadc < 0)) continue;
481 
482  double dqel = fCalorimetryAlg.ElectronsFromADCArea(dqadc, plane);
483 
484  double tdrift = hits[h].PeakTime();
485  double correllifetime = fCalorimetryAlg.LifetimeCorrection(clockData, detProp, tdrift, fT0);
486 
487  double dq = dqel * correllifetime * fElectronsToGeV * 1000;
488  if (!std::isnormal(dq) || (dq < 0)) continue;
489 
490  dqsum += dq;
491  }
492 
493  return dqsum;
494 }
495 
496 double proto::EdepCal::GetEdepHitsMeV(const std::vector< recob::Hit > & hits) const
497 {
498  if (!hits.size()) return 0.0;
499 
500  double dqsum = 0.0;
501  for (size_t h = 0; h < hits.size(); ++h)
502  {
503  unsigned short plane = hits[h].WireID().Plane;
504  if (plane != fBestview) continue;
505 
506  double dqadc = hits[h].Integral();
507  if (!std::isnormal(dqadc) || (dqadc < 0)) continue;
508 
509  double dqel = fCalorimetryAlg.ElectronsFromADCArea(dqadc, plane);
510 
511  double dq = dqel * fElectronsToGeV * 1000;
512  if (!std::isnormal(dq) || (dq < 0)) continue;
513 
514  dqsum += dq;
515  }
516 
517  return dqsum;
518 }
519 
521  detinfo::DetectorPropertiesData const& detProp,
522  const std::vector< art::Ptr<recob::Hit> > & hits) const
523 {
524  if (!hits.size()) return 0.0;
525 
526  double dqsum = 0.0;
527  for (size_t h = 0; h < hits.size(); ++h)
528  {
529  unsigned short plane = hits[h]->WireID().Plane;
530  if (plane != fBestview) continue;
531 
532  double dqadc = hits[h]->Integral();
533  if (!std::isnormal(dqadc) || (dqadc < 0)) continue;
534 
535  double dqel = fCalorimetryAlg.ElectronsFromADCArea(dqadc, plane);
536 
537  double tdrift = hits[h]->PeakTime();
538  double correllifetime = fCalorimetryAlg.LifetimeCorrection(clockData, detProp, tdrift, fT0);
539 
540  double dq = dqel * correllifetime * fElectronsToGeV * 1000;
541  if (!std::isnormal(dq) || (dq < 0)) continue;
542 
543  dqsum += dq;
544  }
545 
546  return dqsum;
547 }
548 
549 bool proto::EdepCal::Has(std::vector<int> v, int i) const
550 {
551  for (auto c : v) if (c == i) return true;
552  return false;
553 }
554 
555 
557 {
558  fEdep = 0.0;
559  fEdepMeV = 0.0;
560  fEnGen = 0.0;
561  fEkGen = 0.0;
562  fEdepCl = 0.0;
563  fEdepMC = 0.0;
564  fEdepAttMC = 0.0;
565 // fEdepMCTotV = 0.0;
566  fEdepEMMC = 0.0;
567  fEdepEMAttMC = 0.0;
568  fRatioTot = 0.0;
569  fRatioEM = 0.0;
570  fRatioHad = 0.0;
571  fT0 = 0.0;
572  fParticleMap.clear();
573 }
574 
double GetEdepAttenuatedMC(art::Event const &e) const
def analyze(root, level, gtrees, gbranches, doprint)
Definition: rootstat.py:69
double GetEdepMC(art::Event const &e) const
Store parameters for running LArG4.
void beginJob() override
int PdgCode() const
Definition: MCParticle.h:212
std::string fSimulationLabel
double fElectronsToGeV
Container of LAr voxel information.
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
std::string string
Definition: nybbler.cc:12
void analyze(art::Event const &e) override
struct vector vector
void reconfigure(fhicl::ParameterSet const &p)
uint8_t channel
Definition: CRTFragment.hh:201
unsigned int Index
Particle class.
art framework interface to geometry description
Definition: Run.h:17
double GetEdepHits(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const std::vector< recob::Hit > &hits) const
void beginRun(const art::Run &run) override
double ElectronsFromADCArea(double area, unsigned short plane) const
double GetEdepEMAttenuated_MC(art::Event const &e) const
const double e
Encapsulates the information we want store for a voxel.
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
void beginJob()
Definition: Breakpoints.cc:14
virtual void reconfigure(fhicl::ParameterSet const &pset)
T get(std::string const &key) const
Definition: ParameterSet.h:271
geo::GeometryCore const * fGeometry
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:441
p
Definition: test.py:223
General LArSoft Utilities.
Definition: search.py:1
RunNumber_t run() const
Definition: DataViewImpl.cc:71
Description of geometry of one entire detector.
bHitInfo(size_t i, double x, double e, int w)
View_t View(geo::PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
EdepCal(fhicl::ParameterSet const &p)
Declaration of signal hit object.
Contains all timing reference information for the detector.
calo::CalorimetryAlg fCalorimetryAlg
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
EventNumber_t event() const
Definition: EventID.h:116
Access the description of detector geometry.
list x
Definition: train.py:276
std::string fClusterModuleLabel
double GetEdepHitsMeV(const std::vector< recob::Hit > &hits) const
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
double LifetimeCorrection(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, double time, double T0=0) const
Collection of Physical constants used in LArSoft.
std::unordered_map< int, const simb::MCParticle * > fParticleMap
double GetEdepEM_MC(art::Event const &e) const
std::string fHitsModuleLabel
bool Has(std::vector< int > v, int i) const
EventID id() const
Definition: Event.cc:34
double GeVToElectrons() const