EnergyStandardCandles_module.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <memory>
3 #include <numeric>
4 #include <random>
5 #include <stdexcept>
6 
7 #include "TGeoMaterial.h"
8 #include "TParameter.h"
9 #include "TROOT.h"
10 #include "TTree.h"
11 #include "TVector3.h"
12 
20 
21 #define ML_DEBUG
23 
24 #include "fhiclcpp/ParameterSet.h"
25 
26 #include "Geometry/Geometry.h"
27 #include "Geometry/LiveGeometry.h"
28 #include "Utilities/AssociationUtil.h"
29 
30 #include "MCCheater/BackTracker.h"
31 #include "MEFinder/MEClusters.h"
32 #include "RecoBase/CellHit.h"
33 #include "RecoBase/RecoHit.h"
34 #include "RecoBase/Track.h"
35 #include "SummaryData/SpillData.h"
36 
37 
38 
39 #include "StdCandlesNtupleVars.h"
40 
41 namespace calib
42 {
43  static const double GeV2MeV = CLHEP::GeV / CLHEP::MeV;
44 
46  {
47  public:
48  explicit EnergyStandardCandles(const fhicl::ParameterSet& pset);
49  virtual ~EnergyStandardCandles() {};
50 
51  void analyze(const art::Event& evt);
52 
53  void reconfigure(const fhicl::ParameterSet& pset);
54 
55  void beginJob();
56  void endJob();
57 
58  private:
59  /// Utility type to hold info about stopping muons
60  struct StoppingMu
61  {
63  const std::vector<art::Ptr<me::TrkME> > & cs,
64  const art::PtrVector<rb::CellHit> & dblCnt,
65  const std::vector<cheat::TrackIDE> & truth)
66  : track(t), michelEClusters(cs), doubleCountedHits(dblCnt), michelTruth(truth)
67  {};
68  const rb::Track & track;
69  const std::vector<art::Ptr<me::TrkME> > michelEClusters;
71  const std::vector<cheat::TrackIDE> michelTruth;
72  };
73 
74  // methods to find stopping muons
75  void FindStoppingMuons(const art::Event& evt);
76  bool IsContained(const rb::Track & trk) const;
77  bool IsMuon(const rb::Track & trk) const;
78 
79  // methods to find pi0s
80  void FindPi0s(const art::Event& evt);
81 
82  // methods to make ntuples
83  void PrepareMichelENtuple();
84  void PrepareMuondEdXNtuple();
85  void PreparePi0Ntuple();
86 
87  void FillMichelENtuple();
88  void FillMichelRecoInfo(const std::shared_ptr<StoppingMu> stoppingMu, MichelCandidate & michelCand);
89  void FillTruthInfo(const std::shared_ptr<StoppingMu> stoppingMu, MichelCandidate & michelCand);
90  void FillTrueMichelInfo(const sim::Particle * trueParticle, TrueMichel * eContrib);
91  void FixPhotonAncestry(const sim::Particle *& trueParticle, const sim::Particle * parentParticle, TrueEContrib * eContrib);
92 
93  void FillMuondEdXNtuple();
94  void FillPi0Ntuple();
95 
96  std::vector<std::shared_ptr<StoppingMu> > fStoppingMus;
97 
98  struct Pi0
99  {
100  };
101  std::vector<std::shared_ptr<Pi0> > fPi0s;
102 
103  double fFrontSize; ///< maximum distance to front of detector to qualify as "front-entering"
104  std::string fGeneratorLabel; ///< name of generator information label in event record
105  std::string fHitLabel; ///< name of hits container in event record
106  std::string fMichelLabel; ///< name of Michel electron container in event record
107  std::string fNumiBeamLabel; ///< name of NuMI spill information container in event record
108  double fMinDistToEdge; ///< minimum distance in cm to edge of detector to qualify as "contained"
109  std::string fScintillatorName; ///< name of geometry TGeoMaterial corresponding to scintillator
110  std::string fTrackLabel; ///< name of tracks container in event record
111 
112  std::map<std::string, int> fHeaderInfo; ///< run/subrun/gate for output ntuple
113  TParameter<double> * fPOTStorage; ///< container for POT that will go to output ntuple
114  double fPOT;
115  MichelCandidate * fMichelCandidate; ///< used for filling output ntuple
116 
117  std::map<std::string, TTree*> fOutTrees;
118 
121  art::ServiceHandle<geo::LiveGeometry> fLiveGeom; // i.e., *instrumented* geometry (which instrumented cells are alive?)
123  };
124 
126  : EDAnalyzer(pset)
127  {
128  reconfigure(pset);
129  }
130 
132  {
133  this->fHeaderInfo = { {"run", 0}, {"subRun", 0}, {"spillNum", 0} };
134  this->fPOTStorage = this->fFileService->make<TParameter<double> >("POT", this->fPOT);
135  gDirectory->Append(this->fPOTStorage); // doesn't make it to the file without this. probably due to ROOT's crazy 'directory' system
136  this->fPOT = 0;
137 
138  this->PrepareMichelENtuple();
139  this->PrepareMuondEdXNtuple();
140  this->PreparePi0Ntuple();
141  }
142 
144  {
145  // store the POT in the file
146  this->fPOTStorage->SetVal(this->fPOT);
147 
148  if (fMichelCandidate)
149  delete fMichelCandidate;
150  } // EnergyStandardCandles::endJob()
151 
153  {
154  fGeneratorLabel = pset.get<std::string> ( "GeneratorInput" );
155  fHitLabel = pset.get<std::string> ( "HitLabel" );
156  fMichelLabel = pset.get<std::string> ( "MichelLabel" );
157  fMinDistToEdge = pset.get<double> ( "MinDistToEdge", 50. ); // TODO: what are the units here?
158  fNumiBeamLabel = pset.get<std::string> ( "NuMIBeamInput" );
159  fScintillatorName = pset.get<std::string> ( "ScintillatorName" );
160  fTrackLabel = pset.get<std::string> ( "TrackLabel" );
161  }
162 
164  {
165 
166  fHeaderInfo["run"] = evt.run();
167  fHeaderInfo["subRun"] = evt.subRun();
168  fHeaderInfo["spillNum"] = evt.id().event();
169 
172  evt.getByLabel(fGeneratorLabel, spillPot);
173  else
174  evt.getByLabel(fNumiBeamLabel, spillPot);
175 
176  if (spillPot.isValid()) // won't be for PC
177  this->fPOT += spillPot->spillpot;
178 
179  fStoppingMus.clear();
180 
181  this->FindStoppingMuons(evt);
182  this->FindPi0s(evt);
183 
184  if (fStoppingMus.size() > 0)
185  {
186  this->FillMichelENtuple();
187  this->FillMuondEdXNtuple();
188  }
189 
190  if (fPi0s.size() > 0)
191  {
192  this->FillPi0Ntuple();
193  }
194  }
195 
197  {
198 
199  // look at all the tracks.
200  // do any of them meet the stopping muon criteria?:
201  // - pass muon PID
202  // - stops inside fiducial volume
203  // - matched Michel electron candidate
204  // if so, store the relevant info for plotting later on.
205 
207  evt.getByLabel(fTrackLabel, tracks);
208  art::FindManyP<me::TrkME> michelEs(tracks, evt, fMichelLabel);
209 
210  // will need these for checking whether Michel finding is vacuuming up
211  // already tracked hits
212  std::set<art::Ptr<rb::CellHit> > hitsOnTracks;
213  for(const auto & track : *tracks)
214  {
215  const auto & hits = track.AllCells();
216  hitsOnTracks.insert(hits.begin(), hits.end());
217  }
218 
219  mf::LogDebug("EnergyStandardCandles") << "There are " << tracks->size() << " tracks in this event.";
220 
221  auto it_track = tracks->begin();
222  for (unsigned long int trk_idx = 0; it_track != tracks->end(); trk_idx++, it_track++)
223  {
224  mf::LogDebug("EnergyStandardCandles") << "Considering track #" << trk_idx << ":";
225  const rb::Track & track = *it_track;
226  if ( !IsContained(track) )
227  {
228  mf::LogDebug("EnergyStandardCandles") << " track not contained. skip.";
229  continue;
230  }
231 
232  if ( !IsMuon(track) )
233  {
234  mf::LogDebug("EnergyStandardCandles") << " track not a muon. skip.";
235  continue;
236  }
237 
238  // because of the association it's kind of messy to package the Michel matching
239  // into a function (requires lots of arguments). it's short so just do it here
240  std::vector<art::Ptr<me::TrkME> > matchedMichelClusters = michelEs.at(trk_idx);
241  if ( matchedMichelClusters.size() == 0 )
242  {
243  mf::LogDebug("EnergyStandardCandles") << " track has no matched Michel clusters. skip.";
244  continue;
245  }
246 
247  // grab all the hits from all the Michel-matched clusters.
248  // we need to pass hits (not clusters) to the BackTracker?
249  art::PtrVector<rb::CellHit> michelHits;
250  for (auto michelClus : matchedMichelClusters)
251  {
252  auto clusHits = michelClus->AllCells();
253  michelHits.insert(michelHits.end(), clusHits.begin(), clusHits.end());
254 
255  }
256  // were any of these hits already owned by some other track?
257  // (need the set because set_intersection() assumes they're unique & sorted;
258  // shoving them into a std::set does that automatically)
259  std::set<art::Ptr<rb::CellHit> > michelHits_set(michelHits.begin(), michelHits.end());
261  std::set_intersection(michelHits_set.begin(),
262  michelHits_set.end(),
263  hitsOnTracks.begin(),
264  hitsOnTracks.end(),
265  ixn.begin());
266 
267  std::vector<cheat::TrackIDE> michelTruthInfo;
269  {
270  michelTruthInfo = fBackTracker->HitsToTrackIDE(michelHits);
271  for (const auto & geantTrackInfo : michelTruthInfo)
272  {
273  const sim::Particle * trueParticle = this->fBackTracker->TrackIDToParticle(geantTrackInfo.trackID);
274  mf::LogDebug("EnergyStandardCandles") << " frac from PID " << trueParticle->PdgCode() << " = " << geantTrackInfo.energyFrac;
275  mf::LogDebug("EnergyStandardCandles") << " (particle energy = " << trueParticle->E() << " MeV)";
276  }
277  }
278 
279  // this is (probably) a stopping muon! store info for histogramming later.
280  fStoppingMus.push_back( std::make_shared<StoppingMu>(track, std::move(matchedMichelClusters), std::move(ixn), michelTruthInfo) );
281  mf::LogDebug("EnergyStandardCandles") << " accepted.";
282  } // for (it_trk, trk_idx)
283  } // StandardCandleSelection::FindStoppingMuons()
284 
286  {
287  // cribbed from StopperSelection in this package
288 
289  const TVector3 v1 = trk.Start();
290  const TVector3 v2 = trk.Stop();
291  const TVector3 stop = (v1.Y() < v2.Y()) ? v1 : v2; // Start and end don't mean much, pick whichever is lowest in y.
292 
293  mf::LogDebug("EnergyStandardCandles") << "Track "
294  << "start pos: (" << v1.x() << "," << v1.y() << "," << v1.z() << "); "
295  << "stop pos: (" << v2.x() << "," << v2.y() << "," << v2.z() << ")";
296  if(fLiveGeom->DistToEdgeXY(stop) < fMinDistToEdge ||
298  return false;
299  if(v1.X() == v2.X() || v1.Y() == v2.Y())
300  return false; // dx or dy = 0 means bad fit or 2d track only
301 
302  return true;
303  }
304 
306  {
307  // just a stub for the moment
308  return true;
309  }
310 
312  {
313 
314  // look at the event.
315  // is there anything that meets the pi0 criteria?:
316  }
317 
319  {
320  this->fOutTrees["Michel"] = fFileService->make<TTree>("MichelE", "Michel electron candidates");
321 
322  for (auto & branchPair : this->fHeaderInfo)
323  {
324  std::string typeCode = branchPair.first + "/I";
325  this->fOutTrees["Michel"]->Branch(branchPair.first.c_str(), &branchPair.second, typeCode.c_str());
326  }
327 
328  this->fMichelCandidate = new MichelCandidate;
329  this->fOutTrees["Michel"]->Branch("michel_candidate", this->fMichelCandidate);
330  }
331 
332  void EnergyStandardCandles::FillMichelRecoInfo(const std::shared_ptr<StoppingMu> stoppingMu, MichelCandidate & michelCand)
333  {
334  double deltaT_Ewgtd = 0;
335  double deltaR_Ewgtd = 0;
336  double MID_Ewgtd = 0;
337  double wgt = 0; // going to do a pulse-height-weighted sum
338  double hit_wgt = 0;
339  TVector3 location;
340  michelCand.hitEnergies.clear(); // this buffer is re-used. gotta empty it
341 
342  double t_mean_running(0), t_M_2_running(0);
343  for (const auto & michelClus : stoppingMu->michelEClusters)
344  {
345  for (unsigned int idx_cell = 0; idx_cell < michelClus->NCell(); idx_cell++)
346  {
347  const auto & recoHit = michelClus->RecoHit(idx_cell);
348  if (!recoHit.IsCalibrated())
349  continue;
350  michelCand.hitEnergies.push_back( recoHit.GeV() * michelClus->Weight(idx_cell) * GeV2MeV );
351  }
352  double energy = michelClus->TotalGeV() * GeV2MeV;
353  deltaT_Ewgtd += energy * michelClus->DeltaT();
354  deltaR_Ewgtd += energy * michelClus->DistToTrk();
355  MID_Ewgtd += energy * michelClus->MID();
356  location += energy * michelClus->MeanXYZ();
357  wgt += energy;
358 
359  // single-pass calculation for variance from https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
360  for (const auto & hit : michelClus->WeightedHits())
361  {
362  double hit_pe = hit.hit->PE() * hit.weight;
363  double old_wgt = hit_wgt;
364  hit_wgt += hit_pe;
365  double delta = hit.hit->TNS() - t_mean_running;
366  double wgtd_delta = delta * hit_pe / hit_wgt;
367  t_mean_running += wgtd_delta;
368  t_M_2_running += old_wgt * delta * wgtd_delta;
369  mf::LogDebug("EnergyStandardCandles") << " calculating hit time variance: hit time, hit PE, hit_wgt, delta, t_mean_running, t_M_2_running = "
370  << hit.hit->TNS() << ", "
371  << hit_pe << ", "
372  << delta << ", "
373  << t_mean_running << ", "
374  << t_M_2_running;
375  }
376  }
377  if (wgt > 0)
378  {
379  michelCand.energy = wgt;
380  michelCand.deltaT = deltaT_Ewgtd / wgt;
381  michelCand.location = location * (1./wgt); // no operator/() in TVector3, strangely
382 // std::cout << " weighted location:" << std::endl;
383 // michelCand.location.Print();
384  michelCand.deltaR = deltaR_Ewgtd / wgt;
385  michelCand.MID = MID_Ewgtd / wgt;
386  }
387  else
388  {
389  michelCand.location = TVector3();
390  michelCand.deltaT = 0;
391  michelCand.deltaR = 0;
392  michelCand.MID = -999;
393  }
394 
395  if (hit_wgt > 0)
396  michelCand.t_RMS = sqrt(t_M_2_running/hit_wgt);
397  else
398  michelCand.t_RMS = 0;
399  mf::LogDebug("EnergyStandardCandles") << " hit time std. dev.: " << michelCand.t_RMS << std::endl;
400 
401  michelCand.nAlreadyUsedHits = stoppingMu->doubleCountedHits.size();
402  michelCand.alreadyUsedEnergy = std::accumulate(stoppingMu->doubleCountedHits.begin(),
403  stoppingMu->doubleCountedHits.end(),
404  0,
405  [] (double runningSum, const art::Ptr<rb::CellHit> & hit)
406  { return runningSum + hit->PE(); }
407  );
408 
409  michelCand.parentFrontZ = stoppingMu->track.Start().Z();
410 
411  } // EnergyStandardCandles::FillMichelRecoInfo
412 
413  /// Sometimes photon ancestries are more complicated than we need them to be
414  void EnergyStandardCandles::FixPhotonAncestry(const sim::Particle *& trueParticle, const sim::Particle * parentParticle, TrueEContrib * eContrib)
415  {
416  if (parentParticle && parentParticle->PdgCode() == 111)
417  trueParticle = parentParticle;
418  else if (parentParticle && trueParticle->Process() == "nCapture")
419  {
420  TrueNeutronCapture * nCapContrib = dynamic_cast<TrueNeutronCapture*>(eContrib);
421 
422  trueParticle = parentParticle;
423 // std::cout << " neutron capture. neutron daughter particles:" << std::endl;
424  for (int i = 0; i < parentParticle->NumberDaughters(); i++)
425  {
426  auto daughterParticle = this->fBackTracker->TrackIDToParticle(parentParticle->Daughter(i));
427  // in GEANT4 evidently the neutron capture byproduct nucleus
428  // is the only daughter of the captured neutron that has a nucleus
429  // PDG code (1000000000 to 10999999999) AND has the 'nCapture' process set.
430  // (the neutron might bump into other nuclei and knock them out via
431  // elastic or inelastic scattering, so we've got to look at the process too.)
432  // note that this tells you the RESULTING nucleus, not the STRUCK one.
433  // Hopefully just subtracting 1 from the atomic mass of the resulting nucleus
434  // yields the correct struck nucleus?
435  auto daughterPdg = daughterParticle->PdgCode();
436  if (daughterPdg > 1000000000 && daughterPdg < 1099999999 && daughterParticle->Process() == "nCapture")
437  {
438  // PDG code for nuclei is: 10LZZZAAAI
439  // (L = strangeness, ZZZ = 3-digit atomic number, AAA = 3-digit mass number, I = isomer number (resonant state))
440  nCapContrib->nCaptureNucleus = daughterPdg - 10;
441 // std::cout << " neutron captured on nucleus: " << nCapContrib->nCaptureNucleus << std::endl;
442  break;
443  }
444 // std::cout << " " << daughterParticle->PdgCode() << " (process: " << daughterParticle->Process() << ")" << std::endl;
445  } // for (i)
446  } // else if (is N capture)
447 
448  /*
449  else
450  {
451  std::cout << " photon with unclassified parent. Process: " << trueParticle->Process();
452  if (parentParticle)
453  std::cout << "; parent PDG: " << parentParticle->PdgCode();
454  std::cout << std::endl;
455  }
456 */
457 
458  } // EnergyStandardCandles::FixPhotonAncestry()
459 
461  {
462  if (!eContrib)
463  return;
464 
465  auto allHits = fBackTracker->ParticleToFLSHit(trueParticle->TrackId());
466  // auto onlyEhits = fBackTracker->ParticleToFLSHit(trueParticle->TrackId(), trueParticle->PdgCode()); // not trustworthy since there might be more particles in the daughter chain with the same PDG as the original e+/e-
467 
468 
469  mf::LogDebug("EnergyStandardCandles") << " True Michel electron:";
470 
471  double deadEnergy = 0; // that is, energy deposited in inactive material
472  double deadRange = 0;
473  double totalRange = 0;
474  const auto & traj = trueParticle->Trajectory();
475  {
476  // these braces to make the LogDebug object go out of scope
477  // after we're done with it and thus write to the stream
478  mf::LogDebug trajLog("EnergyStandardCandles");
479  trajLog << " Trajectory points:\n";
480  double lastE = -1;
481  TVector3 lastPos;
482  for (const auto & trajP : traj)
483  {
484  trajLog << " pos = (" << trajP.first.X() << "," << trajP.first.Y() << "," << trajP.first.Z() << ") mm;"
485  << " fourP = (" << trajP.second.E() << "," << trajP.second.Px() << "," << trajP.second.Py() << "," << trajP.second.Pz() << ") GeV\n";
486 
487  bool haveLast = lastE > 0;
488 
489  // compute how much of the energy in this trajectory step
490  // was lost in inactive ("dead") materials.
491  // can't do this exactly, because not enough information was
492  // retained here from GEANT. to first approximation, though,
493  // ionization energy (which is most of what we're worried about
494  // here, since the rad length is reasonably long) is proportional
495  // to density. so we weight the dE of this step
496  // by the fraction of the range (= density * step size)
497  // that was in "dead" materials. that's our "dead" energy.
498  if (haveLast)
499  {
500  double dE = lastE - trajP.second.E();
501  // now count how much of the energy was deposited in inactive material
502  std::vector<double> pathLengths;
503  std::vector<const TGeoMaterial*> materials;
504  fGeom->MaterialsBetweenPoints(lastPos, trajP.first.Vect(), pathLengths, materials);
505 
506  trajLog << " in this step:\n";
507  double stepTotalRange = 0;
508  double stepDeadRange = 0;
509  for (unsigned int idx_mat = 0; idx_mat < pathLengths.size(); idx_mat++)
510  {
511  const TGeoMaterial * material = materials[idx_mat];
512  double pathLength = pathLengths[idx_mat];
513  trajLog << " material: " << material->GetName()
514  << "; path length = " << pathLength << "mm;\n";
515  double stepRange = material->GetDensity() * (pathLength * CLHEP::mm / CLHEP::cm);
516  stepTotalRange += stepRange;
517  if (material->GetName() != fScintillatorName)
518  stepDeadRange += stepRange;
519  }
520 
521  if (stepTotalRange > 0)
522  deadEnergy += stepDeadRange / stepTotalRange * dE;
523 
524  deadRange += stepDeadRange;
525  totalRange += stepTotalRange;
526  } // if (haveLast)
527 
528  lastE = trajP.second.E();
529  lastPos = trajP.first.Vect();
530  } // for (traj)
531  } // scope for trajLog
532 
533  mf::LogDebug("EnergyStandardCandles") << " Children: " << trueParticle->NumberDaughters() << std::endl;
534 // double childEnergy = 0; // not sure what this quantity really tells us.
535  double bremEnergy = 0;
536  for (int i = 0; i < trueParticle->NumberDaughters(); i++)
537  {
538  auto daughterParticle = this->fBackTracker->TrackIDToParticle(trueParticle->Daughter(i));
539  double daughterE = daughterParticle->E();
540  TLorentzVector daughterPos = daughterParticle->Position();
541  mf::LogDebug("EnergyStandardCandles") << " PID: " << daughterParticle->PdgCode() << ";"
542  << " process name: " << daughterParticle->Process() << ";"
543  << " energy: " << daughterE * GeV2MeV << " MeV;"
544  << " creation point: (" << daughterPos.X() << "," << daughterPos.Y() << "," << daughterPos.Z() << ")";
545 // childEnergy += daughterE;
546  if (daughterParticle->Process() == "eBrem")
547  {
548  bremEnergy += daughterE;
549  // if the Brem. was created in 'dead' material,
550  // don't count the Brem's energy in the "dead energy" -- that double-counts
551  if (fGeom->Material(daughterPos.X(), daughterPos.Y(), daughterPos.Z())->GetName() != fScintillatorName )
552  deadEnergy -= daughterE;
553  }
554  // since the parent is an electron or positron,
555  // this is Moller/Bhaba scattering,
556  // which is delta ray ("knock-on electron") production by electrons.
557  // this will produce prompt, nearby energy that should be
558  // lumped together with the regular ionization from the Michel.
559  else if (daughterParticle->Process() == "eIoni")
560  {
561  auto daughterHits = fBackTracker->ParticleToFLSHit(daughterParticle->TrackId());
562  allHits.insert(allHits.end(), daughterHits.begin(), daughterHits.end());
563  }
564  }
565  eContrib->ElostToBrem = bremEnergy * GeV2MeV;
566  eContrib->ElostToDeadMaterials = deadEnergy * GeV2MeV;
567  eContrib->totalRange = totalRange;
568  if (totalRange > 0)
569  eContrib->deadFraction = deadRange / totalRange;
570 
571 
572  double allHitEnergy = std::accumulate(allHits.begin(), allHits.end(), 0.0, [&](double runningSum, const auto & hit) {return runningSum + hit.GetEdep(); });
573 
574  mf::LogDebug("EnergyStandardCandles") << "Michel E = " << eContrib->trueFourP.E() << " MeV).";
575  mf::LogDebug("EnergyStandardCandles") << " Total energy in " << allHits.size() << " hits: " << GeV2MeV * allHitEnergy;
576 // mf::LogDebug("EnergyStandardCandles") << " Total energy in " << onlyEhits.size() << " hits directly from Michel E: " << GeV2MeV * EhitEnergy << std::endl;
577  mf::LogDebug("EnergyStandardCandles") << " End energy of Michel e: " << trueParticle->EndE() * GeV2MeV << " MeV";
578 // mf::LogDebug("EnergyStandardCandles") << " Total energy of children: " << childEnergy * GeV2MeV << " MeV" << std::endl;
579  mf::LogDebug("EnergyStandardCandles") << " Invisible energy: " << (trueParticle->E() - allHitEnergy) * GeV2MeV << " MeV";
580  mf::LogDebug("EnergyStandardCandles") << " (estimated energy lost in inactive material: " << deadEnergy * GeV2MeV << " MeV)";
581  mf::LogDebug("EnergyStandardCandles") << " (energy lost in Bremsstrahlung: " << bremEnergy * GeV2MeV << " MeV)";
582 
583  } // EnergyStandardCandles::FillTrueMichelInfo
584 
585  void EnergyStandardCandles::FillTruthInfo(const std::shared_ptr<StoppingMu> stoppingMu, MichelCandidate & michelCand)
586  {
587  michelCand.eContribs.clear();
588  for (const auto & geantTrackInfo : stoppingMu->michelTruth)
589  {
590  TrueEContrib * eContrib;
591  const sim::Particle * trueParticle = this->fBackTracker->TrackIDToParticle(geantTrackInfo.trackID);
592  const sim::Particle * parentParticle = NULL;
593  if (trueParticle->Process() != "Primary")
594  parentParticle = this->fBackTracker->TrackIDToParticle(trueParticle->Mother());
595 
596  if ((trueParticle->PdgCode() == 11 && trueParticle->Process() == "muMinusCaptureAtRest")
597  || (abs(trueParticle->PdgCode()) == 11 && trueParticle->Process() == "Decay" && abs(parentParticle->PdgCode()) == 13))
598  eContrib = new TrueMichel;
599  else if (parentParticle && trueParticle->Process() == "nCapture")
600  eContrib = new TrueNeutronCapture;
601  else
602  eContrib = new TrueEContrib;
603 
604  // photons from pi0 decay, neutron capture need special treatment
605  if (trueParticle->PdgCode() == 22)
606  this->FixPhotonAncestry(trueParticle, parentParticle, eContrib);
607 // else if (abs(trueParticle->PdgCode()) == 13)
608 // {
609 // std::cout << " muon. Process: " << trueParticle->Process();
610 // if (parentParticle)
611 // std::cout << "; parent PDG: " << parentParticle->PdgCode();
612 // std::cout << std::endl;
613 // }
614  eContrib->pdg = trueParticle->PdgCode();
615  eContrib->geantCreationProcName = trueParticle->Process();
616  eContrib->trueFourP = TLorentzVector(
617  trueParticle->Px() * GeV2MeV,
618  trueParticle->Py() * GeV2MeV,
619  trueParticle->Pz() * GeV2MeV,
620  trueParticle->E() * GeV2MeV
621  );
622  eContrib->contribEFrac = geantTrackInfo.energyFrac;
623  eContrib->fracCaptured = geantTrackInfo.energy/trueParticle->E();
624 
625  // if it really was a Michel electron, how much of it did we get?
626  // (aside: looks like GEANT4 combines muon capture and decay-in-orbit
627  // into a single process, 'muMinusCaptureAtRest'. in that process
628  // we wind up with the usual electron and two neutrinos, possibly along
629  // with some de-excitation photons, the decay time away from the muon's
630  // final time. this is in contrast to normal muon decay, in which
631  // the muon's lifetime is kept by the muon particle and the delta T
632  // between the electron and the muon is 0.)
633  this->FillTrueMichelInfo(trueParticle, dynamic_cast<TrueMichel*>(eContrib));
634 
635  michelCand.eContribs.push_back(eContrib);
636  } // for (geant particle contributor)
637 
638  } //EnergyStandardCandles::FillTruthInfo()
639 
641  {
642 
643  for (const auto & stoppingMu : this->fStoppingMus)
644  {
645  MichelCandidate & michelCand = *(this->fMichelCandidate); // save on some typing
646 
647  this->FillMichelRecoInfo(stoppingMu, michelCand);
648  this->FillTruthInfo(stoppingMu, michelCand);
649 
650  this->fOutTrees["Michel"]->Fill();
651  } // for (stoppingMu)
652 
653  } // EnergyStandardCandles::FillMichelENtuple()
654 
655  // --------------------------------------------
657  {
658  }
659 
661  {
662  }
663 
664  // --------------------------------------------
666  {
667 
668 
669  }
670 
672  {
673  }
674 
676 
677 } // calib namespace
double E(const int i=0) const
Definition: MCParticle.h:232
std::map< std::string, TTree * > fOutTrees
SubRunNumber_t subRun() const
Definition: Event.h:72
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
std::vector< std::shared_ptr< StoppingMu > > fStoppingMus
EnergyStandardCandles(const fhicl::ParameterSet &pset)
double Py(const int i=0) const
Definition: MCParticle.h:230
void FixPhotonAncestry(const sim::Particle *&trueParticle, const sim::Particle *parentParticle, TrueEContrib *eContrib)
Sometimes photon ancestries are more complicated than we need them to be.
#define location
double EndE() const
Definition: MCParticle.h:243
TGeoMaterial * Material(double x, double y, double z) const
const std::vector< cheat::TrackIDE > michelTruth
static const double GeV2MeV
void FillTruthInfo(const std::shared_ptr< StoppingMu > stoppingMu, MichelCandidate &michelCand)
double delta
Definition: runWimpSim.h:98
const simb::MCTrajectory & Trajectory() const
Definition: MCParticle.h:252
MichelCandidate * fMichelCandidate
used for filling output ntuple
double DistToEdgeXY(TVector3 vertex)
int Mother() const
Definition: MCParticle.h:212
double fFrontSize
maximum distance to front of detector to qualify as "front-entering"
iterator begin()
Definition: PtrVector.h:223
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 Px(const int i=0) const
Definition: MCParticle.h:229
Definition: event.h:19
void abs(TH1 *hist)
StoppingMu(const rb::Track &t, const std::vector< art::Ptr< me::TrkME > > &cs, const art::PtrVector< rb::CellHit > &dblCnt, const std::vector< cheat::TrackIDE > &truth)
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
DEFINE_ART_MODULE(TestTMapFile)
double ElostToBrem
if it was a true Michel electron, how much of its total energy was radiated away in Bremsstrahlung...
virtual TVector3 Start() const
Definition: Prong.h:73
std::string Process() const
Definition: MCParticle.h:214
const art::PtrVector< rb::CellHit > doubleCountedHits
void MaterialsBetweenPoints(const double *p1, const double *p2, std::vector< double > &ds, std::vector< const TGeoMaterial * > &mat) const
double dE
int NumberDaughters() const
Definition: MCParticle.h:216
int TrackId() const
Definition: MCParticle.h:209
int Daughter(const int i) const
Definition: MCParticle.cxx:112
TParameter< double > * fPOTStorage
container for POT that will go to output ntuple
string material
Definition: eplot.py:19
std::string geantCreationProcName
GEANT4&#39;s name for the process that created this particle.
bool isValid() const
Definition: Handle.h:189
std::string fNumiBeamLabel
name of NuMI spill information container in event record
Track finder for cosmic rays.
static constexpr double cm
Definition: SystemOfUnits.h:99
bool IsContained(const rb::Track &trk) const
std::string GetName(int i)
CDPStorage service.
void hits()
Definition: readHits.C:15
double DistToEdgeZ(TVector3 vertex)
double fracCaptured
how much of this particle&#39;s true energy is contained in the reconstructed Michel
std::vector< TrueEContrib * > eContribs
T get(std::string const &key) const
Definition: ParameterSet.h:231
std::string fScintillatorName
name of geometry TGeoMaterial corresponding to scintillator
static constexpr double MeV
int evt
iterator end()
Definition: PtrVector.h:237
double energy
Definition: plottest35.C:25
const std::vector< art::Ptr< me::TrkME > > michelEClusters
void reconfigure(const fhicl::ParameterSet &pset)
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
double contribEFrac
how much of the Michel candidate&#39;s energy is due to this particle
art::ServiceHandle< geo::LiveGeometry > fLiveGeom
const ana::Var wgt
double ElostToDeadMaterials
if it was a true Michel electron, how much of its total energy was deposited in inactive materials (l...
const sim::Particle * TrackIDToParticle(int const &id) const
Returns a pointer to the sim::Particle object corresponding to the given TrackID. ...
When a Michel candidate has a contribution from a true Michel electron, we want more info about it...
std::map< std::string, int > fHeaderInfo
run/subrun/gate for output ntuple
bool IsMuon(const rb::Track &trk) const
std::vector< double > hitEnergies
EventNumber_t event() const
Definition: EventID.h:116
art::ServiceHandle< art::TFileService > fFileService
std::vector< TrackIDE > HitsToTrackIDE(const std::vector< const rb::CellHit * > &hits, bool useBirksE=false) const
Returns vector of TrackIDE structs contributing to the given collection of hits.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
iterator insert(iterator position, Ptr< U > const &p)
T * make(ARGS...args) const
std::vector< std::shared_ptr< Pi0 > > fPi0s
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
long int nCaptureNucleus
what was the target nucleus? (default: 0)
std::string fTrackLabel
name of tracks container in event record
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
Definition: structs.h:12
static constexpr double mm
Definition: SystemOfUnits.h:95
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:133
std::string fHitLabel
name of hits container in event record
double Pz(const int i=0) const
Definition: MCParticle.h:231
void FindStoppingMuons(const art::Event &evt)
void FillMichelRecoInfo(const std::shared_ptr< StoppingMu > stoppingMu, MichelCandidate &michelCand)
Utility type to hold info about stopping muons.
art::ServiceHandle< cheat::BackTracker > fBackTracker
Information about a reconstructed Michel electron candidate.
Information about one of the true particles that contribute to a candidate Michel electron...
double fMinDistToEdge
minimum distance in cm to edge of detector to qualify as "contained"
std::string fMichelLabel
name of Michel electron container in event record
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
RunNumber_t run() const
Definition: Event.h:77
void FillTrueMichelInfo(const sim::Particle *trueParticle, TrueMichel *eContrib)
static constexpr double GeV
Definition: fwd.h:28
double spillpot
POT for spill normalized by 10^12.
Definition: SpillData.h:26
std::string fGeneratorLabel
name of generator information label in event record
EventID id() const
Definition: Event.h:56
Encapsulate the geometry of one entire detector (near, far, ndos)
When a Michel candidate has a contribution from a neutron capture, that&#39;s also interesting.
art::ServiceHandle< geo::Geometry > fGeom