37 , fWireReadoutGeom(wireReadoutGeom)
38 , fDelay(config.Delay())
39 , fG4ModuleLabel(config.G4ModuleLabel())
40 , fG4ModuleLabels(config.G4ModuleLabels())
41 , fOpHitLabel(config.OpHitLabel())
42 , fOpFlashLabel(config.OpFlashLabel())
43 , fMinOpHitEnergyFraction(config.MinOpHitEnergyFraction())
51 , fWireReadoutGeom(wireReadoutGeom)
52 , fDelay(pSet.
get<double>(
"Delay"))
54 , fG4ModuleLabels(pSet.
get<
std::
vector<
art::InputTag>>(
"G4ModuleLabels", {}))
56 , fOpFlashLabel(pSet.get<
art::InputTag>(
"OpFlashLabel",
"opflash"))
57 , fMinOpHitEnergyFraction(pSet.get<
double>(
"MinimumOpHitEnergyFraction", 0.1))
69 priv_OpDetBTRs.clear();
70 priv_OpFlashToOpHits.clear();
76 return !priv_OpDetBTRs.empty();
82 return !priv_OpFlashToOpHits.empty();
88 return priv_OpDetBTRs;
94 std::vector<const sim::SDP*> sdp_Ps;
95 for (
size_t odet = 0; odet < priv_OpDetBTRs.size(); ++odet) {
96 const auto& pdTimeSDPmap = priv_OpDetBTRs[odet]->timePDclockSDPsMap();
97 for (
auto mapitr = pdTimeSDPmap.begin(); mapitr != pdTimeSDPmap.end(); mapitr++) {
98 std::vector<sim::SDP>
const& sdpvec = (*mapitr).second;
99 for (
size_t iv = 0; iv < sdpvec.size(); ++iv) {
100 if (
abs(sdpvec[iv].trackID) == id) sdp_Ps.push_back(&(sdpvec[iv]));
111 for (
size_t sc = 0;
sc < priv_OpDetBTRs.size(); ++
sc) {
114 if (priv_OpDetBTRs.at(
sc)->OpDetNum() == opDetNum) opDet = priv_OpDetBTRs.at(
sc);
117 throw cet::exception(
"PhotonBackTracker2") <<
"No sim:: OpDetBacktrackerRecord corresponding " 118 <<
"to opDetNum: " << opDetNum <<
"\n";
125 double const opHit_start_time,
126 double const opHit_end_time)
const 128 std::vector<sim::TrackSDP> tSDPs;
132 std::vector<sim::SDP> simSDPs =
134 for (
size_t e = 0;
e < simSDPs.size(); ++
e)
136 if (totalE < 1.
e-5) totalE = 1.;
137 for (
size_t e = 0;
e < simSDPs.size(); ++
e) {
142 info.
energy = simSDPs[
e].energy;
143 tSDPs.push_back(info);
158 auto OpDetNum = fWireReadoutGeom->OpDetFromOpChannel(opHit_P->
OpChannel());
159 const double pTime = opHit_P->
PeakTime();
160 const double pWidth = opHit_P->
Width();
161 const double start = (pTime - pWidth) * 1000 - fDelay;
162 const double end = (pTime + pWidth) * 1000 - fDelay;
163 return OpDetToTrackSDPs(OpDetNum, start, end);
169 auto OpDetNum = fWireReadoutGeom->OpDetFromOpChannel(opHit.
OpChannel());
170 const double pTime = opHit.
PeakTime();
171 const double pWidth = opHit.
Width();
172 const double start = (pTime - pWidth) * 1000 - fDelay;
173 const double end = (pTime + pWidth) * 1000 - fDelay;
174 return OpDetToTrackSDPs(OpDetNum, start, end);
180 std::vector<int> retVec;
181 for (
auto const trackSDP : OpHitToTrackSDPs(opHit)) {
182 retVec.push_back(trackSDP.trackID);
190 return OpHitToTrackIds(*opHit);
196 std::vector<int> retVec;
197 for (
auto const trackSDP : OpHitToEveTrackSDPs(opHit)) {
198 retVec.push_back(trackSDP.trackID);
206 return OpHitToEveTrackIds(*opHit_P);
213 return OpHitToEveTrackSDPs(*opHit_P);
219 std::vector<sim::TrackSDP> trackSDPs = OpHitToTrackSDPs(opHit);
223 std::map<int, float> eveIDtoEfrac;
226 for (
size_t t = 0; t < trackSDPs.size(); ++t) {
227 eveIDtoEfrac[(fPartInv->ParticleList()).EveId(trackSDPs[t].trackID)] += trackSDPs[t].energy;
228 totalE += trackSDPs[t].energy;
232 std::vector<sim::TrackSDP> eveSDPs;
233 eveSDPs.reserve(eveIDtoEfrac.size());
234 for (
auto itr = eveIDtoEfrac.begin(); itr != eveIDtoEfrac.end(); itr++) {
238 temp.
energy = (*itr).second;
239 eveSDPs.push_back(std::move(temp));
254 std::vector<int> tkidFake(1, tkId);
255 return TrackIdsToOpHits_Ps(tkidFake, hitsIn).at(0);
260 std::vector<int>
const& tkIds,
263 std::vector<std::pair<int, art::Ptr<recob::OpHit>>> opHitList;
264 for (
auto itr = hitsIn.begin(); itr != hitsIn.end(); ++itr) {
266 auto OpDetNum = fWireReadoutGeom->OpDetFromOpChannel(opHit->
OpChannel());
267 const double pTime = opHit->
PeakTime(), pWidth = opHit->
Width();
268 const double start = (pTime - pWidth) * 1000.0 - fDelay,
269 end = (pTime + pWidth) * 1000.0 - fDelay;
270 std::vector<sim::TrackSDP> tids = OpDetToTrackSDPs(OpDetNum, start,
end);
271 for (
auto itid = tids.begin(); itid != tids.end(); ++itid) {
272 for (
auto itkid = tkIds.begin(); itkid != tkIds.end(); ++itkid) {
273 if (itid->trackID == *itkid and itid->energyFrac > fMinOpHitEnergyFraction)
274 opHitList.push_back(std::make_pair(*itkid, opHit));
279 std::vector<std::vector<art::Ptr<recob::OpHit>>> truOpHits;
281 std::vector<art::Ptr<recob::OpHit>> tmpOpHits;
282 for (
auto itkid = tkIds.begin(); itkid != tkIds.end(); ++itkid) {
284 for (
auto itr = opHitList.begin(); itr != opHitList.end(); ++itr) {
285 if (*itkid == (*itr).first) tmpOpHits.push_back((*itr).second);
287 truOpHits.push_back(tmpOpHits);
295 std::vector<const sim::SDP*> retVec;
296 double fPeakTime = opHit.
PeakTime();
297 double fWidth = opHit.
Width();
299 ((fPeakTime - fWidth) * 1000.0) - fDelay;
301 if (start_time > end_time) {
throw; }
303 auto const& timeSDPMap = FindOpDetBTR(fWireReadoutGeom->OpDetFromOpChannel(opHit.
OpChannel()))
304 ->timePDclockSDPsMap();
306 std::vector<const std::pair<double, std::vector<sim::SDP>>*> timePDclockSDPMap_SortedPointers;
307 for (
auto& pair : timeSDPMap) {
308 timePDclockSDPMap_SortedPointers.push_back(&pair);
310 auto pairSort = [](
auto& a,
auto& b) {
return a->first < b->first; };
311 if (!std::is_sorted(timePDclockSDPMap_SortedPointers.begin(),
312 timePDclockSDPMap_SortedPointers.end(),
315 timePDclockSDPMap_SortedPointers.begin(), timePDclockSDPMap_SortedPointers.end(), pairSort);
319 std::vector<sim::SDP> dummyVec;
320 std::pair<double, std::vector<sim::SDP>> start_timePair = std::make_pair(start_time, dummyVec);
321 std::pair<double, std::vector<sim::SDP>> end_timePair = std::make_pair(end_time, dummyVec);
322 auto start_timePair_P = &start_timePair;
323 auto end_timePair_P = &end_timePair;
325 auto mapFirst = std::lower_bound(timePDclockSDPMap_SortedPointers.begin(),
326 timePDclockSDPMap_SortedPointers.end(),
331 std::upper_bound(mapFirst, timePDclockSDPMap_SortedPointers.end(), end_timePair_P, pairSort);
333 for (
auto& mapitr = mapFirst; mapitr != mapLast; ++mapitr)
334 for (
auto& sdp : (*mapitr)->second)
335 retVec.push_back(&sdp);
344 return OpHitToSimSDPs_Ps(*opHit_P);
350 std::vector<double> xyz(3, -999.);
356 for (
auto const& sdp : sdps) {
357 double weight = sdp.numPhotons;
367 <<
"No sim::SDPs providing non-zero number of photons" 368 <<
" can't determine originating location from truth\n";
377 std::vector<const sim::SDP*>
const& sdps_Ps)
const&
379 std::vector<double> xyz(3, -999.);
385 for (
const sim::SDP* sdp_P : sdps_Ps) {
387 double weight = sdp.numPhotons;
397 <<
"No sim::SDPs providing non-zero number of photons" 398 <<
" can't determine originating location from truth\n";
408 return SimSDPsToXYZ(OpHitToSimSDPs_Ps(opHit));
414 return SimSDPsToXYZ(OpHitToSimSDPs_Ps(*opHit));
421 std::vector<const sim::SDP*> sdps_Ps;
422 for (
auto opHit_P : opHits_Ps) {
423 std::vector<const sim::SDP*> to_add_sdps_Ps = OpHitToSimSDPs_Ps(opHit_P);
424 sdps_Ps.insert(sdps_Ps.end(), to_add_sdps_Ps.begin(), to_add_sdps_Ps.end());
433 const std::vector<const sim::SDP*> SDPs_Ps = OpHitsToSimSDPs_Ps(opHits_Ps);
434 return SimSDPsToXYZ(SDPs_Ps);
441 const std::vector<int> ids = OpHitToEveTrackIds(opHit_P);
442 std::unordered_set<const sim::SDP*> sdps;
443 for (
auto const&
id : ids) {
444 std::vector<const sim::SDP*> tmp_sdps = TrackIdToSimSDPs_Ps(
id);
445 for (
const sim::SDP* tmp_sdp : tmp_sdps) {
446 sdps.insert(tmp_sdp);
456 const std::vector<int> ids = OpHitToEveTrackIds(opHit);
457 std::unordered_set<const sim::SDP*> sdps;
458 for (
auto const&
id : ids) {
459 std::vector<const sim::SDP*> tmp_sdps = TrackIdToSimSDPs_Ps(
id);
460 for (
const sim::SDP* tmp_sdp : tmp_sdps) {
461 sdps.insert(tmp_sdp);
470 return fPartInv->GetSetOfEveIds();
476 return fPartInv->GetSetOfTrackIds();
483 std::set<int> eveIds;
484 for (
auto const& opHit_P : opHits_Ps) {
485 const std::vector<sim::TrackSDP> sdps = OpHitToEveTrackSDPs(opHit_P);
486 for (
auto const& sdp : sdps) {
487 eveIds.insert(sdp.trackID);
496 std::set<int> eveIds;
497 for (
auto const& opHit : opHits) {
498 const std::vector<sim::TrackSDP> sdps = OpHitToEveTrackSDPs(opHit);
499 for (
auto const& sdp : sdps) {
500 eveIds.insert(sdp.trackID);
511 for (
auto const& opHit : opHits) {
512 for (
auto const& sdp : OpHitToTrackSDPs(opHit)) {
513 tids.insert(sdp.trackID);
523 for (
auto const& opHit : opHits) {
524 for (
auto const& sdp : OpHitToTrackSDPs(opHit)) {
525 tids.insert(sdp.trackID);
537 float total = 1. * opHits.size();
539 for (
size_t h = 0; h < opHits.size(); ++h) {
541 std::vector<sim::TrackSDP> opHitTrackSDPs = OpHitToTrackSDPs(opHit);
544 for (
size_t e = 0;
e < opHitTrackSDPs.size(); ++
e) {
545 if (tkIds.find(opHitTrackSDPs[
e].trackID) != tkIds.end()) {
552 if (total > 0) purity = desired / total;
558 std::set<int>
const& tkIds,
567 for (
size_t h = 0; h < opHits.size(); ++h) {
569 std::vector<sim::TrackSDP> opHitTrackIDs = OpHitToTrackSDPs(opHit);
570 total += opHit->
Area();
573 for (
size_t e = 0;
e < opHitTrackIDs.size(); ++
e) {
574 if (tkIds.find(opHitTrackIDs[
e].trackID) != tkIds.end()) {
575 desired += opHit->
Area();
581 if (total > 0) purity = desired / total;
587 std::set<int>
const& tkIds,
593 for (
size_t h = 0; h < opHits.size(); ++h) {
595 std::vector<sim::TrackSDP> opHitTrackSDPs = OpHitToTrackSDPs(opHit);
598 for (
size_t e = 0;
e < opHitTrackSDPs.size(); ++
e) {
599 if (tkIds.find(opHitTrackSDPs[
e].trackID) != tkIds.end() &&
600 opHitTrackSDPs[
e].energyFrac >= fMinOpHitEnergyFraction) {
607 for (
size_t h = 0; h < opHitsIn.size(); ++h) {
609 std::vector<sim::TrackSDP> opHitTrackSDPs = OpHitToTrackSDPs(opHit);
610 for (
size_t e = 0;
e < opHitTrackSDPs.size(); ++
e) {
614 if (tkIds.find(opHitTrackSDPs[
e].trackID) != tkIds.end() &&
615 opHitTrackSDPs[
e].energyFrac >= fMinOpHitEnergyFraction) {
621 double efficiency = 0.;
622 if (total > 0.) efficiency = desired / total;
628 std::set<int>
const& tkIds,
637 for (
size_t h = 0; h < opHits.size(); ++h) {
640 std::vector<sim::TrackSDP> opHitTrackIDs = OpHitToTrackSDPs(opHit);
645 for (
size_t e = 0;
e < opHitTrackIDs.size(); ++
e) {
646 if (tkIds.find(opHitTrackIDs[
e].trackID) != tkIds.end() &&
647 opHitTrackIDs[
e].energyFrac >= fMinOpHitEnergyFraction) {
648 desired += opHit->
Area();
653 for (
size_t h = 0; h < opHitsIn.size(); ++h) {
657 std::vector<sim::TrackSDP> opHitTrackIDs = OpHitToTrackSDPs(opHit);
658 for (
size_t e = 0;
e < opHitTrackIDs.size(); ++
e) {
662 if (tkIds.find(opHitTrackIDs[
e].trackID) != tkIds.end() &&
663 opHitTrackIDs[
e].energyFrac >= fMinOpHitEnergyFraction) {
664 total += opHit->
Area();
669 double efficiency = 0.;
670 if (total > 0.) efficiency = desired / total;
680 return priv_OpFlashToOpHits.at(flash_P);
686 return OpHitsToXYZ(OpFlashToOpHits_Ps(flash_P));
692 std::vector<art::Ptr<recob::OpHit>> opHits_Ps = OpFlashToOpHits_Ps(flash_P);
694 for (
auto& opHit_P : opHits_Ps) {
695 for (
const int&
id : OpHitToTrackIds(opHit_P)) {
std::vector< double > OpHitsToXYZ(std::vector< art::Ptr< recob::OpHit >> const &opHits_Ps) const
double OpHitLightCollectionPurity(std::set< int > const &tkIds, std::vector< art::Ptr< recob::OpHit >> const &opHits)
std::vector< double > OpHitToXYZ(art::Ptr< recob::OpHit > const &hit)
std::vector< double > OpFlashToXYZ(art::Ptr< recob::OpFlash > &flash_P) const
std::vector< sim::TrackSDP > OpDetToTrackSDPs(int OpDetNum, double opHit_start_time, double opHit_end_time) const
std::vector< const sim::SDP * > TrackIdToSimSDPs_Ps(int id) const
art::Ptr< sim::OpDetBacktrackerRecord > FindOpDetBTR(int opDetNum) const
std::vector< double > SimSDPsToXYZ(std::vector< sim::SDP > const &sdps)
constexpr auto abs(T v)
Returns the absolute value of the argument.
std::vector< int > OpHitToEveTrackIds(recob::OpHit const &opHit)
double OpHitLightCollectionEfficiency(std::set< int > const &tkIds, std::vector< art::Ptr< recob::OpHit >> const &opHits, std::vector< art::Ptr< recob::OpHit >> const &opHitsIn)
std::set< int > OpFlashToTrackIds(art::Ptr< recob::OpFlash > &flash_P) const
std::vector< art::Ptr< sim::OpDetBacktrackerRecord > > const & OpDetBTRs() const
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
std::vector< const sim::SDP * > OpHitToSimSDPs_Ps(recob::OpHit const &opHit) const
int trackID
Geant4 supplied trackID.
std::vector< const sim::SDP * > OpHitsToSimSDPs_Ps(std::vector< art::Ptr< recob::OpHit >> const &opHits_Ps) const
static const int NoParticleId
Interface for a class providing readout channel mapping to geometry.
float energyFrac
fraction of OpHit energy from the particle with this trackID
double timePDclock_t
Type for iTimePDclock tick used in the interface.
std::set< int > GetSetOfTrackIds() const
back track the reconstruction to the simulation
art::InputTag fG4ModuleLabel
Ionization photons from a Geant4 track.
std::vector< sim::TrackSDP > OpHitToTrackSDPs(art::Ptr< recob::OpHit > const &hit)
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
PhotonBackTracker(fhicl::ParameterSet const &pset, art::ActivityRegistry ®)
code to link reconstructed objects back to the MC truth information
bool OpFlashToOpHitsReady() const
std::vector< int > OpHitToTrackIds(recob::OpHit const &opHit) const
double OpHitCollectionEfficiency(std::set< int > trackIDs, std::vector< art::Ptr< recob::OpHit >> const &hits, std::vector< art::Ptr< recob::OpHit >> const &allhits)
std::vector< std::vector< art::Ptr< recob::OpHit > > > TrackIdsToOpHits_Ps(std::vector< int > const &tkIds, std::vector< art::Ptr< recob::OpHit >> const &hitsIn) const
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Header for the ParticleInvenotry Service Provider.
double OpHitCollectionPurity(std::set< int > trackIDs, std::vector< art::Ptr< recob::OpHit >> const &hits)
std::unordered_set< const sim::SDP * > OpHitToEveSimSDPs_Ps(recob::OpHit const &opHit)
std::vector< sim::SDP > TrackIDsAndEnergies(timePDclock_t startTimePDclock, timePDclock_t endTimePDclock) const
Return all the recorded energy deposition within a time interval.
std::set< int > GetSetOfEveIds() const
float energy
energy from the particle with this trackID [MeV]
std::vector< art::Ptr< recob::OpHit > > OpFlashToOpHits_Ps(art::Ptr< recob::OpFlash > &flash_P) const
std::vector< art::Ptr< recob::OpHit > > TrackIdToOpHits_Ps(int tkId, std::vector< art::Ptr< recob::OpHit >> const &hitsIn) const
std::vector< sim::TrackSDP > OpHitToEveTrackSDPs(art::Ptr< recob::OpHit > const &opHit_P) const
Tools and modules for checking out the basics of the Monte Carlo.
Interface to geometry for wire readouts .
cet::coded_exception< error, detail::translate > exception