Public 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 Member Functions

 Calorimetry (fhicl::ParameterSet const &pset)
 
- Public Member Functions inherited from art::EDProducer
 EDProducer (fhicl::ParameterSet const &pset)
 
template<typename Config >
 EDProducer (Table< Config > const &config)
 
std::string workerType () const
 
- Public Member Functions inherited from art::detail::Producer
virtual ~Producer () noexcept
 
 Producer (fhicl::ParameterSet const &)
 
 Producer (Producer const &)=delete
 
 Producer (Producer &&)=delete
 
Produceroperator= (Producer const &)=delete
 
Produceroperator= (Producer &&)=delete
 
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)
 
- Public Member Functions inherited from art::Modifier
 ~Modifier () noexcept
 
 Modifier ()
 
 Modifier (Modifier const &)=delete
 
 Modifier (Modifier &&)=delete
 
Modifieroperator= (Modifier const &)=delete
 
Modifieroperator= (Modifier &&)=delete
 
- Public Member Functions inherited from art::ModuleBase
virtual ~ModuleBase () noexcept
 
 ModuleBase ()
 
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)
 
template<typename T , BranchType BT>
ViewToken< T > consumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< T > mayConsumeView (InputTag const &tag)
 

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
 

Additional Inherited Members

- Public Types inherited from art::EDProducer
using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
- Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using Table = Modifier::Table< UserConfig, KeysToIgnore >
 
- Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >
 
- Static Public Member Functions inherited from art::EDProducer
static void commitEvent (EventPrincipal &ep, Event &e)
 
- Protected Member Functions inherited from art::ModuleBase
ConsumesCollectorconsumesCollector ()
 
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 ()
 

Detailed Description

Estimates the energy deposited by reconstructed tracks.

Output

Configuration

Note
This documentation is grossly incomplete.

Definition at line 88 of file Calorimetry_module.cc.

Constructor & Destructor Documentation

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

Definition at line 140 of file Calorimetry_module.cc.

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 }
std::optional< double > fNotOnTrackZcut
std::string string
Definition: nybbler.cc:12
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
std::string fSpacePointModuleLabel
CalorimetryAlg caloAlg
std::string fTrackModuleLabel

Member Function Documentation

bool calo::Calorimetry::BeginsOnBoundary ( art::Ptr< recob::Track lar_track)
private
bool calo::Calorimetry::EndsOnBoundary ( art::Ptr< recob::Track lar_track)
private
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 686 of file Calorimetry_module.cc.

