LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
Public Types | Public Member Functions | Protected Member Functions | Private Member Functions | Private Attributes | List of all members
calo::Calorimetry Class Reference

Estimates the energy deposited by reconstructed tracks. More...

Inheritance diagram for calo::Calorimetry:
art::EDProducer art::detail::Producer art::detail::LegacyModule art::Modifier art::ModuleBase art::ProductRegistryHelper

Public Types

using ModuleType = EDProducer
 
template<typename UserConfig , typename KeysToIgnore = void>
using Table = Modifier::Table< UserConfig, KeysToIgnore >
 

Public Member Functions

 Calorimetry (fhicl::ParameterSet const &pset)
 
void doBeginJob (SharedResources const &resources)
 
void doEndJob ()
 
void doRespondToOpenInputFile (FileBlock const &fb)
 
void doRespondToCloseInputFile (FileBlock const &fb)
 
void doRespondToOpenOutputFiles (FileBlock const &fb)
 
void doRespondToCloseOutputFiles (FileBlock const &fb)
 
bool doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doEndRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
 
void fillProductDescriptions ()
 
void registerProducts (ProductDescriptions &productsToRegister)
 
ModuleDescription const & moduleDescription () const
 
void setModuleDescription (ModuleDescription const &)
 
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables () const
 
void sortConsumables (std::string const &current_process_name)
 
std::unique_ptr< Worker > makeWorker (WorkerParams const &wp)
 
template<typename T , BranchType BT>
ViewToken< T > consumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< T > mayConsumeView (InputTag const &tag)
 

Protected Member Functions

ConsumesCollector & consumesCollector ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 

Private Member Functions

void produce (art::Event &evt) override
 
void ReadCaloTree ()
 
bool BeginsOnBoundary (art::Ptr< recob::Track > lar_track)
 
bool EndsOnBoundary (art::Ptr< recob::Track > lar_track)
 
void GetPitch (detinfo::DetectorPropertiesData const &det_prop, art::Ptr< recob::Hit > const &hit, std::vector< double > const &trkx, std::vector< double > const &trky, std::vector< double > const &trkz, std::vector< double > const &trkw, std::vector< double > const &trkx0, double *xyz3d, double &pitch, double TickT0)
 

Private Attributes

std::string fTrackModuleLabel
 
std::string fSpacePointModuleLabel
 
std::string fT0ModuleLabel
 
bool fUseArea
 
bool fSCE
 
std::optional< double > fNotOnTrackZcut
 
bool fFlipTrack_dQdx
 
CalorimetryAlg caloAlg
 
int fnsps
 
std::vector< int > fwire
 
std::vector< double > ftime
 
std::vector< double > fstime
 
std::vector< double > fetime
 
std::vector< double > fMIPs
 
std::vector< double > fdQdx
 
std::vector< double > fdEdx
 
std::vector< double > fResRng
 
std::vector< float > fpitch
 
std::vector< TVector3 > fXYZ
 
std::vector< size_t > fHitIndex
 

Detailed Description

Estimates the energy deposited by reconstructed tracks.

Output

Configuration

Note
This documentation is grossly incomplete.
  • NotOnTrackZcut (real, optional): if specified, hits associated to all trajectory points whose z coordinate is below NotOnTrackZcut value (including electric field distortion correction if enabled) are excluded from the calorimetry. The value is specified as absolute z coordinate in world reference frame, in centimeters. The legacy value of this cut was hard coded to -100.0 cm.

Definition at line 88 of file Calorimetry_module.cc.

Member Typedef Documentation

using art::EDProducer::ModuleType = EDProducer
inherited

Definition at line 17 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>
using art::detail::Producer::Table = Modifier::Table<UserConfig, KeysToIgnore>
inherited

Definition at line 26 of file Producer.h.

Constructor & Destructor Documentation

calo::Calorimetry::Calorimetry ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 140 of file Calorimetry_module.cc.

References caloAlg, fFlipTrack_dQdx, fNotOnTrackZcut, fSCE, fSpacePointModuleLabel, fT0ModuleLabel, fTrackModuleLabel, and fUseArea.

141  : EDProducer{pset}
142  , fTrackModuleLabel(pset.get<std::string>("TrackModuleLabel"))
143  , fSpacePointModuleLabel(pset.get<std::string>("SpacePointModuleLabel"))
144  , fT0ModuleLabel(pset.get<std::string>("T0ModuleLabel"))
145  , fUseArea(pset.get<bool>("UseArea"))
146  , fSCE(pset.get<bool>("CorrectSCE"))
147  , fFlipTrack_dQdx(pset.get<bool>("FlipTrack_dQdx", true))
148  , caloAlg(pset.get<fhicl::ParameterSet>("CaloAlg"))
149 {
150 
151  if (pset.has_key("NotOnTrackZcut")) fNotOnTrackZcut = pset.get<double>("NotOnTrackZcut");
152 
153  produces<std::vector<anab::Calorimetry>>();
154  produces<art::Assns<recob::Track, anab::Calorimetry>>();
155 }
calo::Calorimetry::fNotOnTrackZcut
std::optional< double > fNotOnTrackZcut
Definition: Calorimetry_module.cc:116
art::EDProducer::EDProducer
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
calo::Calorimetry::fSpacePointModuleLabel
std::string fSpacePointModuleLabel
Definition: Calorimetry_module.cc:112
calo::Calorimetry::caloAlg
CalorimetryAlg caloAlg
Definition: Calorimetry_module.cc:120
calo::Calorimetry::fTrackModuleLabel
std::string fTrackModuleLabel
Definition: Calorimetry_module.cc:111

Member Function Documentation

bool calo::Calorimetry::BeginsOnBoundary ( art::Ptr< recob::Track >  lar_track)
private
template<typename T , BranchType BT>
ProductToken< T > art::ModuleBase::consumes ( InputTag const &  tag)
protectedinherited

Definition at line 61 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumes().

62  {
63  return collector_.consumes<T, BT>(tag);
64  }
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:56
art::ConsumesCollector::consumes
ProductToken< T > consumes(InputTag const &)
Definition: ConsumesCollector.h:50
ConsumesCollector & art::ModuleBase::consumesCollector ( )
protectedinherited

Definition at line 57 of file ModuleBase.cc.

References art::ModuleBase::collector_.

58  {
59  return collector_;
60  }
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:56
template<typename T , BranchType BT>
void art::ModuleBase::consumesMany ( )
protectedinherited

Definition at line 75 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumesMany().

76  {
77  collector_.consumesMany<T, BT>();
78  }
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:56
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::ModuleBase::consumesView ( InputTag const &  )
protectedinherited
template<typename T , BranchType BT>
ViewToken<T> art::ModuleBase::consumesView ( InputTag const &  tag)
inherited

Definition at line 68 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumesView().

69  {
70  return collector_.consumesView<T, BT>(tag);
71  }
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:56
art::ConsumesCollector::consumesView
ViewToken< Element > consumesView(InputTag const &)
void art::detail::Producer::doBeginJob ( SharedResources const &  resources)
inherited

Definition at line 22 of file Producer.cc.

References art::detail::Producer::beginJobWithFrame(), and art::detail::Producer::setupQueues().

23  {
24  setupQueues(resources);
25  ProcessingFrame const frame{ScheduleID{}};
26  beginJobWithFrame(frame);
27  }
art::detail::Producer::setupQueues
virtual void setupQueues(SharedResources const &)=0
art::detail::Producer::beginJobWithFrame
virtual void beginJobWithFrame(ProcessingFrame const &)=0
bool art::detail::Producer::doBeginRun ( RunPrincipal &  rp,
ModuleContext const &  mc 
)
inherited

Definition at line 65 of file Producer.cc.

References art::detail::Producer::beginRunWithFrame(), art::RangeSet::forRun(), art::RunPrincipal::makeRun(), r, art::RunPrincipal::runID(), and art::ModuleContext::scheduleID().

