All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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:

Public Member Functions

 Calorimetry (fhicl::ParameterSet const &pset)
 

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.

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 }
calo::Calorimetry::fNotOnTrackZcut
std::optional< double > fNotOnTrackZcut
Definition: Calorimetry_module.cc:116
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
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
701  art::ServiceHandle<geo::Geometry const> geom;
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 }
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:164
hit
process_name hit
Definition: cheaterreco.fcl:51
make_calib_df.t0
int t0
Definition: make_calib_df.py:42
icarus::distance
double distance(geo::Point_t const &point, CathodeDesc_t const &cathode)
Returns the distance of a point from the cathode.
Definition: CathodeCrossingUtils.cxx:35
dropbox.dir
tuple dir
Definition: dropbox.py:28
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:184
void calo::Calorimetry::produce ( art::Event &  evt)
overrideprivate

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 =
163  art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clock_data);
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
172  art::ServiceHandle<geo::Geometry const> geom;
173 
174  // channel quality
175  lariov::ChannelStatusProvider const& channelStatus =
176  art::ServiceHandle<lariov::ChannelStatusService const>()->GetProvider();
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(
184  new art::Assns<recob::Track, anab::Calorimetry>);
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,
258  util::kBogusD,
259  fpitch,
260  recob::tracking::convertCollToPoint(vXYZ),
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,
525  util::kBogusD,
526  fpitch,
527  recob::tracking::convertCollToPoint(vXYZ),
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,
673  recob::tracking::convertCollToPoint(vXYZ),
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 }
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:164
pmtana::std
double std(const std::vector< short > &wf, const double ped_mean, size_t start, size_t nsample)
Definition: UtilFunc.cxx:42
if
then if[["$THISISATEST"==1]]
Definition: neoSmazza.sh:95
calo::Calorimetry::fdEdx
std::vector< double > fdEdx
Definition: Calorimetry_module.cc:129
calo::Calorimetry::fNotOnTrackZcut
std::optional< double > fNotOnTrackZcut
Definition: Calorimetry_module.cc:116
Plane
BEGIN_PROLOG true icarus_rawdigitfilter FilterTools FilterPlane1 Plane
Definition: rawdigitfilter_icarus.fcl:8
x
process_name opflash particleana ie x
Definition: run_opflash_electron.fcl:90
lariov::ChannelStatusProvider::IsBad
virtual bool IsBad(raw::ChannelID_t channel) const =0
Returns whether the specified channel is bad in the current run.
geo::PlaneID
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
geo::CryostatID::isValid
bool isValid
Whether this ID points to a valid element.
Definition: geo_types.h:211
geo::CryostatID::Cryostat
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:212
icarus::ns::util::size
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
calo::Calorimetry::fXYZ
std::vector< TVector3 > fXYZ
Definition: Calorimetry_module.cc:132
X
then echo echo For and will not be changed by echo further linking echo echo B echo The symbol is in the uninitialized data multiple common symbols may appear with the echo same name If the symbol is defined the common echo symbols are treated as undefined references For more echo details on common see the discussion of warn common echo in *Note Linker see the discussion of warn common echo in *Note Linker such as a global int variable echo as opposed to a large global array echo echo I echo The symbol is an indirect reference to another symbol This echo is a GNU extension to the a out object file format which is echo rarely used echo echo N echo The symbol is a debugging symbol echo echo R echo The symbol is in a read only data section echo echo S echo The symbol is in an uninitialized data section for small echo objects echo echo T echo The symbol is in the the normal defined echo symbol is used with no error When a weak undefined symbol echo is linked and the symbol is not the value of the echo weak symbol becomes zero with no error echo echo W echo The symbol is a weak symbol that has not been specifically echo tagged as a weak object symbol When a weak defined symbol echo is linked with a normal defined the normal defined echo symbol is used with no error When a weak undefined symbol echo is linked and the symbol is not the value of the echo weak symbol becomes zero with no error echo echo echo The symbol is a stabs symbol in an a out object file In echo this the next values printed are the stabs other echo the stabs desc and the stab type Stabs symbols are echo used to hold debugging information For more echo see *Note or object file format specific echo echo For Mac OS X
Definition: find_global_symbol.sh:107
TPC
BEGIN_PROLOG TPC
Definition: sbnd_flashalgo.fcl:14
calo::Calorimetry::fResRng
std::vector< double > fResRng
Definition: Calorimetry_module.cc:130
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
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:42
geo::fhicl::WireID
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Definition: geo_types_fhicl.h:407
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:55
lariov::ChannelStatusProvider
Class providing information about the quality of channels.
Definition: ChannelStatusProvider.h:49
geo::TPCID
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
tca::dEdx
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
Definition: PFPUtils.cxx:2687
calo::Calorimetry::fMIPs
std::vector< double > fMIPs
Definition: Calorimetry_module.cc:127
gen_crt_frags.w
tuple w
Definition: gen_crt_frags.py:1416
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:686
util::absDiff
constexpr auto absDiff(A const &a, B const &b)
Returns the absolute value of the difference between two values.
Definition: NumericUtils.h:43
calo::Calorimetry::caloAlg
CalorimetryAlg caloAlg
Definition: Calorimetry_module.cc:120
dropbox.dir
tuple dir
Definition: dropbox.py:28
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: lardata/lardata/Utilities/AssociationUtil.h:586
calo::Calorimetry::fdQdx
std::vector< double > fdQdx
Definition: Calorimetry_module.cc:128
recob::Track::Point_t
tracking::Point_t Point_t
Definition: bj/lardataobj/RecoBase/Track.h:53
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:79
calo::CalorimetryAlg::dEdx_AREA
double dEdx_AREA(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
Definition: CalorimetryAlg.cxx:131
e
do i e
Definition: updateVersionList.sh:68
calo::Calorimetry::fwire
std::vector< int > fwire
Definition: Calorimetry_module.cc:123
calo::Calorimetry::fTrackModuleLabel
std::string fTrackModuleLabel
Definition: Calorimetry_module.cc:111
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
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
geo::TPCID::TPC
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
calo::Calorimetry::fpitch
std::vector< float > fpitch
Definition: Calorimetry_module.cc:131
detinfo::sampling_rate
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
Definition: DetectorClocksData.h:643
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:184
T0
process_name opdaq physics producers generator physics producers generator physics producers generator physics producers generator physics producers generator physics producers generator physics producers generator physics producers generator T0
Definition: gen_protons.fcl:45
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:
Generated by   doxygen 1.8.5