696 {
697  // Get 3d coordinates and track pitch for each hit
698  // Find 5 nearest space points and determine xyz and curvature->track pitch
699 
700  // Get services
702  auto const* sce = lar::providerFrom<spacecharge::SpaceChargeService>();
703 
704  //save distance to each spacepoint sorted by distance
705  std::map<double, size_t> sptmap;
706  //save the sign of distance
707  std::map<size_t, int> sptsignmap;
708 
709  double wire_pitch = geom->WirePitch(0);
710 
711  double t0 = hit->PeakTime() - TickT0;
712  double x0 =
713  det_prop.ConvertTicksToX(t0, hit->WireID().Plane, hit->WireID().TPC, hit->WireID().Cryostat);
714  double w0 = hit->WireID().Wire;
715 
716  for (size_t i = 0; i < trkx.size(); ++i) {
717  double distance = cet::sum_of_squares((trkw[i] - w0) * wire_pitch, trkx0[i] - x0);
718  if (distance > 0) distance = sqrt(distance);
719  sptmap.insert(std::pair<double, size_t>(distance, i));
720  if (w0 - trkw[i] > 0)
721  sptsignmap.emplace(i, 1);
722  else
723  sptsignmap.emplace(i, -1);
724  }
725 
726  //x,y,z vs distance
727  std::vector<double> vx;
728  std::vector<double> vy;
729  std::vector<double> vz;
730  std::vector<double> vs;
731 
732  double kx = 0, ky = 0, kz = 0;
733 
734  int np = 0;
735  for (auto isp = sptmap.begin(); isp != sptmap.end(); isp++) {
736  double xyz[3];
737  xyz[0] = trkx[isp->second];
738  xyz[1] = trky[isp->second];
739  xyz[2] = trkz[isp->second];
740 
741  double distancesign = sptsignmap[isp->second];
742  if (np == 0 && isp->first > 30) { // hit not on track
743  xyz3d[0] = std::numeric_limits<double>::lowest();
744  xyz3d[1] = std::numeric_limits<double>::lowest();
745  xyz3d[2] = std::numeric_limits<double>::lowest();
746  pitch = -1;
747  return;
748  }
749  if (np < 5) {
750  vx.push_back(xyz[0]);
751  vy.push_back(xyz[1]);
752  vz.push_back(xyz[2]);
753  vs.push_back(isp->first * distancesign);
754  }
755  else {
756  break;
757  }
758  np++;
759  }
760  if (np >= 2) { // at least two points
761  TGraph* xs = new TGraph(np, &vs[0], &vx[0]);
762  try {
763  if (np > 2) { xs->Fit("pol2", "Q"); }
764  else {
765  xs->Fit("pol1", "Q");
766  }
767  TF1* pol = 0;
768  if (np > 2)
769  pol = (TF1*)xs->GetFunction("pol2");
770  else
771  pol = (TF1*)xs->GetFunction("pol1");
772  xyz3d[0] = pol->Eval(0);
773  kx = pol->GetParameter(1);
774  }
775  catch (...) {
776  mf::LogWarning("Calorimetry::GetPitch") << "Fitter failed";
777  xyz3d[0] = vx[0];
778  }
779  delete xs;
780  TGraph* ys = new TGraph(np, &vs[0], &vy[0]);
781  try {
782  if (np > 2) { ys->Fit("pol2", "Q"); }
783  else {
784  ys->Fit("pol1", "Q");
785  }
786  TF1* pol = 0;
787  if (np > 2)
788  pol = (TF1*)ys->GetFunction("pol2");
789  else
790  pol = (TF1*)ys->GetFunction("pol1");
791  xyz3d[1] = pol->Eval(0);
792  ky = pol->GetParameter(1);
793  }
794  catch (...) {
795  mf::LogWarning("Calorimetry::GetPitch") << "Fitter failed";
796  xyz3d[1] = vy[0];
797  }
798  delete ys;
799  TGraph* zs = new TGraph(np, &vs[0], &vz[0]);
800  try {
801  if (np > 2) { zs->Fit("pol2", "Q"); }
802  else {
803  zs->Fit("pol1", "Q");
804  }
805  TF1* pol = 0;
806  if (np > 2)
807  pol = (TF1*)zs->GetFunction("pol2");
808  else
809  pol = (TF1*)zs->GetFunction("pol1");
810  xyz3d[2] = pol->Eval(0);
811  kz = pol->GetParameter(1);
812  }
813  catch (...) {
814  mf::LogWarning("Calorimetry::GetPitch") << "Fitter failed";
815  xyz3d[2] = vz[0];
816  }
817  delete zs;
818  }
819  else if (np) {
820  xyz3d[0] = vx[0];
821  xyz3d[1] = vy[0];
822  xyz3d[2] = vz[0];
823  }
824  else {
825  xyz3d[0] = std::numeric_limits<double>::lowest();
826  xyz3d[1] = std::numeric_limits<double>::lowest();
827  xyz3d[2] = std::numeric_limits<double>::lowest();
828  pitch = -1;
829  return;
830  }
831  pitch = -1;
832  if (kx * kx + ky * ky + kz * kz) {
833  double tot = sqrt(kx * kx + ky * ky + kz * kz);
834  kx /= tot;
835  ky /= tot;
836  kz /= tot;
837  //get pitch
838  double wirePitch =
839  geom->WirePitch(hit->WireID().Plane, hit->WireID().TPC, hit->WireID().Cryostat);
840  double angleToVert = geom->Plane(hit->WireID().Plane, hit->WireID().TPC, hit->WireID().Cryostat)
841  .Wire(0)
842  .ThetaZ(false) -
843  0.5 * TMath::Pi();
844  double cosgamma = TMath::Abs(TMath::Sin(angleToVert) * ky + TMath::Cos(angleToVert) * kz);
845  if (cosgamma > 0) pitch = wirePitch / cosgamma;
846 
847  //Correct for SCE
848  geo::Vector_t posOffsets = {0., 0., 0.};
849  geo::Vector_t dirOffsets = {0., 0., 0.};
850  if (sce->EnableCalSpatialSCE() && fSCE)
851  posOffsets =
852  sce->GetCalPosOffsets(geo::Point_t{xyz3d[0], xyz3d[1], xyz3d[2]}, hit->WireID().TPC);
853  if (sce->EnableCalSpatialSCE() && fSCE)
854  dirOffsets = sce->GetCalPosOffsets(
855  geo::Point_t{xyz3d[0] + pitch * kx, xyz3d[1] + pitch * ky, xyz3d[2] + pitch * kz},
856  hit->WireID().TPC);
857 
858  xyz3d[0] = xyz3d[0] - posOffsets.X();
859  xyz3d[1] = xyz3d[1] + posOffsets.Y();
860  xyz3d[2] = xyz3d[2] + posOffsets.Z();
861 
862  // Correct pitch for SCE
863  TVector3 dir = {pitch * kx - dirOffsets.X() + posOffsets.X(),
864  pitch * ky + dirOffsets.Y() - posOffsets.Y(),
865  pitch * kz + dirOffsets.Z() - posOffsets.Z()};
866  pitch = dir.Mag();
867  }
868 }
code to link reconstructed objects back to the MC truth information
PlaneGeo const & Plane(unsigned int const p, unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified wire.
geo::WireID WireID() const
Definition: Hit.h:233
constexpr T sum_of_squares(T x, T y)
Definition: pow.h:139
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:212
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:164
string dir
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
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:184
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:218
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
static constexpr double zs
Definition: Units.h:102
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
static constexpr double ys
Definition: Units.h:103
void calo::Calorimetry::produce ( art::Event evt)
overrideprivatevirtual

Implements art::EDProducer.

Definition at line 159 of file Calorimetry_module.cc.

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

Member Data Documentation

CalorimetryAlg calo::Calorimetry::caloAlg
private

Definition at line 120 of file Calorimetry_module.cc.

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

Definition at line 129 of file Calorimetry_module.cc.

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

Definition at line 128 of file Calorimetry_module.cc.

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

Definition at line 126 of file Calorimetry_module.cc.

bool calo::Calorimetry::fFlipTrack_dQdx
private

Definition at line 118 of file Calorimetry_module.cc.

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

Definition at line 133 of file Calorimetry_module.cc.

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

Definition at line 127 of file Calorimetry_module.cc.

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.

int calo::Calorimetry::fnsps
private

Definition at line 122 of file Calorimetry_module.cc.

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

Definition at line 131 of file Calorimetry_module.cc.

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

Definition at line 130 of file Calorimetry_module.cc.

bool calo::Calorimetry::fSCE
private

Definition at line 115 of file Calorimetry_module.cc.

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

Definition at line 112 of file Calorimetry_module.cc.

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

Definition at line 125 of file Calorimetry_module.cc.

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

Definition at line 113 of file Calorimetry_module.cc.

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

Definition at line 124 of file Calorimetry_module.cc.

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

Definition at line 111 of file Calorimetry_module.cc.

bool calo::Calorimetry::fUseArea
private

Definition at line 114 of file Calorimetry_module.cc.

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

Definition at line 123 of file Calorimetry_module.cc.

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

Definition at line 132 of file Calorimetry_module.cc.


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