66  {
67  auto r = rp.makeRun(mc, RangeSet::forRun(rp.runID()));
68  ProcessingFrame const frame{mc.scheduleID()};
69  beginRunWithFrame(r, frame);
70  r.commitProducts();
71  return true;
72  }
r
TRandom r
Definition: spectrum.C:23
art::detail::Producer::beginRunWithFrame
virtual void beginRunWithFrame(Run &, ProcessingFrame const &)=0
art::RangeSet::forRun
static RangeSet forRun(RunID)
Definition: RangeSet.cc:51
bool art::detail::Producer::doBeginSubRun ( SubRunPrincipal &  srp,
ModuleContext const &  mc 
)
inherited

Definition at line 85 of file Producer.cc.

References art::detail::Producer::beginSubRunWithFrame(), art::RangeSet::forSubRun(), art::SubRunPrincipal::makeSubRun(), art::ModuleContext::scheduleID(), and art::SubRunPrincipal::subRunID().

86  {
87  auto sr = srp.makeSubRun(mc, RangeSet::forSubRun(srp.subRunID()));
88  ProcessingFrame const frame{mc.scheduleID()};
89  beginSubRunWithFrame(sr, frame);
90  sr.commitProducts();
91  return true;
92  }
art::detail::Producer::beginSubRunWithFrame
virtual void beginSubRunWithFrame(SubRun &, ProcessingFrame const &)=0
art::RangeSet::forSubRun
static RangeSet forSubRun(SubRunID)
Definition: RangeSet.cc:57
void art::detail::Producer::doEndJob ( )
inherited

Definition at line 30 of file Producer.cc.

References art::detail::Producer::endJobWithFrame().

31  {
32  ProcessingFrame const frame{ScheduleID{}};
33  endJobWithFrame(frame);
34  }
art::detail::Producer::endJobWithFrame
virtual void endJobWithFrame(ProcessingFrame const &)=0
bool art::detail::Producer::doEndRun ( RunPrincipal &  rp,
ModuleContext const &  mc 
)
inherited

Definition at line 75 of file Producer.cc.

References art::detail::Producer::endRunWithFrame(), art::RunPrincipal::makeRun(), r, art::ModuleContext::scheduleID(), and art::Principal::seenRanges().

76  {
77  auto r = rp.makeRun(mc, rp.seenRanges());
78  ProcessingFrame const frame{mc.scheduleID()};
79  endRunWithFrame(r, frame);
80  r.commitProducts();
81  return true;
82  }
r
TRandom r
Definition: spectrum.C:23
art::detail::Producer::endRunWithFrame
virtual void endRunWithFrame(Run &, ProcessingFrame const &)=0
bool art::detail::Producer::doEndSubRun ( SubRunPrincipal &  srp,
ModuleContext const &  mc 
)
inherited

Definition at line 95 of file Producer.cc.

References art::detail::Producer::endSubRunWithFrame(), art::SubRunPrincipal::makeSubRun(), art::ModuleContext::scheduleID(), and art::Principal::seenRanges().

96  {
97  auto sr = srp.makeSubRun(mc, srp.seenRanges());
98  ProcessingFrame const frame{mc.scheduleID()};
99  endSubRunWithFrame(sr, frame);
100  sr.commitProducts();
101  return true;
102  }
art::detail::Producer::endSubRunWithFrame
virtual void endSubRunWithFrame(SubRun &, ProcessingFrame const &)=0
bool art::detail::Producer::doEvent ( EventPrincipal &  ep,
ModuleContext const &  mc,
std::atomic< std::size_t > &  counts_run,
std::atomic< std::size_t > &  counts_passed,
std::atomic< std::size_t > &  counts_failed 
)
inherited

Definition at line 105 of file Producer.cc.

References art::detail::Producer::checkPutProducts_, e, art::EventPrincipal::makeEvent(), art::detail::Producer::produceWithFrame(), and art::ModuleContext::scheduleID().

110  {
111  auto e = ep.makeEvent(mc);
112  ++counts_run;
113  ProcessingFrame const frame{mc.scheduleID()};
114  produceWithFrame(e, frame);
115  e.commitProducts(checkPutProducts_, &expectedProducts<InEvent>());
116  ++counts_passed;
117  return true;
118  }
art::detail::Producer::checkPutProducts_
bool const checkPutProducts_
Definition: Producer.h:70
e
Float_t e
Definition: plot.C:35
art::detail::Producer::produceWithFrame
virtual void produceWithFrame(Event &, ProcessingFrame const &)=0
void art::detail::Producer::doRespondToCloseInputFile ( FileBlock const &  fb)
inherited

Definition at line 44 of file Producer.cc.

References art::detail::Producer::respondToCloseInputFileWithFrame().

45  {
46  ProcessingFrame const frame{ScheduleID{}};
47  respondToCloseInputFileWithFrame(fb, frame);
48  }
art::detail::Producer::respondToCloseInputFileWithFrame
virtual void respondToCloseInputFileWithFrame(FileBlock const &, ProcessingFrame const &)=0
fb
TFile fb("Li6.root")
void art::detail::Producer::doRespondToCloseOutputFiles ( FileBlock const &  fb)
inherited

Definition at line 58 of file Producer.cc.

References art::detail::Producer::respondToCloseOutputFilesWithFrame().

59  {
60  ProcessingFrame const frame{ScheduleID{}};
61  respondToCloseOutputFilesWithFrame(fb, frame);
62  }
art::detail::Producer::respondToCloseOutputFilesWithFrame
virtual void respondToCloseOutputFilesWithFrame(FileBlock const &, ProcessingFrame const &)=0
fb
TFile fb("Li6.root")
void art::detail::Producer::doRespondToOpenInputFile ( FileBlock const &  fb)
inherited

Definition at line 37 of file Producer.cc.

References art::detail::Producer::respondToOpenInputFileWithFrame().

38  {
39  ProcessingFrame const frame{ScheduleID{}};
40  respondToOpenInputFileWithFrame(fb, frame);
41  }
art::detail::Producer::respondToOpenInputFileWithFrame
virtual void respondToOpenInputFileWithFrame(FileBlock const &, ProcessingFrame const &)=0
fb
TFile fb("Li6.root")
void art::detail::Producer::doRespondToOpenOutputFiles ( FileBlock const &  fb)
inherited

Definition at line 51 of file Producer.cc.

References art::detail::Producer::respondToOpenOutputFilesWithFrame().

52  {
53  ProcessingFrame const frame{ScheduleID{}};
54  respondToOpenOutputFilesWithFrame(fb, frame);
55  }
art::detail::Producer::respondToOpenOutputFilesWithFrame
virtual void respondToOpenOutputFilesWithFrame(FileBlock const &, ProcessingFrame const &)=0
fb
TFile fb("Li6.root")
bool calo::Calorimetry::EndsOnBoundary ( art::Ptr< recob::Track >  lar_track)
private
void art::Modifier::fillProductDescriptions ( )
inherited

Definition at line 10 of file Modifier.cc.

References art::ProductRegistryHelper::fillDescriptions(), and art::ModuleBase::moduleDescription().

11  {
12  ProductRegistryHelper::fillDescriptions(moduleDescription());
13  }
art::ProductRegistryHelper::fillDescriptions
void fillDescriptions(ModuleDescription const &md)
Definition: ProductRegistryHelper.cc:48
art::ModuleBase::moduleDescription
ModuleDescription const & moduleDescription() const
Definition: ModuleBase.cc:13
std::array< std::vector< ProductInfo >, NumBranchTypes > const & art::ModuleBase::getConsumables ( ) const
inherited

Definition at line 43 of file ModuleBase.cc.

References art::ModuleBase::collector_, and art::ConsumesCollector::getConsumables().

44  {
45  return collector_.getConsumables();
46  }
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:56
art::ConsumesCollector::getConsumables
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables() const
Definition: ConsumesCollector.cc:13
void calo::Calorimetry::GetPitch ( detinfo::DetectorPropertiesData const &  det_prop,
art::Ptr< recob::Hit > const &  hit,
std::vector< double > const &  trkx,
std::vector< double > const &  trky,
std::vector< double > const &  trkz,
std::vector< double > const &  trkw,
std::vector< double > const &  trkx0,
double *  xyz3d,
double &  pitch,
double  TickT0 
)
private

Definition at line 675 of file Calorimetry_module.cc.

References geo::PlaneID::asPlaneID(), detinfo::DetectorPropertiesData::ConvertTicksToX(), DEFINE_ART_MODULE, dir, fSCE, recob::Hit::PeakTime(), geo::GeometryCore::Plane(), geo::PlaneGeo::ThetaZ(), geo::TPCID::TPC, geo::WireID::Wire, recob::Hit::WireID(), and geo::GeometryCore::WirePitch().

Referenced by produce().

685 {
686  // Get 3d coordinates and track pitch for each hit
687  // Find 5 nearest space points and determine xyz and curvature->track pitch
688 
689  // Get services
690  art::ServiceHandle<geo::Geometry const> geom;
691  auto const* sce = lar::providerFrom<spacecharge::SpaceChargeService>();
692 
693  //save distance to each spacepoint sorted by distance
694  std::map<double, size_t> sptmap;
695  //save the sign of distance
696  std::map<size_t, int> sptsignmap;
697 
698  double wire_pitch = geom->WirePitch(geo::WireID(0, 0, 0, 0));
699 
700  double t0 = hit->PeakTime() - TickT0;
701  double x0 = det_prop.ConvertTicksToX(t0, hit->WireID().asPlaneID());
702  double w0 = hit->WireID().Wire;
703 
704  for (size_t i = 0; i < trkx.size(); ++i) {
705  double distance = cet::sum_of_squares((trkw[i] - w0) * wire_pitch, trkx0[i] - x0);
706  if (distance > 0) distance = sqrt(distance);
707  sptmap.insert(std::pair<double, size_t>(distance, i));
708  if (w0 - trkw[i] > 0)
709  sptsignmap.emplace(i, 1);
710  else
711  sptsignmap.emplace(i, -1);
712  }
713 
714  //x,y,z vs distance
715  std::vector<double> vx;
716  std::vector<double> vy;
717  std::vector<double> vz;
718  std::vector<double> vs;
719 
720  double kx = 0, ky = 0, kz = 0;
721 
722  int np = 0;
723  for (auto isp = sptmap.begin(); isp != sptmap.end(); isp++) {
724  double xyz[3];
725  xyz[0] = trkx[isp->second];
726  xyz[1] = trky[isp->second];
727  xyz[2] = trkz[isp->second];
728 
729  double distancesign = sptsignmap[isp->second];
730  if (np == 0 && isp->first > 30) { // hit not on track
731  xyz3d[0] = std::numeric_limits<double>::lowest();
732  xyz3d[1] = std::numeric_limits<double>::lowest();
733  xyz3d[2] = std::numeric_limits<double>::lowest();
734  pitch = -1;
735  return;
736  }
737  if (np < 5) {
738  vx.push_back(xyz[0]);
739  vy.push_back(xyz[1]);
740  vz.push_back(xyz[2]);
741  vs.push_back(isp->first * distancesign);
742  }
743  else {
744  break;
745  }
746  np++;
747  }
748  if (np >= 2) { // at least two points
749  TGraph* xs = new TGraph(np, &vs[0], &vx[0]);
750  try {
751  if (np > 2) { xs->Fit("pol2", "Q"); }
752  else {
753  xs->Fit("pol1", "Q");
754  }
755  TF1* pol = 0;
756  if (np > 2)
757  pol = (TF1*)xs->GetFunction("pol2");
758  else
759  pol = (TF1*)xs->GetFunction("pol1");
760  xyz3d[0] = pol->Eval(0);
761  kx = pol->GetParameter(1);
762  }
763  catch (...) {
764  mf::LogWarning("Calorimetry::GetPitch") << "Fitter failed";
765  xyz3d[0] = vx[0];
766  }
767  delete xs;
768  TGraph* ys = new TGraph(np, &vs[0], &vy[0]);
769  try {
770  if (np > 2) { ys->Fit("pol2", "Q"); }
771  else {
772  ys->Fit("pol1", "Q");
773  }
774  TF1* pol = 0;
775  if (np > 2)
776  pol = (TF1*)ys->GetFunction("pol2");
777  else
778  pol = (TF1*)ys->GetFunction("pol1");
779  xyz3d[1] = pol->Eval(0);
780  ky = pol->GetParameter(1);
781  }
782  catch (...) {
783  mf::LogWarning("Calorimetry::GetPitch") << "Fitter failed";
784  xyz3d[1] = vy[0];
785  }
786  delete ys;
787  TGraph* zs = new TGraph(np, &vs[0], &vz[0]);
788  try {
789  if (np > 2) { zs->Fit("pol2", "Q"); }
790  else {
791  zs->Fit("pol1", "Q");
792  }
793  TF1* pol = 0;
794  if (np > 2)
795  pol = (TF1*)zs->GetFunction("pol2");
796  else
797  pol = (TF1*)zs->GetFunction("pol1");
798  xyz3d[2] = pol->Eval(0);
799  kz = pol->GetParameter(1);
800  }
801  catch (...) {
802  mf::LogWarning("Calorimetry::GetPitch") << "Fitter failed";
803  xyz3d[2] = vz[0];
804  }
805  delete zs;
806  }
807  else if (np) {
808  xyz3d[0] = vx[0];
809  xyz3d[1] = vy[0];
810  xyz3d[2] = vz[0];
811  }
812  else {
813  xyz3d[0] = std::numeric_limits<double>::lowest();
814  xyz3d[1] = std::numeric_limits<double>::lowest();
815  xyz3d[2] = std::numeric_limits<double>::lowest();
816  pitch = -1;
817  return;
818  }
819  pitch = -1;
820  if (kx * kx + ky * ky + kz * kz) {
821  double tot = sqrt(kx * kx + ky * ky + kz * kz);
822  kx /= tot;
823  ky /= tot;
824  kz /= tot;
825  //get pitch
826  double wirePitch = geom->WirePitch(hit->WireID().asPlaneID());
827  double angleToVert =
828  geom->Plane(hit->WireID().asPlaneID()).Wire(0).ThetaZ(false) - 0.5 * TMath::Pi();
829  double cosgamma = TMath::Abs(TMath::Sin(angleToVert) * ky + TMath::Cos(angleToVert) * kz);
830  if (cosgamma > 0) pitch = wirePitch / cosgamma;
831 
832  //Correct for SCE
833  geo::Vector_t posOffsets = {0., 0., 0.};
834  geo::Vector_t dirOffsets = {0., 0., 0.};
835  if (sce->EnableCalSpatialSCE() && fSCE)
836  posOffsets =
837  sce->GetCalPosOffsets(geo::Point_t{xyz3d[0], xyz3d[1], xyz3d[2]}, hit->WireID().TPC);
838  if (sce->EnableCalSpatialSCE() && fSCE)
839  dirOffsets = sce->GetCalPosOffsets(
840  geo::Point_t{xyz3d[0] + pitch * kx, xyz3d[1] + pitch * ky, xyz3d[2] + pitch * kz},
841  hit->WireID().TPC);
842 
843  xyz3d[0] = xyz3d[0] - posOffsets.X();
844  xyz3d[1] = xyz3d[1] + posOffsets.Y();
845  xyz3d[2] = xyz3d[2] + posOffsets.Z();
846 
847  // Correct pitch for SCE
848  TVector3 dir = {pitch * kx - dirOffsets.X() + posOffsets.X(),
849  pitch * ky + dirOffsets.Y() - posOffsets.Y(),
850  pitch * kz + dirOffsets.Z() - posOffsets.Z()};
851  pitch = dir.Mag();
852  }
853 }
t0
code to link reconstructed objects back to the MC truth information
Definition: DirectHitParticleAssns_tool.cc:18
geo::Vector_t
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:160
geo::WireID::Wire
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:563
recob::Hit::WireID
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:280
geo::PlaneGeo::ThetaZ
double ThetaZ() const
Angle of the wires from positive z axis; .
Definition: PlaneGeo.cxx:646
geo::GeometryCore::Plane
PlaneGeo const & Plane(PlaneID const &planeid) const
Returns the specified wire.
Definition: GeometryCore.h:1034
geo::PlaneID::asPlaneID
constexpr PlaneID const & asPlaneID() const
Conversion to PlaneID (for convenience of notation).
Definition: geo_types.h:520
geo::Point_t
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:180
recob::Hit::PeakTime
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:220
dir
TDirectory * dir
Definition: macro.C:5
mf::LogWarning
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Definition: MessageLogger.h:210
geo::TPCID::TPC
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:399
geo::GeometryCore::WirePitch
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
Definition: GeometryCore.cxx:763
std::unique_ptr< Worker > art::ModuleBase::makeWorker ( WorkerParams const &  wp)
inherited

Definition at line 37 of file ModuleBase.cc.

References art::ModuleBase::doMakeWorker(), and art::NumBranchTypes.

38  {
39  return doMakeWorker(wp);
40  }
art::ModuleBase::doMakeWorker
virtual std::unique_ptr< Worker > doMakeWorker(WorkerParams const &wp)=0
template<typename T , BranchType BT>
ProductToken< T > art::ModuleBase::mayConsume ( InputTag const &  tag)
protectedinherited

Definition at line 82 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsume().

83  {
84  return collector_.mayConsume<T, BT>(tag);
85  }
art::ConsumesCollector::mayConsume
ProductToken< T > mayConsume(InputTag const &)
Definition: ConsumesCollector.h:75
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:56
template<typename T , BranchType BT>
void art::ModuleBase::mayConsumeMany ( )
protectedinherited

Definition at line 96 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsumeMany().

97  {
98  collector_.mayConsumeMany<T, BT>();
99  }
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:56
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::ModuleBase::mayConsumeView ( InputTag const &  )
protectedinherited
template<typename T , BranchType BT>
ViewToken<T> art::ModuleBase::mayConsumeView ( InputTag const &  tag)
inherited

Definition at line 89 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsumeView().

90  {
91  return collector_.mayConsumeView<T, BT>(tag);
92  }
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:56
art::ConsumesCollector::mayConsumeView
ViewToken< Element > mayConsumeView(InputTag const &)
ModuleDescription const & art::ModuleBase::moduleDescription ( ) const
inherited

Definition at line 13 of file ModuleBase.cc.

References art::errors::LogicError.

Referenced by art::OutputModule::doRespondToOpenInputFile(), art::OutputModule::doWriteEvent(), art::Modifier::fillProductDescriptions(), art::OutputModule::makePlugins_(), art::OutputWorker::OutputWorker(), reco::shower::LArPandoraModularShowerCreation::produce(), art::Modifier::registerProducts(), and art::OutputModule::registerProducts().

14  {
15  if (md_.has_value()) {
16  return *md_;
17  }
18 
19  throw Exception{errors::LogicError,
20  "There was an error while calling moduleDescription().\n"}
21  << "The moduleDescription() base-class member function cannot be called\n"
22  "during module construction. To determine which module is "
23  "responsible\n"
24  "for calling it, find the '<module type>:<module "
25  "label>@Construction'\n"
26  "tag in the message prefix above. Please contact artists@fnal.gov\n"
27  "for guidance.\n";
28  }
art::Exception
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
art::ModuleBase::md_
std::optional< ModuleDescription > md_
Definition: ModuleBase.h:55
void calo::Calorimetry::produce ( art::Event &  evt)
overrideprivatevirtual

Implements art::EDProducer.

Definition at line 158 of file Calorimetry_module.cc.

References util::abs(), lar::util::absDiff(), caloAlg, recob::tracking::convertCollToPoint(), util::CreateAssn(), geo::CryostatID::Cryostat, art::FindManyP< ProdB, Data >::data(), tca::dEdx(), calo::CalorimetryAlg::dEdx_AMP(), calo::CalorimetryAlg::dEdx_AREA(), dir, e, fdEdx, fdQdx, fetime, fFlipTrack_dQdx, fHitIndex, art::fill_ptr_vector(), geo::GeometryCore::FindTPCAtPosition(), fMIPs, fNotOnTrackZcut, fnsps, fpitch, fResRng, fSCE, fSpacePointModuleLabel, fstime, fT0ModuleLabel, ftime, fTrackModuleLabel, fUseArea, fwire, fXYZ, art::ProductRetriever::getByLabel(), GetPitch(), hits(), geo::CryostatID::isValid, util::kBogusD, MF_LOG_DEBUG, geo::GeometryCore::Nplanes(), geo::PlaneID::Plane, geo::GeometryCore::Plane(), geo::GeometryCore::PlaneWireToChannel(), art::Event::put(), detinfo::sampling_rate(), util::size(), geo::TPCID::TPC, lar::util::TrackPitchInView(), geo::PlaneGeo::View(), w, geo::GeometryCore::WireAngleToVertical(), geo::GeometryCore::WirePitch(), x, X, xx, Y, Z, and zz.

159 {
160  auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
161  auto const det_prop =
162  art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clock_data);
163  auto const* sce = lar::providerFrom<spacecharge::SpaceChargeService>();
164 
165  art::Handle<std::vector<recob::Track>> trackListHandle;
166  std::vector<art::Ptr<recob::Track>> tracklist;
167  if (evt.getByLabel(fTrackModuleLabel, trackListHandle))
168  art::fill_ptr_vector(tracklist, trackListHandle);
169 
170  // Get Geometry
171  art::ServiceHandle<geo::Geometry const> geom;
172 
173  // channel quality
174  lariov::ChannelStatusProvider const& channelStatus =
175  art::ServiceHandle<lariov::ChannelStatusService const>()->GetProvider();
176 
177  size_t nplanes = geom->Nplanes();
178 
179  //create anab::Calorimetry objects and make association with recob::Track
180  std::unique_ptr<std::vector<anab::Calorimetry>> calorimetrycol(
181  new std::vector<anab::Calorimetry>);
182  std::unique_ptr<art::Assns<recob::Track, anab::Calorimetry>> assn(
183  new art::Assns<recob::Track, anab::Calorimetry>);
184 
185  art::FindManyP<recob::Hit> fmht(trackListHandle, evt, fTrackModuleLabel);
186  art::FindManyP<recob::Hit, recob::TrackHitMeta> fmthm(
187  trackListHandle,
188  evt,
189  fTrackModuleLabel); //this has more information about hit-track association, only available in PMA for now
190  art::FindManyP<anab::T0> fmt0(trackListHandle, evt, fT0ModuleLabel);
191 
192  for (size_t trkIter = 0; trkIter < tracklist.size(); ++trkIter) {
193 
194  decltype(auto) larEnd = tracklist[trkIter]->Trajectory().End();
195 
196  // Some variables for the hit
197  float time; //hit time at maximum
198  float stime; //hit start time
199  float etime; //hit end time
200  uint32_t channel = 0; //channel number
201  unsigned int cstat = 0; //hit cryostat number
202  unsigned int tpc = 0; //hit tpc number
203  unsigned int wire = 0; //hit wire number
204  unsigned int plane = 0; //hit plane number
205 
206  std::vector<art::Ptr<recob::Hit>> allHits = fmht.at(trkIter);
207  double T0 = 0;
208  double TickT0 = 0;
209  if (fmt0.isValid()) {
210  std::vector<art::Ptr<anab::T0>> allT0 = fmt0.at(trkIter);
211  if (allT0.size()) T0 = allT0[0]->Time();
212  TickT0 = T0 / sampling_rate(clock_data);
213  }
214 
215  std::vector<std::vector<unsigned int>> hits(nplanes);
216 
217  art::FindManyP<recob::SpacePoint> fmspts(allHits, evt, fSpacePointModuleLabel);
218  for (size_t ah = 0; ah < allHits.size(); ++ah) {
219  hits[allHits[ah]->WireID().Plane].push_back(ah);
220  }
221  //get hits in each plane
222  for (unsigned int ipl = 0; ipl < nplanes; ++ipl) { //loop over all wire planes
223 
224  geo::PlaneID planeID; //(cstat,tpc,ipl);
225 
226  fwire.clear();
227  ftime.clear();
228  fstime.clear();
229  fetime.clear();
230  fMIPs.clear();
231  fdQdx.clear();
232  fdEdx.clear();
233  fpitch.clear();
234  fResRng.clear();
235  fXYZ.clear();
236  fHitIndex.clear();
237 
238  float Kin_En = 0.;
239  float Trk_Length = 0.;
240  std::vector<float> vdEdx;
241  std::vector<float> vresRange;
242  std::vector<float> vdQdx;
243  std::vector<float> deadwire; //residual range for dead wires
244  std::vector<TVector3> vXYZ;
245 
246  // Require at least 2 hits in this view
247  if (hits[ipl].size() < 2) {
248  if (hits[ipl].size() == 1) {
249  mf::LogWarning("Calorimetry")
250  << "Only one hit in plane " << ipl << " associated with track id " << trkIter;
251  }
252  calorimetrycol->emplace_back(util::kBogusD,
253  vdEdx,
254  vdQdx,
255  vresRange,
256  deadwire,
257  util::kBogusD,
258  fpitch,
259  recob::tracking::convertCollToPoint(vXYZ),
260  planeID);
261  util::CreateAssn(evt, *calorimetrycol, tracklist[trkIter], *assn);
262  continue;
263  }
264 
265  //range of wire signals
266  unsigned int wire0 = 100000;
267  unsigned int wire1 = 0;
268  double PIDA = 0;
269  int nPIDA = 0;
270 
271  // determine track direction. Fill residual range array
272  bool GoingDS = true;
273  // find the track direction by comparing US and DS charge BB
274  double USChg = 0;
275  double DSChg = 0;
276  // temp array holding distance betweeen space points
277  std::vector<double> spdelta;
278  fnsps = 0; // number of space points
279  std::vector<double> ChargeBeg;
280  std::stack<double> ChargeEnd;
281 
282  // find track pitch
283  double fTrkPitch = 0;
284  for (size_t itp = 0; itp < tracklist[trkIter]->NumberTrajectoryPoints(); ++itp) {
285 
286  const auto& pos = tracklist[trkIter]->LocationAtPoint(itp);
287  const auto& dir = tracklist[trkIter]->DirectionAtPoint(itp);
288 
289  geo::TPCID const tpcid = geom->FindTPCAtPosition(pos);
290  if (!tpcid.isValid) continue;
291 
292  try {
293  fTrkPitch =
294  lar::util::TrackPitchInView(*tracklist[trkIter], geom->Plane({tpcid, ipl}).View(), itp);
295 
296  //Correct for SCE
297  geo::Vector_t posOffsets = {0., 0., 0.};
298  geo::Vector_t dirOffsets = {0., 0., 0.};
299  if (sce->EnableCalSpatialSCE() && fSCE) {
300  posOffsets = sce->GetCalPosOffsets(pos, tpcid.TPC);
301  dirOffsets = sce->GetCalPosOffsets(pos + fTrkPitch * dir, tpcid.TPC);
302  }
303  TVector3 dir_corr = {fTrkPitch * dir.X() - dirOffsets.X() + posOffsets.X(),
304  fTrkPitch * dir.Y() + dirOffsets.Y() - posOffsets.Y(),
305  fTrkPitch * dir.Z() + dirOffsets.Z() - posOffsets.Z()};
306 
307  fTrkPitch = dir_corr.Mag();
308  }
309  catch (cet::exception& e) {
310  mf::LogWarning("Calorimetry")
311  << "caught exception " << e << "\n setting pitch (C) to " << util::kBogusD;
312  fTrkPitch = 0;
313  }
314  break;
315  }
316 
317  // find the separation between all space points
318  double xx = 0., yy = 0., zz = 0.;
319 
320  //save track 3d points
321  std::vector<double> trkx;
322  std::vector<double> trky;
323  std::vector<double> trkz;
324  std::vector<double> trkw;
325  std::vector<double> trkx0;
326  for (size_t i = 0; i < hits[ipl].size(); ++i) {
327  //Get space points associated with the hit
328  std::vector<art::Ptr<recob::SpacePoint>> sptv = fmspts.at(hits[ipl][i]);
329  for (size_t j = 0; j < sptv.size(); ++j) {
330 
331  double t = allHits[hits[ipl][i]]->PeakTime() -
332  TickT0; // Want T0 here? Otherwise ticks to x is wrong?
333  double x = det_prop.ConvertTicksToX(t,
334  allHits[hits[ipl][i]]->WireID().Plane,
335  allHits[hits[ipl][i]]->WireID().TPC,
336  allHits[hits[ipl][i]]->WireID().Cryostat);
337  double w = allHits[hits[ipl][i]]->WireID().Wire;
338  if (TickT0) {
339  trkx.push_back(sptv[j]->XYZ()[0] -
340  det_prop.ConvertTicksToX(TickT0,
341  allHits[hits[ipl][i]]->WireID().Plane,
342  allHits[hits[ipl][i]]->WireID().TPC,
343  allHits[hits[ipl][i]]->WireID().Cryostat));
344  }
345  else {
346  trkx.push_back(sptv[j]->XYZ()[0]);
347  }
348  trky.push_back(sptv[j]->XYZ()[1]);
349  trkz.push_back(sptv[j]->XYZ()[2]);
350  trkw.push_back(w);
351  trkx0.push_back(x);
352  }
353  }
354  for (size_t ihit = 0; ihit < hits[ipl].size();
355  ++ihit) { // loop over all hits on each wire plane
356 
357  if (!planeID.isValid) {
358  plane = allHits[hits[ipl][ihit]]->WireID().Plane;
359  tpc = allHits[hits[ipl][ihit]]->WireID().TPC;
360  cstat = allHits[hits[ipl][ihit]]->WireID().Cryostat;
361  planeID.Cryostat = cstat;
362  planeID.TPC = tpc;
363  planeID.Plane = plane;
364  planeID.isValid = true;
365  }
366 
367  wire = allHits[hits[ipl][ihit]]->WireID().Wire;
368  time = allHits[hits[ipl][ihit]]->PeakTime(); // What about here? T0
369  stime = allHits[hits[ipl][ihit]]->PeakTimeMinusRMS();
370  etime = allHits[hits[ipl][ihit]]->PeakTimePlusRMS();
371  const size_t& hitIndex = allHits[hits[ipl][ihit]].key();
372 
373  double charge = allHits[hits[ipl][ihit]]->PeakAmplitude();
374  if (fUseArea) charge = allHits[hits[ipl][ihit]]->Integral();
375  //get 3d coordinate and track pitch for the current hit
376  //not all hits are associated with space points, the method uses neighboring spacepts to interpolate
377  double xyz3d[3];
378  double pitch;
379  bool fBadhit = false;
380  if (fmthm.isValid()) {
381  auto vhit = fmthm.at(trkIter);
382  auto vmeta = fmthm.data(trkIter);
383  for (size_t ii = 0; ii < vhit.size(); ++ii) {
384  if (vhit[ii].key() == allHits[hits[ipl][ihit]].key()) {
385  if (vmeta[ii]->Index() == int_max_as_unsigned_int) {
386  fBadhit = true;
387  continue;
388  }
389  if (vmeta[ii]->Index() >= tracklist[trkIter]->NumberTrajectoryPoints()) {
390  throw cet::exception("Calorimetry_module.cc")
391  << "Requested track trajectory index " << vmeta[ii]->Index()
392  << " exceeds the total number of trajectory points "
393  << tracklist[trkIter]->NumberTrajectoryPoints() << " for track index " << trkIter
394  << ". Something is wrong with the track reconstruction. Please contact "
395  "tjyang@fnal.gov";
396  }
397  if (!tracklist[trkIter]->HasValidPoint(vmeta[ii]->Index())) {
398  fBadhit = true;
399  continue;
400  }
401 
402  //Correct location for SCE
403  geo::Point_t const loc = tracklist[trkIter]->LocationAtPoint(vmeta[ii]->Index());
404  geo::Vector_t locOffsets = {0., 0., 0.};
405  if (sce->EnableCalSpatialSCE() && fSCE)
406  locOffsets = sce->GetCalPosOffsets(loc, vhit[ii]->WireID().TPC);
407  xyz3d[0] = loc.X() - locOffsets.X();
408  xyz3d[1] = loc.Y() + locOffsets.Y();
409  xyz3d[2] = loc.Z() + locOffsets.Z();
410 
411  double angleToVert =
412  geom->WireAngleToVertical(vhit[ii]->View(), vhit[ii]->WireID().asPlaneID()) -
413  0.5 * ::util::pi<>();
414  const geo::Vector_t& dir = tracklist[trkIter]->DirectionAtPoint(vmeta[ii]->Index());
415  double cosgamma =
416  std::abs(std::sin(angleToVert) * dir.Y() + std::cos(angleToVert) * dir.Z());
417  if (cosgamma) { pitch = geom->WirePitch(vhit[ii]->View()) / cosgamma; }
418  else {
419  pitch = 0;
420  }
421 
422  //Correct pitch for SCE
423  geo::Vector_t dirOffsets = {0., 0., 0.};
424  if (sce->EnableCalSpatialSCE() && fSCE)
425  dirOffsets = sce->GetCalPosOffsets(geo::Point_t{loc.X() + pitch * dir.X(),
426  loc.Y() + pitch * dir.Y(),
427  loc.Z() + pitch * dir.Z()},
428  vhit[ii]->WireID().TPC);
429  const TVector3& dir_corr = {pitch * dir.X() - dirOffsets.X() + locOffsets.X(),
430  pitch * dir.Y() + dirOffsets.Y() - locOffsets.Y(),
431  pitch * dir.Z() + dirOffsets.Z() - locOffsets.Z()};
432 
433  pitch = dir_corr.Mag();
434 
435  break;
436  }
437  }
438  }
439  else
440  GetPitch(det_prop,
441  allHits[hits[ipl][ihit]],
442  trkx,
443  trky,
444  trkz,
445  trkw,
446  trkx0,
447  xyz3d,
448  pitch,
449  TickT0);
450 
451  if (fBadhit) continue;
452  if (fNotOnTrackZcut && (xyz3d[2] < fNotOnTrackZcut.value())) continue; //hit not on track
453  if (pitch <= 0) pitch = fTrkPitch;
454  if (!pitch) continue;
455 
456  if (fnsps == 0) {
457  xx = xyz3d[0];
458  yy = xyz3d[1];
459  zz = xyz3d[2];
460  spdelta.push_back(0);
461  }
462  else {
463  double dx = xyz3d[0] - xx;
464  double dy = xyz3d[1] - yy;
465  double dz = xyz3d[2] - zz;
466  spdelta.push_back(sqrt(dx * dx + dy * dy + dz * dz));
467  Trk_Length += spdelta.back();
468  xx = xyz3d[0];
469  yy = xyz3d[1];
470  zz = xyz3d[2];
471  }
472 
473  ChargeBeg.push_back(charge);
474  ChargeEnd.push(charge);
475 
476  double MIPs = charge;
477  double dQdx = MIPs / pitch;
478  double dEdx = 0;
479  if (fUseArea)
480  dEdx = caloAlg.dEdx_AREA(clock_data, det_prop, *allHits[hits[ipl][ihit]], pitch, T0);
481  else
482  dEdx = caloAlg.dEdx_AMP(clock_data, det_prop, *allHits[hits[ipl][ihit]], pitch, T0);
483 
484  Kin_En = Kin_En + dEdx * pitch;
485 
486  if (allHits[hits[ipl][ihit]]->WireID().Wire < wire0)
487  wire0 = allHits[hits[ipl][ihit]]->WireID().Wire;
488  if (allHits[hits[ipl][ihit]]->WireID().Wire > wire1)
489  wire1 = allHits[hits[ipl][ihit]]->WireID().Wire;
490 
491  fMIPs.push_back(MIPs);
492  fdEdx.push_back(dEdx);
493  fdQdx.push_back(dQdx);
494  fwire.push_back(wire);
495  ftime.push_back(time);
496  fstime.push_back(stime);
497  fetime.push_back(etime);
498  fpitch.push_back(pitch);
499  TVector3 v(xyz3d[0], xyz3d[1], xyz3d[2]);
500  fXYZ.push_back(v);
501  fHitIndex.push_back(hitIndex);
502  ++fnsps;
503  }
504  if (fnsps < 2) {
505  vdEdx.clear();
506  vdQdx.clear();
507  vresRange.clear();
508  deadwire.clear();
509  fpitch.clear();
510  calorimetrycol->push_back(anab::Calorimetry(util::kBogusD,
511  vdEdx,
512  vdQdx,
513  vresRange,
514  deadwire,
515  util::kBogusD,
516  fpitch,
517  recob::tracking::convertCollToPoint(vXYZ),
518  planeID));
519  util::CreateAssn(evt, *calorimetrycol, tracklist[trkIter], *assn);
520  continue;
521  }
522  for (int isp = 0; isp < fnsps; ++isp) {
523  if (isp > 3) break;
524  USChg += ChargeBeg[isp];
525  }
526  int countsp = 0;
527  while (!ChargeEnd.empty()) {
528  if (countsp > 3) break;
529  DSChg += ChargeEnd.top();
530  ChargeEnd.pop();
531  ++countsp;
532  }
533  if (fFlipTrack_dQdx) {
534  // Going DS if charge is higher at the end
535  GoingDS = (DSChg > USChg);
536  }
537  else {
538  // Use the track direction to determine the residual range
539  if (!fXYZ.empty()) {
540  TVector3 track_start(tracklist[trkIter]->Trajectory().Vertex().X(),
541  tracklist[trkIter]->Trajectory().Vertex().Y(),
542  tracklist[trkIter]->Trajectory().Vertex().Z());
543  TVector3 track_end(tracklist[trkIter]->Trajectory().End().X(),
544  tracklist[trkIter]->Trajectory().End().Y(),
545  tracklist[trkIter]->Trajectory().End().Z());
546 
547  if ((fXYZ[0] - track_start).Mag() + (fXYZ.back() - track_end).Mag() <
548  (fXYZ[0] - track_end).Mag() + (fXYZ.back() - track_start).Mag()) {
549  GoingDS = true;
550  }
551  else {
552  GoingDS = false;
553  }
554  }
555  }
556 
557  // determine the starting residual range and fill the array
558  fResRng.resize(fnsps);
559  if (fResRng.size() < 2 || spdelta.size() < 2) {
560  mf::LogWarning("Calorimetry")
561  << "fResrng.size() = " << fResRng.size() << " spdelta.size() = " << spdelta.size();
562  }
563  if (GoingDS) {
564  fResRng[fnsps - 1] = spdelta[fnsps - 1] / 2;
565  for (int isp = fnsps - 2; isp > -1; isp--) {
566  fResRng[isp] = fResRng[isp + 1] + spdelta[isp + 1];
567  }
568  }
569  else {
570  fResRng[0] = spdelta[1] / 2;
571  for (int isp = 1; isp < fnsps; isp++) {
572  fResRng[isp] = fResRng[isp - 1] + spdelta[isp];
573  }
574  }
575 
576  MF_LOG_DEBUG("CaloPrtHit") << " pt wire time ResRng MIPs pitch dE/dx Ai X Y Z\n";
577 
578  double Ai = -1;
579  for (int i = 0; i < fnsps; ++i) { //loop over all 3D points
580  vresRange.push_back(fResRng[i]);
581  vdEdx.push_back(fdEdx[i]);
582  vdQdx.push_back(fdQdx[i]);
583  vXYZ.push_back(fXYZ[i]);
584  if (i != 0 && i != fnsps - 1) { // ignore the first and last point
585  // Calculate PIDA
586  Ai = fdEdx[i] * pow(fResRng[i], 0.42);
587  nPIDA++;
588  PIDA += Ai;
589  }
590 
591  MF_LOG_DEBUG("CaloPrtHit")
592  << std::setw(4) << trkIter << std::setw(4) << ipl << std::setw(4) << i << std::setw(4)
593  << fwire[i] << std::setw(6) << (int)ftime[i]
594  << std::setiosflags(std::ios::fixed | std::ios::showpoint) << std::setprecision(2)
595  << std::setw(8) << fResRng[i] << std::setprecision(1) << std::setw(8) << fMIPs[i]
596  << std::setprecision(2) << std::setw(8) << fpitch[i] << std::setw(8) << fdEdx[i]
597  << std::setw(8) << Ai << std::setw(8) << fXYZ[i].x() << std::setw(8) << fXYZ[i].y()
598  << std::setw(8) << fXYZ[i].z() << "\n";
599  } // end looping over 3D points
600  if (nPIDA > 0) { PIDA = PIDA / (double)nPIDA; }
601  else {
602  PIDA = -1;
603  }
604  MF_LOG_DEBUG("CaloPrtTrk") << "Plane # " << ipl << "TrkPitch= " << std::setprecision(2)
605  << fTrkPitch << " nhits= " << fnsps << "\n"
606  << std::setiosflags(std::ios::fixed | std::ios::showpoint)
607  << "Trk Length= " << std::setprecision(1) << Trk_Length << " cm,"
608  << " KE calo= " << std::setprecision(1) << Kin_En << " MeV,"
609  << " PIDA= " << PIDA << "\n";
610 
611  // look for dead wires
612  for (unsigned int iw = wire0; iw < wire1 + 1; ++iw) {
613  plane = allHits[hits[ipl][0]]->WireID().Plane;
614  tpc = allHits[hits[ipl][0]]->WireID().TPC;
615  cstat = allHits[hits[ipl][0]]->WireID().Cryostat;
616  channel = geom->PlaneWireToChannel(geo::WireID{cstat, tpc, plane, iw});
617  if (channelStatus.IsBad(channel)) {
618  MF_LOG_DEBUG("Calorimetry") << "Found dead wire at Plane = " << plane << " Wire =" << iw;
619  unsigned int closestwire = 0;
620  unsigned int endwire = 0;
621  unsigned int dwire = 100000;
622  double mindis = 100000;
623  double goodresrange = 0;
624  for (size_t ihit = 0; ihit < hits[ipl].size(); ++ihit) {
625  channel = allHits[hits[ipl][ihit]]->Channel();
626  if (channelStatus.IsBad(channel)) continue;
627  // grab the space points associated with this hit
628  std::vector<art::Ptr<recob::SpacePoint>> sppv = fmspts.at(hits[ipl][ihit]);
629  if (sppv.size() < 1) continue;
630  // only use the first space point in the collection, really each hit
631  // should only map to 1 space point
632  const recob::Track::Point_t xyz{
633  sppv[0]->XYZ()[0], sppv[0]->XYZ()[1], sppv[0]->XYZ()[2]};
634  double dis1 = (larEnd - xyz).Mag2();
635  if (dis1) dis1 = std::sqrt(dis1);
636  if (dis1 < mindis) {
637  endwire = allHits[hits[ipl][ihit]]->WireID().Wire;
638  mindis = dis1;
639  }
640  if (lar::util::absDiff(wire, iw) < dwire) {
641  closestwire = allHits[hits[ipl][ihit]]->WireID().Wire;
642  dwire = lar::util::absDiff(allHits[hits[ipl][ihit]]->WireID().Wire, iw);
643  goodresrange = dis1;
644  }
645  }
646  if (closestwire) {
647  if (iw < endwire) {
648  deadwire.push_back(goodresrange + (int(closestwire) - int(iw)) * fTrkPitch);
649  }
650  else {
651  deadwire.push_back(goodresrange + (int(iw) - int(closestwire)) * fTrkPitch);
652  }
653  }
654  }
655  }
656  calorimetrycol->push_back(anab::Calorimetry(Kin_En,
657  vdEdx,
658  vdQdx,
659  vresRange,
660  deadwire,
661  Trk_Length,
662  fpitch,
663  recob::tracking::convertCollToPoint(vXYZ),
664  fHitIndex,
665  planeID));
666  util::CreateAssn(evt, *calorimetrycol, tracklist[trkIter], *assn);
667 
668  } //end looping over planes
669  } //end looping over tracks
670 
671  evt.put(std::move(calorimetrycol));
672  evt.put(std::move(assn));
673 }
x
Float_t x
Definition: compare.C:6
calo::Calorimetry::fdEdx
std::vector< double > fdEdx
Definition: Calorimetry_module.cc:129
calo::Calorimetry::fNotOnTrackZcut
std::optional< double > fNotOnTrackZcut
Definition: Calorimetry_module.cc:116
xx
Double_t xx
Definition: macro.C:12
recob
Reconstruction base classes.
Definition: FlashHypothesisCreator.h:18
geo::Vector_t
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:160
Y
Float_t Y
Definition: plot.C:37
geo::PlaneID
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
geo::CryostatID::isValid
bool isValid
Whether this ID points to a valid element.
Definition: geo_types.h:210
util::abs
constexpr auto abs(T v)
Returns the absolute value of the argument.
Definition: constexpr_math.h:40
geo::CryostatID::Cryostat
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:211
std
STL namespace.
zz
Double_t zz
Definition: plot.C:276
calo::Calorimetry::fXYZ
std::vector< TVector3 > fXYZ
Definition: Calorimetry_module.cc:132
calo::Calorimetry::fResRng
std::vector< double > fResRng
Definition: Calorimetry_module.cc:130
Z
Float_t Z
Definition: plot.C:37
art::Event::put
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
geo::GeometryCore::Plane
PlaneGeo const & Plane(PlaneID const &planeid) const
Returns the specified wire.
Definition: GeometryCore.h:1034
geo::GeometryCore::WireAngleToVertical
double WireAngleToVertical(View_t view, TPCID const &tpcid) const
Returns the angle of the wires in the specified view from vertical.
Definition: GeometryCore.cxx:781
geo::PlaneGeo::View
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:166
art::Ptr
Ptr(H, T) -> Ptr< detail::not_map_vector_t< typename H::element_type >>
hits
void hits()
Definition: readHits.C:15
geo::GeometryCore::FindTPCAtPosition
TPCID FindTPCAtPosition(Point_t const &point) const
Returns the ID of the TPC at specified location.
Definition: GeometryCore.cxx:326
if
if(nlines<=0)
Definition: comparison_ascii.C:19
calo::CalorimetryAlg::dEdx_AMP
double dEdx_AMP(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
Definition: CalorimetryAlg.cxx:60
geo::fhicl::WireID
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Definition: geo_types_fhicl.h:383
calo::Calorimetry::fHitIndex
std::vector< size_t > fHitIndex
Definition: Calorimetry_module.cc:133
calo::Calorimetry::fSpacePointModuleLabel
std::string fSpacePointModuleLabel
Definition: Calorimetry_module.cc:112
calo::Calorimetry::ftime
std::vector< double > ftime
Definition: Calorimetry_module.cc:124
calo::Calorimetry::fetime
std::vector< double > fetime
Definition: Calorimetry_module.cc:126
recob::tracking::convertCollToPoint
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
Definition: TrackingTypes.h:68
geo::TPCID
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
tca::dEdx
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
Definition: PFPUtils.cxx:2675
calo::Calorimetry::fMIPs
std::vector< double > fMIPs
Definition: Calorimetry_module.cc:127
calo::Calorimetry::GetPitch
void GetPitch(detinfo::DetectorPropertiesData const &det_prop, art::Ptr< recob::Hit > const &hit, std::vector< double > const &trkx, std::vector< double > const &trky, std::vector< double > const &trkz, std::vector< double > const &trkw, std::vector< double > const &trkx0, double *xyz3d, double &pitch, double TickT0)
Definition: Calorimetry_module.cc:675
calo::Calorimetry::caloAlg
CalorimetryAlg caloAlg
Definition: Calorimetry_module.cc:120
util::CreateAssn
bool CreateAssn(art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
Definition: AssociationUtil.h:558
geo::Point_t
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:180
geo::GeometryCore::PlaneWireToChannel
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
Definition: GeometryCore.cxx:1124
art::ProductRetriever::getByLabel
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: ProductRetriever.h:444
dir
TDirectory * dir
Definition: macro.C:5
lar::util::absDiff
constexpr auto absDiff(A const &a, B const &b)
Returns the absolute value of the difference between two values.
Definition: NumericUtils.h:69
calo::Calorimetry::fdQdx
std::vector< double > fdQdx
Definition: Calorimetry_module.cc:128
mf::LogWarning
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Definition: MessageLogger.h:210
recob::Track::Point_t
tracking::Point_t Point_t
Definition: Track.h:52
lar::util::TrackPitchInView
double TrackPitchInView(recob::Track const &track, geo::View_t view, size_t trajectory_point=0U)
Returns the projected length of track on a wire pitch step [cm].
Definition: TrackUtils.cxx:75
MF_LOG_DEBUG
#define MF_LOG_DEBUG(id)
Definition: MessageLogger.h:294
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:127
calo::Calorimetry::fwire
std::vector< int > fwire
Definition: Calorimetry_module.cc:123
art
Definition: MVAAlg.h:12
calo::Calorimetry::fTrackModuleLabel
std::string fTrackModuleLabel
Definition: Calorimetry_module.cc:111
geo::GeometryCore::Nplanes
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
Definition: GeometryCore.h:977
art::detail::Vertex
boost::graph_traits< ModuleGraph >::vertex_descriptor Vertex
Definition: ModuleGraph.h:25
util::kBogusD
constexpr double kBogusD
obviously bogus double value
Definition: PhysicalConstants.h:71
calo::Calorimetry::fstime
std::vector< double > fstime
Definition: Calorimetry_module.cc:125
art::fill_ptr_vector
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
geo::TPCID::TPC
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:399
calo::Calorimetry::fpitch
std::vector< float > fpitch
Definition: Calorimetry_module.cc:131
e
Float_t e
Definition: plot.C:35
detinfo::sampling_rate
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
Definition: DetectorClocksData.h:575
trkf::Plane
recob::tracking::Plane Plane
Definition: TrackState.h:17
X
Float_t X
Definition: plot.C:37
w
Float_t w
Definition: plot.C:20
geo::GeometryCore::WirePitch
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
Definition: GeometryCore.cxx:763
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void calo::Calorimetry::ReadCaloTree ( )
private
void art::Modifier::registerProducts ( ProductDescriptions &  productsToRegister)
inherited

Definition at line 16 of file Modifier.cc.

References art::ModuleBase::moduleDescription(), and art::ProductRegistryHelper::registerProducts().

17  {
18  ProductRegistryHelper::registerProducts(productsToRegister,
19  moduleDescription());
20  }
art::ProductRegistryHelper::registerProducts
void registerProducts(ProductDescriptions &productsToRegister, ModuleDescription const &md)
Definition: ProductRegistryHelper.cc:24
art::ModuleBase::moduleDescription
ModuleDescription const & moduleDescription() const
Definition: ModuleBase.cc:13
void art::ModuleBase::setModuleDescription ( ModuleDescription const &  md)
inherited

Definition at line 31 of file ModuleBase.cc.

References art::ModuleBase::md_.

32  {
33  md_ = md;
34  }
art::ModuleBase::md_
std::optional< ModuleDescription > md_
Definition: ModuleBase.h:55
void art::ModuleBase::sortConsumables ( std::string const &  current_process_name)
inherited

Definition at line 49 of file ModuleBase.cc.

References art::ModuleBase::collector_, and art::ConsumesCollector::sortConsumables().

50  {
51  // Now that we know we have seen all the consumes declarations,
52  // sort the results for fast lookup later.
53  collector_.sortConsumables(current_process_name);
54  }
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:56
art::ConsumesCollector::sortConsumables
void sortConsumables(std::string const &current_process_name)
Definition: ConsumesCollector.cc:19

Member Data Documentation

CalorimetryAlg calo::Calorimetry::caloAlg
private

Definition at line 120 of file Calorimetry_module.cc.

Referenced by Calorimetry(), and produce().

std::vector<double> calo::Calorimetry::fdEdx
private

Definition at line 129 of file Calorimetry_module.cc.

Referenced by produce().

std::vector<double> calo::Calorimetry::fdQdx
private

Definition at line 128 of file Calorimetry_module.cc.

Referenced by produce().

std::vector<double> calo::Calorimetry::fetime
private

Definition at line 126 of file Calorimetry_module.cc.

Referenced by produce().

bool calo::Calorimetry::fFlipTrack_dQdx
private

Definition at line 118 of file Calorimetry_module.cc.

Referenced by Calorimetry(), and produce().

std::vector<size_t> calo::Calorimetry::fHitIndex
private

Definition at line 133 of file Calorimetry_module.cc.

Referenced by produce().

std::vector<double> calo::Calorimetry::fMIPs
private

Definition at line 127 of file Calorimetry_module.cc.

Referenced by produce().

std::optional<double> calo::Calorimetry::fNotOnTrackZcut
private

Exclude trajectory points with z lower than this [cm]

Definition at line 116 of file Calorimetry_module.cc.

Referenced by Calorimetry(), and produce().

int calo::Calorimetry::fnsps
private

Definition at line 122 of file Calorimetry_module.cc.

Referenced by produce().

std::vector<float> calo::Calorimetry::fpitch
private

Definition at line 131 of file Calorimetry_module.cc.

Referenced by produce().

std::vector<double> calo::Calorimetry::fResRng
private

Definition at line 130 of file Calorimetry_module.cc.

Referenced by produce().

bool calo::Calorimetry::fSCE
private

Definition at line 115 of file Calorimetry_module.cc.

Referenced by Calorimetry(), GetPitch(), and produce().

std::string calo::Calorimetry::fSpacePointModuleLabel
private

Definition at line 112 of file Calorimetry_module.cc.

Referenced by Calorimetry(), and produce().

std::vector<double> calo::Calorimetry::fstime
private

Definition at line 125 of file Calorimetry_module.cc.

Referenced by produce().

std::string calo::Calorimetry::fT0ModuleLabel
private

Definition at line 113 of file Calorimetry_module.cc.

Referenced by Calorimetry(), and produce().

std::vector<double> calo::Calorimetry::ftime
private

Definition at line 124 of file Calorimetry_module.cc.

Referenced by produce().

std::string calo::Calorimetry::fTrackModuleLabel
private

Definition at line 111 of file Calorimetry_module.cc.

Referenced by Calorimetry(), and produce().

bool calo::Calorimetry::fUseArea
private

Definition at line 114 of file Calorimetry_module.cc.

Referenced by Calorimetry(), and produce().

std::vector<int> calo::Calorimetry::fwire
private

Definition at line 123 of file Calorimetry_module.cc.

Referenced by produce().

std::vector<TVector3> calo::Calorimetry::fXYZ
private

Definition at line 132 of file Calorimetry_module.cc.

Referenced by produce().

The documentation for this class was generated from the following file: