All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Classes | Public Types | Public Member Functions | Private Types | Private Member Functions | Private Attributes | Static Private Attributes | List of all members
FlashPredict Class Reference

#include <FlashPredict.hh>

Inheritance diagram for FlashPredict:

Classes

struct  BookKeeping
 
struct  ChargeDigest
 
struct  ChargeMetrics
 
struct  Fits
 
struct  FlashMetrics
 
struct  SimpleFlash
 

Public Types

using ChargeDigestMap = std::map< double, ChargeDigest, std::greater< double >>
 
using OpHitIt = std::vector< recob::OpHit >::iterator
 

Public Member Functions

 FlashPredict (fhicl::ParameterSet const &p)
 
 FlashPredict (FlashPredict const &)=delete
 
 FlashPredict (FlashPredict &&)=delete
 
FlashPredictoperator= (FlashPredict const &)=delete
 
FlashPredictoperator= (FlashPredict &&)=delete
 
void produce (art::Event &evt) override
 
void beginJob () override
 
void endJob () override
 

Private Types

using sFM = sbn::SimpleFlashMatch
 
using Charge = sbn::SimpleFlashMatch::Charge
 
using Flash = sbn::SimpleFlashMatch::Flash
 
using Score = sbn::SimpleFlashMatch::Score
 

Private Member Functions

void initTree (void)
 
void loadMetrics (void)
 
double cheatMCT0 (const std::vector< art::Ptr< recob::Hit >> &hits, const std::vector< art::Ptr< simb::MCParticle >> &mcParticles)
 
ChargeMetrics computeChargeMetrics (const flashmatch::QCluster_t &qClusters) const
 
FlashMetrics computeFlashMetrics (const SimpleFlash &simpleFlash) const
 
Score computeScore (const ChargeMetrics &charge, const FlashMetrics &flash, const std::set< unsigned > &tpcWithHits, const int pdgc) const
 
std::tuple< double, double,
double, double > 
hypoFlashX_fits (double flash_rr, double flash_ratio) const
 
std::tuple< double, double,
double, double > 
hypoFlashX_H2 (double flash_rr, double flash_ratio) const
 
std::tuple< double, double > xEstimateAndRMS (double metric_value, const TH2D *metric_h2) const
 
ChargeDigestMap makeChargeDigest (const art::Event &evt, const art::ValidHandle< std::vector< recob::PFParticle >> &pfps_h)
 
void addDaughters (const std::unordered_map< size_t, size_t > &pfpMap, const art::Ptr< recob::PFParticle > &pfp_ptr, const art::ValidHandle< std::vector< recob::PFParticle >> &pfps_h, std::vector< art::Ptr< recob::PFParticle >> &pfp_v) const
 
unsigned trueNus (art::Event &evt) const
 
void updateChargeMetrics (const ChargeMetrics &chargeMetrics)
 
void updateFlashMetrics (const FlashMetrics &flashMetrics)
 
void updateScore (const Score &score)
 
double scoreTerm (const double m, const double n, const double mean, const double spread) const
 
double scoreTerm (const double m, const double mean, const double spread) const
 
bool pfpNeutrinoOnEvent (const art::ValidHandle< std::vector< recob::PFParticle >> &pfps_h) const
 
void copyOpHitsInBeamWindow (std::vector< recob::OpHit > &opHits, const art::Handle< std::vector< recob::OpHit >> &ophit_h) const
 
std::vector< SimpleFlashmakeSimpleFlashes (std::vector< recob::OpHit > &opHits, std::vector< recob::OpHit > &opHitsRght, std::vector< recob::OpHit > &opHitsLeft) const
 
std::vector< SimpleFlashmakeSimpleFlashes (std::vector< recob::OpHit > &opHits) const
 
bool createOpHitsTimeHist (const std::vector< recob::OpHit > &opHits, std::vector< recob::OpHit > &opHitsRght, std::vector< recob::OpHit > &opHitsLeft, std::unique_ptr< TH1D > &opHitsTimeHist, std::unique_ptr< TH1D > &opHitsTimeHistRght, std::unique_ptr< TH1D > &opHitsTimeHistLeft) const
 
unsigned createOpHitsTimeHist (const std::vector< recob::OpHit > &opHits, std::unique_ptr< TH1D > &opHitsTimeHist) const
 
bool findSimpleFlashes (std::vector< SimpleFlash > &simpleFlashes, std::vector< recob::OpHit > &opHits, const unsigned ophsInVolume, std::unique_ptr< TH1D > &opHitsTimeHist) const
 
std::string detectorName (const std::string detName) const
 
bool isPDInCryo (const int pdChannel) const
 
bool isSBNDPDRelevant (const int pdChannel, const std::set< unsigned > &tpcWithHits) const
 
unsigned sbndPDinTPC (const int pdChannel) const
 
unsigned icarusPDinTPC (const int pdChannel) const
 
double wallXWithMaxPE (const OpHitIt opH_beg, const OpHitIt opH_end) const
 
std::list< double > wiresXGl () const
 
double driftDistance () const
 
double flashXGl (const double hypo_x, const double flash_x) const
 
double foldXGl (const double x_gl) const
 
unsigned driftVolume (const double x) const
 
template<typename Stream >
void printBookKeeping (Stream &&out)
 
void updateBookKeeping ()
 
template<typename Stream >
void printMetrics (const std::string metric, const ChargeMetrics &charge, const FlashMetrics &flash, const int pdgc, const std::set< unsigned > &tpcWithHits, const double term, Stream &&out) const
 

Private Attributes

const art::InputTag fPandoraProducer
 
const art::InputTag fSpacePointProducer
 
const art::InputTag fOpHitProducer
 
const art::InputTag fOpHitARAProducer
 
detinfo::DetectorClocksData const fClockData
 
const double fTickPeriod
 
const double fBeamWindowStart
 
const double fBeamWindowEnd
 
const double fFlashStart
 
const double fFlashEnd
 
const unsigned fTimeBins
 
const bool fSelectNeutrino
 
const bool fOnlyCollectionWires
 
const bool fForceConcurrence
 
const bool fUseUncoatedPMT
 
const bool fUseOppVolMetric
 
const bool fUseARAPUCAS
 
const bool fStoreTrueNus
 
const bool fStoreCheatMCT0
 
const std::string fInputFilename
 
const bool fNoAvailableMetrics
 
const bool fMakeTree
 
const double fChargeToNPhotonsShower
 
const double fChargeToNPhotonsTrack
 
const double fMinHitQ
 
const double fMinSpacePointQ
 
const double fMinParticleQ
 
const double fMinSliceQ
 
const unsigned fMaxFlashes
 
const double fMinOpHPE
 
const double fMinFlashPE
 
const art::ServiceHandle
< geo::Geometry
fGeometry
 
const std::string fDetector
 
const bool fSBND
 
const bool fICARUS
 
const std::unique_ptr
< opdet::PDMapAlg
fPDMapAlgPtr
 
const size_t fNTPC
 
const int fCryostat
 
const std::unique_ptr
< geo::CryostatGeo
fGeoCryo
 
const std::list< double > fWiresX_gl
 
const double fDriftDistance
 
const int fXBins
 
const double fXBinWidth
 
const std::string fRR_TF1_fit
 
const std::string fRatio_TF1_fit
 
const unsigned fYBins
 
const unsigned fZBins
 
const double fYLow
 
const double fYHigh
 
const double fZLow
 
const double fZHigh
 
const double fYBiasSlope
 
const double fZBiasSlope
 
unsigned fDriftVolumes
 
unsigned fTPCPerDriftVolume
 
const unsigned fOpDetNormalizer
 
const double fTermThreshold
 
TTree * _flashmatch_nuslice_tree
 
TH2D * fRRH2
 
TH2D * fRatioH2
 
std::array< Fits, 3 > fRRFits
 
std::array< Fits, 3 > fRatioFits
 
const std::array< std::string, 3 > kSuffixes {"l", "h", "m"}
 
double _charge_x_gl
 
double _charge_x
 
double _charge_y
 
double _charge_z
 
double _charge_q
 
double _flash_x
 
double _flash_x_gl
 
double _flash_y
 
double _flash_yb
 
double _flash_z
 
double _flash_zb
 
double _flash_rr
 
double _flash_pe
 
double _flash_unpe
 
double _flash_ratio
 
double _flash_time
 
double _hypo_x
 
double _hypo_x_err
 
double _hypo_x_rr
 
double _hypo_x_ratio
 
double _y_skew
 
double _z_skew
 
double _score
 
double _scr_y
 
double _scr_z
 
double _scr_rr
 
double _scr_ratio
 
unsigned _evt
 
unsigned _run
 
unsigned _sub
 
unsigned _slices = -1
 
unsigned _true_nus = -1
 
double _mcT0 = -9999.
 
std::vector< double > fdYMeans
 
std::vector< double > fdZMeans
 
std::vector< double > fRRMeans
 
std::vector< double > fRatioMeans
 
std::vector< double > fdYSpreads
 
std::vector< double > fdZSpreads
 
std::vector< double > fRRSpreads
 
std::vector< double > fRatioSpreads
 
BookKeeping bk
 

Static Private Attributes

static constexpr unsigned kRght = 0
 
static constexpr unsigned kLeft = 1
 
static constexpr unsigned kActivityInRght = 100
 
static constexpr unsigned kActivityInLeft = 200
 
static constexpr unsigned kActivityInBoth = 300
 
static constexpr unsigned kMinEntriesInProjection = 100
 
static constexpr double kEps = 1e-4
 
static constexpr bool kNoScr = false
 
static constexpr double kNoScrTime = -9999.
 
static constexpr double kNoScrQ = -9999.
 
static constexpr double kNoScrPE = -9999.
 
static constexpr int kQNoOpHScr = -1
 
static constexpr int kNoChrgScr = -2
 
static constexpr int k0VUVPEScr = -3
 
static constexpr int kNotANuScr = -4
 
static constexpr int kNoOpHInEvt = -11
 
static constexpr int kNoPFPInEvt = -12
 
static constexpr int kNoSlcInEvt = -13
 

Detailed Description

Definition at line 64 of file FlashPredict.hh.

Member Typedef Documentation

Definition at line 199 of file FlashPredict.hh.

using FlashPredict::ChargeDigestMap = std::map<double, ChargeDigest, std::greater<double>>

Definition at line 106 of file FlashPredict.hh.

Definition at line 200 of file FlashPredict.hh.

using FlashPredict::OpHitIt = std::vector<recob::OpHit>::iterator

Definition at line 108 of file FlashPredict.hh.

Definition at line 201 of file FlashPredict.hh.

Definition at line 198 of file FlashPredict.hh.

Constructor & Destructor Documentation

FlashPredict::FlashPredict ( fhicl::ParameterSet const &  p)
explicit

Definition at line 12 of file FlashPredict_module.cc.

13  : EDProducer{p}
14  , fPandoraProducer(p.get<std::string>("PandoraProducer"))
15  , fSpacePointProducer(p.get<std::string>("SpacePointProducer"))
16  , fOpHitProducer(p.get<std::string>("OpHitProducer"))
17  , fOpHitARAProducer(p.get<std::string>("OpHitARAProducer", ""))
18  // , fCaloProducer(p.get<std::string>("CaloProducer"))
19  // , fTrackProducer(p.get<std::string>("TrackProducer"))
20  , fClockData(art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob())
22  , fBeamWindowStart(p.get<double>("BeamWindowStart")) //us
23  , fBeamWindowEnd(p.get<double>("BeamWindowEnd"))// us
24  , fFlashStart(p.get<double>("FlashStart")) // in us w.r.t. flash time
25  , fFlashEnd(p.get<double>("FlashEnd")) // in us w.r.t flash time
27  , fSelectNeutrino(p.get<bool>("SelectNeutrino", true)) // only attempt to match potential neutrino slices
28  , fOnlyCollectionWires(p.get<bool>("OnlyCollectionWires", true))
29  , fForceConcurrence(p.get<bool>("ForceConcurrence", false)) // require light and charge to coincide, different requirements for SBND and ICARUS
30  , fUseUncoatedPMT(p.get<bool>("UseUncoatedPMT", false))
31  , fUseOppVolMetric(p.get<bool>("UseOppVolMetric", false))
32  , fUseARAPUCAS(p.get<bool>("UseARAPUCAS", false))
33  , fStoreTrueNus(p.get<bool>("StoreTrueNus", false))
34  , fStoreCheatMCT0(p.get<bool>("StoreCheatMCT0", false))
35  // , fUseCalo(p.get<bool>("UseCalo", false))
36  , fInputFilename(p.get<std::string>("InputFileName")) // root file with score metrics
37  , fNoAvailableMetrics(p.get<bool>("NoAvailableMetrics", false))
38  , fMakeTree(p.get<bool>("MakeTree", false))
39  , fChargeToNPhotonsShower(p.get<double>("ChargeToNPhotonsShower", 1.0)) // ~40000/1600
40  , fChargeToNPhotonsTrack(p.get<double>("ChargeToNPhotonsTrack", 1.0)) // ~40000/1600
41  , fMinHitQ(p.get<double>("MinHitQ", 0.0))
42  , fMinSpacePointQ(p.get<double>("MinSpacePointQ", 0.0))
43  , fMinParticleQ(p.get<double>("MinParticleQ", 0.0))
44  , fMinSliceQ(p.get<double>("MinSliceQ", 0.0))
45  , fMaxFlashes(p.get<unsigned>("MaxFlashes", 1))
46  , fMinOpHPE(p.get<double>("MinOpHPE", 0.0))
47  , fMinFlashPE(p.get<double>("MinFlashPE", 0.0))
48  , fDetector(detectorName(fGeometry->DetectorName()))
49  , fSBND((fDetector == "SBND") ? true : false )
50  , fICARUS((fDetector == "ICARUS") ? true : false )
51  , fPDMapAlgPtr(art::make_tool<opdet::PDMapAlg>(p.get<fhicl::ParameterSet>("PDMapAlg")))
52  , fNTPC(fGeometry->NTPC())
53  , fCryostat(p.get<int>("Cryostat", 0)) //set =0 or =1 for ICARUS to match reco chain selection
54  , fGeoCryo(std::make_unique<geo::CryostatGeo>(fGeometry->Cryostat(fCryostat)))
55  , fWiresX_gl(wiresXGl())
57  , fXBins(p.get<double>("XBins"))
59  , fRR_TF1_fit(p.get<std::string>("rr_TF1_fit", "pol3"))
60  , fRatio_TF1_fit(p.get<std::string>("ratio_TF1_fit", "pol3"))
61  , fYBins(p.get<unsigned>("YBins", 0.)) // roughly match the rows of opdets
62  , fZBins(p.get<unsigned>("ZBins", 0.)) // roughly match the columns of opdets
63  , fYLow(p.get<double>("YLow", 0.)) // lowest opdet position in cm
64  , fYHigh(p.get<double>("YHigh", 0.)) // highest opdet position in cm
65  , fZLow(p.get<double>("ZLow", 0.)) // most upstream opdet position in cm
66  , fZHigh(p.get<double>("ZHigh", 0.)) // most downstream opdet position in cm
67  , fYBiasSlope(p.get<double>("YBiasSlope", 0.)) // correcting Y factor
68  , fZBiasSlope(p.get<double>("ZBiasSlope", 0.)) // correcting Z factor
69  , fOpDetNormalizer((fSBND) ? 4 : 1)
70  , fTermThreshold(p.get<double>("ThresholdTerm", 30.))
71 {
72  produces< std::vector<sbn::SimpleFlashMatch> >();
73  produces< art::Assns <recob::PFParticle, sbn::SimpleFlashMatch> >();
74 
75  // TODO: check params are sane:
76  if(fFlashStart > 0. || fFlashEnd < 0.){
77  throw cet::exception("FlashPredict")
78  << "fFlashStart has to be non-positive, "
79  << "and fFlashEnd has to be non-negative.";
80  }
81 
82  if(p.get<double>("DriftDistance") != driftDistance()){
83  mf::LogError("FlashPredict")
84  << "Provided driftDistance: " << p.get<double>("DriftDistance")
85  << " is different than expected: " << driftDistance();
86  }
87 
88  if(fSBND && !fICARUS) {
89  if(fUseOppVolMetric) {
90  throw cet::exception("FlashPredict")
91  << "UseOppVolMetric: " << std::boolalpha << fUseOppVolMetric << "\n"
92  << "Not supported on SBND. Stopping.";
93  }
96  }
97  else if(fICARUS && !fSBND) {
99  throw cet::exception("FlashPredict")
100  << "UseUncoatedPMT: " << std::boolalpha << fUseUncoatedPMT << ",\n"
101  << "UseARAPUCAS: " << std::boolalpha << fUseARAPUCAS << "\n"
102  << "Not supported on ICARUS. Stopping.\n";
103  }
104  fTPCPerDriftVolume = 2;
106  }
107  else {
108  throw cet::exception("FlashPredict")
109  << "Detector: " << fDetector
110  << ", not supported. Stopping.\n";
111  }
112 
113  if (fSBND && fCryostat != 0) {
114  throw cet::exception("FlashPredict")
115  << "SBND has only one cryostat. \n"
116  << "Check Detector and Cryostat parameter." << std::endl;
117  }
118  else if (fICARUS && fCryostat > 1) {
119  throw cet::exception("FlashPredict")
120  << "ICARUS has only two cryostats. \n"
121  << "Check Detector and Cryostat parameter." << std::endl;
122  }
123 
124  if (fMakeTree) initTree();
125 
126  loadMetrics();
127 
128  consumes<std::vector<recob::PFParticle>>(fPandoraProducer);
129  consumes<std::vector<recob::Slice>>(fPandoraProducer);
130  consumes<art::Assns<recob::SpacePoint, recob::PFParticle>>(fPandoraProducer);
131  consumes<std::vector<recob::SpacePoint>>(fSpacePointProducer);
132  consumes<art::Assns<recob::Hit, recob::SpacePoint>>(fSpacePointProducer);
133  consumes<std::vector<recob::OpHit>>(fOpHitProducer);
134  if(fUseARAPUCAS && !fOpHitARAProducer.empty())
135  consumes<std::vector<recob::OpHit>>(fOpHitARAProducer);
136 } // FlashPredict::FlashPredict(fhicl::ParameterSet const& p)
const bool fForceConcurrence
const double fZHigh
const bool fSBND
const int fXBins
detinfo::DetectorClocksData const fClockData
const double fMinSliceQ
const std::string fDetector
const double fBeamWindowStart
const double fTermThreshold
const art::ServiceHandle< geo::Geometry > fGeometry
pdgs p
Definition: selectors.fcl:22
const art::InputTag fOpHitARAProducer
const art::InputTag fPandoraProducer
const bool fICARUS
std::string detectorName(const std::string detName) const
const bool fMakeTree
const double fZLow
const double fXBinWidth
const double fZBiasSlope
constexpr double TickPeriod() const noexcept
A single tick period in microseconds.
Definition: ElecClock.h:352
const double fYHigh
const std::list< double > fWiresX_gl
const double fBeamWindowEnd
const bool fUseARAPUCAS
const double fMinParticleQ
const bool fUseUncoatedPMT
const std::string fInputFilename
const double fMinHitQ
const double fYBiasSlope
const std::string fRR_TF1_fit
const unsigned fTimeBins
const double fTickPeriod
void loadMetrics(void)
unsigned fTPCPerDriftVolume
double driftDistance() const
const double fMinOpHPE
const unsigned fYBins
ElecClock const & OpticalClock() const noexcept
Borrow a const Optical clock with time set to Trigger time [us].
std::list< double > wiresXGl() const
const double fFlashEnd
const bool fSelectNeutrino
const double fMinSpacePointQ
const double fDriftDistance
unsigned fDriftVolumes
const double fFlashStart
const art::InputTag fSpacePointProducer
const double fMinFlashPE
const double fYLow
const double fChargeToNPhotonsTrack
const art::InputTag fOpHitProducer
const unsigned fOpDetNormalizer
const bool fOnlyCollectionWires
const std::string fRatio_TF1_fit
const bool fStoreTrueNus
const unsigned fZBins
const bool fStoreCheatMCT0
const std::unique_ptr< opdet::PDMapAlg > fPDMapAlgPtr
const size_t fNTPC
const bool fNoAvailableMetrics
const std::unique_ptr< geo::CryostatGeo > fGeoCryo
const bool fUseOppVolMetric
const double fChargeToNPhotonsShower
const unsigned fMaxFlashes
const int fCryostat
FlashPredict::FlashPredict ( FlashPredict const &  )
delete
FlashPredict::FlashPredict ( FlashPredict &&  )
delete

Member Function Documentation

void FlashPredict::addDaughters ( const std::unordered_map< size_t, size_t > &  pfpMap,
const art::Ptr< recob::PFParticle > &  pfp_ptr,
const art::ValidHandle< std::vector< recob::PFParticle >> &  pfps_h,
std::vector< art::Ptr< recob::PFParticle >> &  pfp_v 
) const
private

Definition at line 1106 of file FlashPredict_module.cc.

1111 {
1112  mothers_daughters.push_back(pfp_ptr);
1113  const auto daughters = pfp_ptr->Daughters();
1114  for(const auto daughter : daughters) {
1115  const art::Ptr<recob::PFParticle> daughter_ptr(pfps_h, pfpMap.at(daughter));
1116  addDaughters(pfpMap, daughter_ptr, pfps_h, mothers_daughters);
1117  }
1118  return;
1119 }
void addDaughters(const std::unordered_map< size_t, size_t > &pfpMap, const art::Ptr< recob::PFParticle > &pfp_ptr, const art::ValidHandle< std::vector< recob::PFParticle >> &pfps_h, std::vector< art::Ptr< recob::PFParticle >> &pfp_v) const
void FlashPredict::beginJob ( )
override

Definition at line 1636 of file FlashPredict_module.cc.

1637 {
1638  bk = BookKeeping();
1639 }
BookKeeping bk
double FlashPredict::cheatMCT0 ( const std::vector< art::Ptr< recob::Hit >> &  hits,
const std::vector< art::Ptr< simb::MCParticle >> &  mcParticles 
)
private

Definition at line 627 of file FlashPredict_module.cc.

629 {
630  auto clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();
631  int pidMaxEnergy =
633  for(auto& mcp: mcParticles){
634  if(mcp->TrackId() == pidMaxEnergy){
635  return mcp->Position().T();
636  }
637  }
638  return -9999999.;
639 }
G4ID TrueParticleIDFromTotalTrueEnergy(detinfo::DetectorClocksData const &clockData, const std::vector< art::Ptr< recob::Hit >> &pHits, const bool rollupUnsavedIDs)
The G4 ID of the true particle which deposits the most energy in a vector of recob::Hit.
FlashPredict::ChargeMetrics FlashPredict::computeChargeMetrics ( const flashmatch::QCluster_t qClusters) const
private

Definition at line 642 of file FlashPredict_module.cc.

644 {
645  double xave = 0.; double yave = 0.;
646  double zave = 0.; double norm = 0.;
647  for (auto& qp : qClusters) {
648  xave += qp.q * qp.x;
649  yave += qp.q * qp.y;
650  zave += qp.q * qp.z;
651  norm += qp.q;
652  }
653  if (norm > 0.) {
654  double charge_q = norm;
655  double charge_x_gl = xave / norm;
656  double charge_x = foldXGl(charge_x_gl);
657  double charge_y = yave / norm;
658  double charge_z = zave / norm;
659  return {charge_x, charge_x_gl, charge_y, charge_z, charge_q, true};
660  }
661  return {};
662 }
double foldXGl(const double x_gl) const
auto norm(Vector const &v)
Return norm of the specified vector.
FlashPredict::FlashMetrics FlashPredict::computeFlashMetrics ( const SimpleFlash simpleFlash) const
private

Definition at line 665 of file FlashPredict_module.cc.

667 {
668  const OpHitIt opH_beg = simpleFlash.opH_beg;
669  const OpHitIt opH_end = simpleFlash.opH_end;
670 
671  std::unique_ptr<TH1F> ophY = std::make_unique<TH1F>("ophY", "", fYBins, fYLow, fYHigh);
672  std::unique_ptr<TH1F> ophZ = std::make_unique<TH1F>("ophZ", "", fZBins, fZLow, fZHigh);
673 
674  double peSumMax_wallX = wallXWithMaxPE(opH_beg, opH_end);
675 
676  double sum = 0.;
677  double sum_PE = 0.;
678  double sum_PE2 = 0.;
679  double sum_unPE = 0.;
680  double sum_visARA_PE = 0.;
681  double sum_PE2Y = 0.; double sum_PE2Z = 0.;
682  double sum_PE2Y2 = 0.; double sum_PE2Z2 = 0.;
683 
684  for(auto oph=opH_beg; oph!=opH_end; ++oph){
685  int opChannel = oph->OpChannel();
686  auto& opDet = fGeometry->OpDetGeoFromOpChannel(opChannel);
687  auto opDetXYZ = opDet.GetCenter();
688 
689  bool is_pmt_vis = false, is_ara_vis = false;
690  if(fSBND){// because VIS light
691  auto op_type = fPDMapAlgPtr->pdType(opChannel);
692  if(op_type == "pmt_uncoated") {
693  if(!fUseUncoatedPMT) continue;
694  is_pmt_vis = true, is_ara_vis = false;
695  }
696  else if(op_type == "xarapuca_vis" || op_type == "arapuca_vis") {
697  is_pmt_vis = false, is_ara_vis = true;
698  // if !fUseArapucas, they weren't loaded at all
699  }
700  }
701 
702  double ophPE = oph->PE();
703  double ophPE2 = ophPE * ophPE;
704  sum += 1.0;
705  sum_PE += ophPE;
706  sum_PE2 += ophPE2;
707  sum_PE2Y += ophPE2 * opDetXYZ.Y();
708  sum_PE2Z += ophPE2 * opDetXYZ.Z();
709  sum_PE2Y2 += ophPE2 * opDetXYZ.Y() * opDetXYZ.Y();
710  sum_PE2Z2 += ophPE2 * opDetXYZ.Z() * opDetXYZ.Z();
711 
712  ophY->Fill(opDetXYZ.Y(), ophPE);
713  ophZ->Fill(opDetXYZ.Z(), ophPE);
714 
715  if(fICARUS){
716  if(fUseOppVolMetric &&
717  std::abs(peSumMax_wallX-opDetXYZ.X()) > 5.) sum_unPE += ophPE;
718  }
719  else {// fSBND
720  if(fUseUncoatedPMT && is_pmt_vis) {
721  sum_unPE += ophPE;
722  }
723  else if(fUseARAPUCAS && is_ara_vis) {
724  sum_visARA_PE += ophPE;
725  }
726  }
727  } // for opHits
728 
729  if (sum_PE > 0.) {
730  FlashMetrics flash;
731  flash.metric_ok = true;
732  flash.pe = sum_PE;
733  flash.unpe = sum_unPE;
734  flash.ratio = fOpDetNormalizer * flash.unpe / flash.pe;
735  if(fUseARAPUCAS) {
736  flash.unpe += sum_visARA_PE;
737  flash.ratio = (fOpDetNormalizer * sum_unPE + sum_visARA_PE ) / flash.pe;
738  }
739  flash.yb = sum_PE2Y / sum_PE2;
740  flash.zb = sum_PE2Z / sum_PE2;
741  flash.rr = std::sqrt(
742  std::abs(sum_PE2Y2 + sum_PE2Z2 + sum_PE2 * (flash.yb * flash.yb + flash.zb * flash.zb)
743  - 2.0 * (flash.yb * sum_PE2Y + flash.zb * sum_PE2Z) ) / sum_PE2);
744  std::tie(flash.h_x, flash.h_xerr, flash.h_xrr, flash.h_xratio) =
745  hypoFlashX_H2(flash.rr, flash.ratio);
746 
747  // TODO: using _hypo_x make further corrections to _flash_time to
748  // account for light transport time and/or rising edge
749  // double flash_time = timeCorrections(simpleFlash.maxpeak_time, hypo_x);
750  flash.time = simpleFlash.maxpeak_time;
751  flash.x = peSumMax_wallX;
752  flash.x_gl = flashXGl(flash.h_x, flash.x);
753 
754  flash.y_skew = ophY->GetSkewness();
755  flash.z_skew = ophZ->GetSkewness();
756  if(!std::isnan(flash.h_x)){
757  double y_correction = 0.;
758  double z_correction = 0.;
759  const double skew_threshold = 10.; // TODO: figure out best value
760  if(fSBND) {
761  y_correction = (std::abs(flash.y_skew)<skew_threshold) ?
762  flash.y_skew * flash.h_x*flash.h_x * fYBiasSlope : 0.;
763  z_correction = (std::abs(flash.z_skew)<skew_threshold) ?
764  flash.z_skew * flash.h_x*flash.h_x * fZBiasSlope : 0.;
765  }
766  else{// fICARUS
767  y_correction = (std::abs(flash.y_skew)<skew_threshold) ?
768  flash.y_skew * flash.h_x * fYBiasSlope : 0.;
769  z_correction = (std::abs(flash.z_skew)<skew_threshold) ?
770  flash.z_skew * flash.h_x * fZBiasSlope : 0.;
771  }
772  flash.y = flash.yb - y_correction;
773  flash.z = flash.zb - z_correction;
774  }
775 
776  return flash;
777  }
778  else {
779  std::string channels;
780  for(auto oph=opH_beg; oph!=opH_end; ++oph) channels += std::to_string(oph->OpChannel()) + ' ';
781  // std::string tpcs;
782  // for(auto itpc: tpcWithHits) tpcs += std::to_string(itpc) + ' ';
783  mf::LogError("FlashPredict")
784  << "Really odd that I landed here, this shouldn't had happen.\n"
785  << "sum: \t" << sum << "\n"
786  << "sum_PE: \t" << sum_PE << "\n"
787  << "sum_unPE: \t" << sum_unPE << "\n"
788  // << "tpcWithHits: \t" << tpcs << "\n"
789  << "opHits size: \t" << std::distance(opH_beg, opH_end) << "\n"
790  << "channels: \t" << channels << std::endl;
791  return {};
792  }
793 }
const double fZHigh
const bool fSBND
const art::ServiceHandle< geo::Geometry > fGeometry
const bool fICARUS
const double fZLow
const double fZBiasSlope
const double fYHigh
T abs(T value)
const bool fUseARAPUCAS
double distance(geo::Point_t const &point, CathodeDesc_t const &cathode)
Returns the distance of a point from the cathode.
const bool fUseUncoatedPMT
const double fYBiasSlope
std::tuple< double, double, double, double > hypoFlashX_H2(double flash_rr, double flash_ratio) const
const unsigned fYBins
double wallXWithMaxPE(const OpHitIt opH_beg, const OpHitIt opH_end) const
const double fYLow
std::string to_string(WindowPattern const &pattern)
const unsigned fOpDetNormalizer
double flashXGl(const double hypo_x, const double flash_x) const
const unsigned fZBins
std::vector< recob::OpHit >::iterator OpHitIt
const std::unique_ptr< opdet::PDMapAlg > fPDMapAlgPtr
const bool fUseOppVolMetric
FlashPredict::Score FlashPredict::computeScore ( const ChargeMetrics charge,
const FlashMetrics flash,
const std::set< unsigned > &  tpcWithHits,
const int  pdgc 
) const
private

Definition at line 796 of file FlashPredict_module.cc.

801 {
802  double score = 0.;
803  unsigned tcount = 0;
804  int xbin = static_cast<int>(fXBins * (charge.x / fDriftDistance));
805  auto out = mf::LogWarning("FlashPredict");
806 
807  double scr_y = scoreTerm(flash.y, charge.y, fdYMeans[xbin], fdYSpreads[xbin]);
808  if(scr_y > fTermThreshold) printMetrics("Y", charge, flash, pdgc, tpcWithHits, scr_y, out);
809  score += scr_y;
810  tcount++;
811  double scr_z = scoreTerm(flash.z, charge.z, fdZMeans[xbin], fdZSpreads[xbin]);
812  if(scr_z > fTermThreshold) printMetrics("Z", charge, flash, pdgc, tpcWithHits, scr_z, out);
813  score += scr_z;
814  tcount++;
815  double scr_rr = scoreTerm(flash.rr, fRRMeans[xbin], fRRSpreads[xbin]);
816  if(scr_rr > fTermThreshold) printMetrics("RR", charge, flash, pdgc, tpcWithHits, scr_rr, out);
817  score += scr_rr;
818  tcount++;
819  double scr_ratio = 0.;
821  scr_ratio = scoreTerm(flash.ratio, fRatioMeans[xbin], fRatioSpreads[xbin]);
822  if(fICARUS && !std::isnan(flash.h_x)){
823  // HACK to penalise matches with flash and charge on opposite volumes
824  double x_gl_diff = std::abs(flash.x_gl-charge.x_gl);
825  double x_diff = std::abs(flash.h_x-charge.x);
826  double cathode_tolerance = 40.;
827  if(x_gl_diff > x_diff + cathode_tolerance) // ok if close to the cathode
828  scr_ratio += scoreTerm((flash.pe-flash.unpe)/flash.pe,
829  fRatioMeans[xbin], fRatioSpreads[xbin]);
830  }
831  if(scr_ratio > fTermThreshold) printMetrics("RATIO", charge, flash, pdgc, tpcWithHits, scr_ratio, out);
832  score += scr_ratio;
833  tcount++;
834  }
835  mf::LogDebug("FlashPredict")
836  << "score:\t" << score << "using " << tcount << " terms";
837  return {score, scr_y, scr_z, scr_rr, scr_ratio};
838 }
std::vector< double > fdZMeans
std::vector< double > fRRSpreads
const int fXBins
void printMetrics(const std::string metric, const ChargeMetrics &charge, const FlashMetrics &flash, const int pdgc, const std::set< unsigned > &tpcWithHits, const double term, Stream &&out) const
const double fTermThreshold
BEGIN_PROLOG or score(default)}sbnd_crttrackmatchingalg_crID
const bool fICARUS
std::vector< double > fRatioMeans
T abs(T value)
std::vector< double > fdYMeans
const bool fUseUncoatedPMT
double scoreTerm(const double m, const double n, const double mean, const double spread) const
std::vector< double > fdZSpreads
const double fDriftDistance
std::vector< double > fRRMeans
std::vector< double > fdYSpreads
std::vector< double > fRatioSpreads
const bool fUseOppVolMetric
void FlashPredict::copyOpHitsInBeamWindow ( std::vector< recob::OpHit > &  opHits,
const art::Handle< std::vector< recob::OpHit >> &  ophit_h 
) const
private

Definition at line 1199 of file FlashPredict_module.cc.

1202 {
1203  double s = fBeamWindowStart;
1204  double e = fBeamWindowEnd;
1205  double m = fMinOpHPE;
1206  // copy ophits that are inside the time window and with PEs
1207  auto opHitInWindow =
1208  [s, e, m, this](const recob::OpHit& oph)-> bool
1209  {return ((oph.PeakTime() > s) &&
1210  (oph.PeakTime() < e) &&
1211  (oph.PE() > m) &&
1212  isPDInCryo(oph.OpChannel())); };
1213  auto it = std::copy_if(ophits_h->begin(), ophits_h->end(), opHits.begin(),
1214  opHitInWindow);
1215  opHits.resize(std::distance(opHits.begin(), it));
1216 }
bool isPDInCryo(const int pdChannel) const
const double fBeamWindowStart
tuple m
now if test mode generate materials, CRT shell, world, gdml header else just generate CRT shell for u...
const double fBeamWindowEnd
double distance(geo::Point_t const &point, CathodeDesc_t const &cathode)
Returns the distance of a point from the cathode.
const double fMinOpHPE
then echo File list $list not found else cat $list while read file do echo $file sed s
Definition: file_to_url.sh:60
do i e
bool FlashPredict::createOpHitsTimeHist ( const std::vector< recob::OpHit > &  opHits,
std::vector< recob::OpHit > &  opHitsRght,
std::vector< recob::OpHit > &  opHitsLeft,
std::unique_ptr< TH1D > &  opHitsTimeHist,
std::unique_ptr< TH1D > &  opHitsTimeHistRght,
std::unique_ptr< TH1D > &  opHitsTimeHistLeft 
) const
private

Definition at line 1282 of file FlashPredict_module.cc.

1289 {
1290  for(const auto& oph : opHits) {
1291  auto ch = oph.OpChannel();
1292  if(!fUseUncoatedPMT &&
1293  fPDMapAlgPtr->isPDType(ch, "pmt_uncoated")) continue;
1294  opHitsTimeHist->Fill(oph.PeakTime(), oph.PE());
1295  if(sbndPDinTPC(ch) == kRght){
1296  opHitsRght.emplace_back(oph);
1297  opHitsTimeHistRght->Fill(oph.PeakTime(), oph.PE());
1298  }
1299  else{// sbndPDinTPC(ch) == kLeft
1300  opHitsLeft.emplace_back(oph);
1301  opHitsTimeHistLeft->Fill(oph.PeakTime(), oph.PE());
1302  }
1303  }
1304  if(opHitsTimeHist->GetEntries() <= 0 ||
1305  opHitsTimeHist->Integral() <= 0. ||
1306  opHitsTimeHist->Integral() <= fMinFlashPE) return false;
1307  if(opHitsTimeHistRght->GetEntries() <= 0 ||
1308  opHitsTimeHistRght->Integral() <= 0. ||
1309  opHitsTimeHistRght->Integral() <= fMinFlashPE)
1310  opHitsTimeHistRght->Reset();
1311  if(opHitsTimeHistLeft->GetEntries() <= 0 ||
1312  opHitsTimeHistLeft->Integral() <= 0. ||
1313  opHitsTimeHistLeft->Integral() <= fMinFlashPE)
1314  opHitsTimeHistLeft->Reset();
1315  return true;
1316 }
static constexpr unsigned kRght
const bool fUseUncoatedPMT
const double fMinFlashPE
const std::unique_ptr< opdet::PDMapAlg > fPDMapAlgPtr
unsigned sbndPDinTPC(const int pdChannel) const
unsigned FlashPredict::createOpHitsTimeHist ( const std::vector< recob::OpHit > &  opHits,
std::unique_ptr< TH1D > &  opHitsTimeHist 
) const
private

Definition at line 1320 of file FlashPredict_module.cc.

1323 {
1324  bool in_right = false, in_left = false;
1325  for(auto const& oph : opHits) {
1326  auto ch = oph.OpChannel();
1327  auto opDetXYZ = fGeometry->OpDetGeoFromOpChannel(ch).GetCenter();
1328  if(!fGeoCryo->ContainsPosition(opDetXYZ)) continue;
1329  opHitsTimeHist->Fill(oph.PeakTime(), oph.PE());
1330  unsigned t = icarusPDinTPC(ch);
1331  if(t/fTPCPerDriftVolume == kRght) in_right = true;
1332  else if(t/fTPCPerDriftVolume == kLeft) in_left = true;
1333  }
1334  if(opHitsTimeHist->GetEntries() <= 0 ||
1335  opHitsTimeHist->Integral() <= 0. ||
1336  opHitsTimeHist->Integral() <= fMinFlashPE) return 0;
1337  if(in_right && in_left) return kActivityInBoth;
1338  else if(in_right && !in_left) return kActivityInRght;
1339  else if(!in_right && in_left) return kActivityInLeft;
1340  else return 0;
1341 }
static constexpr unsigned kActivityInRght
static constexpr unsigned kRght
const art::ServiceHandle< geo::Geometry > fGeometry
static constexpr unsigned kLeft
static constexpr unsigned kActivityInBoth
unsigned icarusPDinTPC(const int pdChannel) const
unsigned fTPCPerDriftVolume
static constexpr unsigned kActivityInLeft
const double fMinFlashPE
const std::unique_ptr< geo::CryostatGeo > fGeoCryo
std::string FlashPredict::detectorName ( const std::string  detName) const
inlineprivate

Definition at line 1386 of file FlashPredict_module.cc.

1387 {
1388  if(detName.find("sbnd") != std::string::npos) return "SBND";
1389  if(detName.find("icarus") != std::string::npos) return "ICARUS";
1390  return "";
1391 }
double FlashPredict::driftDistance ( ) const
private

Definition at line 1476 of file FlashPredict_module.cc.

1477 {
1478  auto wit = fWiresX_gl.begin();
1479  auto wite = fWiresX_gl.begin();
1480  wite++;
1481  return std::abs(*wite - *wit)/2.;
1482 }
T abs(T value)
const std::list< double > fWiresX_gl
unsigned FlashPredict::driftVolume ( const double  x) const
private

Definition at line 1520 of file FlashPredict_module.cc.

1521 {
1522  auto wit = fWiresX_gl.begin();
1523  for(size_t i=0; i<fWiresX_gl.size()-1; i++){
1524  double wxl = *wit;
1525  wit++;
1526  double wxh = *wit;
1527  if(wxl < x && x<= wxh){
1528  double mid = (wxl + wxh)/2.;
1529  return (x<=mid) ? i*fCryostat : (i*fCryostat)+1;
1530  }
1531  }
1532  return 100;
1533 }
process_name opflash particleana ie x
const std::list< double > fWiresX_gl
const int fCryostat
void FlashPredict::endJob ( )
override

Definition at line 1641 of file FlashPredict_module.cc.

1642 {
1643  printBookKeeping(mf::LogWarning("FlashPredict"));
1644 }
void printBookKeeping(Stream &&out)
bool FlashPredict::findSimpleFlashes ( std::vector< SimpleFlash > &  simpleFlashes,
std::vector< recob::OpHit > &  opHits,
const unsigned  ophsInVolume,
std::unique_ptr< TH1D > &  opHitsTimeHist 
) const
private

Definition at line 1344 of file FlashPredict_module.cc.

1349 {
1350  OpHitIt opH_beg = opHits.begin();
1351  for(unsigned flashId=0; flashId<fMaxFlashes; ++flashId){
1352  int ibin = opHitsTimeHist->GetMaximumBin();
1353  double maxpeak_time = opHitsTimeHist->GetBinCenter(ibin);
1354  double lowedge = maxpeak_time + fFlashStart;
1355  double highedge = maxpeak_time + fFlashEnd;
1356  mf::LogDebug("FlashPredict")
1357  << "light window " << lowedge << " " << highedge << std::endl;
1358  int lowedge_bin = opHitsTimeHist->FindBin(lowedge);
1359  int highedge_bin = opHitsTimeHist->FindBin(highedge);
1360  // check if flash has enough PEs, return if is the first one
1361  if (opHitsTimeHist->Integral(lowedge_bin, highedge_bin) <= fMinFlashPE ||
1362  opHitsTimeHist->Integral(lowedge_bin, highedge_bin) <= 0. ){
1363  if(flashId == 0) return false;
1364  break;
1365  }
1366  // clear this peak to enforce non-overlapping flashes
1367  for(int i=lowedge_bin; i<highedge_bin; ++i){
1368  opHitsTimeHist->SetBinContent(i, 0.);
1369  }
1370  auto peakInsideEdges =
1371  [&lowedge, &highedge](const recob::OpHit& oph)-> bool
1372  { return ((lowedge <= oph.PeakTime()) && (oph.PeakTime() <= highedge)); };
1373  // partition container to move the hits of the flash
1374  // the iterators point to the boundaries of the partition
1375  OpHitIt opH_end = std::partition(opH_beg, opHits.end(),
1376  peakInsideEdges);
1377  simpleFlashes.emplace_back
1378  (SimpleFlash(flashId, ophsInVolume,
1379  opH_beg, opH_end, maxpeak_time));
1380  opH_beg = opH_end;
1381  }
1382  return true;
1383 }
const double fFlashEnd
const double fFlashStart
const double fMinFlashPE
std::vector< recob::OpHit >::iterator OpHitIt
const unsigned fMaxFlashes
double FlashPredict::flashXGl ( const double  hypo_x,
const double  flash_x 
) const
private

Definition at line 1485 of file FlashPredict_module.cc.

1487 {
1488  double min = 10000;
1489  unsigned wId = 0;
1490  auto wIt = fWiresX_gl.begin();
1491  auto w = wIt;
1492  for(size_t i=0; i<fWiresX_gl.size(); i++){
1493  if(((*w<0) == (flash_x<0)) && std::abs(*w - flash_x) < min) {
1494  wIt = w;
1495  wId = i;
1496  min = std::abs(*w - flash_x);
1497  }
1498  w++;;
1499  }
1500  return (wId % 2) ? (*wIt - hypo_x) : (*wIt + hypo_x);
1501 }
T abs(T value)
const std::list< double > fWiresX_gl
double FlashPredict::foldXGl ( const double  x_gl) const
private

Definition at line 1504 of file FlashPredict_module.cc.

1505 {
1506  auto wit = fWiresX_gl.begin();
1507  for(size_t i=0; i<fWiresX_gl.size()-1; i++){
1508  double wxl = *wit;
1509  wit++;
1510  double wxh = *wit;
1511  if(wxl < x_gl && x_gl<= wxh){
1512  double mid = (wxl + wxh)/2.;
1513  return (x_gl<=mid) ? std::abs(x_gl-wxl) : std::abs(x_gl-wxh);
1514  }
1515  }
1516  return -10.;
1517 }
double std(const std::vector< short > &wf, const double ped_mean, size_t start, size_t nsample)
Definition: UtilFunc.cxx:42
T abs(T value)
const std::list< double > fWiresX_gl
std::tuple< double, double, double, double > FlashPredict::hypoFlashX_fits ( double  flash_rr,
double  flash_ratio 
) const
private

Definition at line 841 of file FlashPredict_module.cc.

843 {
844  std::vector<double> rrXs;
845  double rr_hypoX, rr_hypoXWgt;
846  for(const auto& rrF : fRRFits){
847  if(rrF.min < flash_rr && flash_rr < rrF.max){
848  try{
849  rrXs.emplace_back(rrF.f->GetX(flash_rr, 0., fDriftDistance, kEps));
850  }catch (...) {
851  mf::LogWarning("FlashPredict")
852  << "Couldn't find root with fRRFits.\n"
853  << "min/_flash_rr/max:"
854  << rrF.min << "/" << flash_rr << "/" << rrF.max;
855  }
856  }
857  }
858  if(rrXs.size() > 1){//between: [l,h], [l,m], or [h,m]
859  rr_hypoX = (rrXs[0] + rrXs[1])/2.;
860  double half_interval = (rrXs[1] - rrXs[0])/2.;
861  rr_hypoXWgt = 1./(half_interval*half_interval);
862  }
863  else if(rrXs.size() == 0){//can't estimate
864  rr_hypoX = 0.;
865  rr_hypoXWgt = 0.;
866  }
867  else{//(rrXs.size() == 1)
868  if(flash_rr < fRRFits[2].min){//between: [l, m)
869  rr_hypoX = rrXs[0]/2.;
870  double half_interval = rrXs[0]/2.;
871  rr_hypoXWgt = 1./(half_interval*half_interval);
872  }
873  else{//between: (m, h]
874  rr_hypoX = (rrXs[0] + fDriftDistance)/2.;
875  double half_interval = (fDriftDistance - rrXs[0]);
876  rr_hypoXWgt = 1./(half_interval*half_interval);
877  }
878  }
879 
881  return {rr_hypoX, rr_hypoXWgt, rr_hypoX, 0.};
882 
883  std::vector<double> ratioXs;
884  double ratio_hypoX, ratio_hypoXWgt;
885  for(const auto& ratioF : fRatioFits){
886  if(ratioF.min < flash_ratio && flash_ratio < ratioF.max){
887  try{
888  ratioXs.emplace_back(ratioF.f->GetX(flash_ratio, 0., fDriftDistance, kEps));
889  }catch (...) {
890  mf::LogWarning("FlashPredict")
891  << "Couldn't find root with fRatioFits.\n"
892  << "min/flash_ratio/max:"
893  << ratioF.min << "/" << flash_ratio << "/" << ratioF.max;
894  }
895  }
896  }
897  if(ratioXs.size() > 1){//between: [l,h], [l,m], or [h,m]
898  ratio_hypoX = (ratioXs[0] + ratioXs[1])/2.;
899  double half_interval = (ratioXs[1] - ratioXs[0])/2.;
900  ratio_hypoXWgt = 1./(half_interval*half_interval);
901  }
902  else if(ratioXs.size() == 0){//can't estimate
903  ratio_hypoX = 0.;
904  ratio_hypoXWgt = 0.;
905  }
906  else{//(ratioXs.size() == 1)
907  if(flash_ratio < fRatioFits[2].min){//between: [l, m)
908  ratio_hypoX = ratioXs[0]/2.;
909  double half_interval = ratioXs[0]/2.;
910  ratio_hypoXWgt = 1./(half_interval*half_interval);
911  }
912  else{//between: (m, h]
913  ratio_hypoX = (ratioXs[0] + fDriftDistance)/2.;
914  double half_interval = (fDriftDistance - ratioXs[0]);
915  ratio_hypoXWgt = 1./(half_interval*half_interval);
916  }
917  }
918 
919  double sum_weights = rr_hypoXWgt + ratio_hypoXWgt;
920  double hypo_x =
921  (rr_hypoX*rr_hypoXWgt + ratio_hypoX*ratio_hypoXWgt) / sum_weights;
922  double hypo_x_err = std::sqrt(sum_weights) / sum_weights;
923  return {hypo_x, hypo_x_err, rr_hypoX, ratio_hypoX};
924 }
const bool fUseUncoatedPMT
std::array< Fits, 3 > fRRFits
static constexpr double kEps
const double fDriftDistance
std::array< Fits, 3 > fRatioFits
const bool fUseOppVolMetric
std::tuple< double, double, double, double > FlashPredict::hypoFlashX_H2 ( double  flash_rr,
double  flash_ratio 
) const
private

Definition at line 927 of file FlashPredict_module.cc.

929 {
930  auto[rr_hypoX, rr_hypoXWgt] = xEstimateAndRMS(flash_rr, fRRH2);
931  auto[ratio_hypoX, ratio_hypoXWgt] = xEstimateAndRMS(flash_ratio, fRatioH2);
932 
933  double sum_weights = rr_hypoXWgt + ratio_hypoXWgt;
934  double hypo_x =
935  (rr_hypoX*rr_hypoXWgt + ratio_hypoX*ratio_hypoXWgt) / sum_weights;
936  double hypo_x_err = std::sqrt(sum_weights) / sum_weights;
937  return {hypo_x, hypo_x_err, rr_hypoX, ratio_hypoX};
938 }
std::tuple< double, double > xEstimateAndRMS(double metric_value, const TH2D *metric_h2) const
unsigned FlashPredict::icarusPDinTPC ( const int  pdChannel) const
private

Definition at line 1409 of file FlashPredict_module.cc.

1410 {
1411  auto p = fGeometry->OpDetGeoFromOpChannel(pdChannel).GetCenter();
1412  if(fCryostat == 0) p.SetX((p.X() + 222.)/2. - 222.);//OpDets are outside the TPCs
1413  if(fCryostat == 1) p.SetX((p.X() - 222.)/2. + 222.);//OpDets are outside the TPCs
1414  return (fGeometry->PositionToTPCID(p)).TPC;
1415 }
const art::ServiceHandle< geo::Geometry > fGeometry
pdgs p
Definition: selectors.fcl:22
BEGIN_PROLOG TPC
const int fCryostat
void FlashPredict::initTree ( void  )
private

Definition at line 422 of file FlashPredict_module.cc.

423 {
424  art::ServiceHandle<art::TFileService> tfs;
425  _flashmatch_nuslice_tree = tfs->make<TTree>("nuslicetree", "nu FlashPredict tree");
426  _flashmatch_nuslice_tree->Branch("evt", &_evt, "evt/I");
427  _flashmatch_nuslice_tree->Branch("run", &_run, "run/I");
428  _flashmatch_nuslice_tree->Branch("sub", &_sub, "sub/I");
429  _flashmatch_nuslice_tree->Branch("slices", &_slices, "slices/I");
430  _flashmatch_nuslice_tree->Branch("flash_time", &_flash_time, "flash_time/D");
431  _flashmatch_nuslice_tree->Branch("flash_x_gl", &_flash_x_gl, "flash_x_gl/D");
432  _flashmatch_nuslice_tree->Branch("flash_x", &_flash_x, "flash_x/D");
433  _flashmatch_nuslice_tree->Branch("flash_y", &_flash_y, "flash_y/D");
434  _flashmatch_nuslice_tree->Branch("flash_yb", &_flash_yb, "flash_yb/D");
435  _flashmatch_nuslice_tree->Branch("flash_z", &_flash_z, "flash_z/D");
436  _flashmatch_nuslice_tree->Branch("flash_zb", &_flash_zb, "flash_zb/D");
437  _flashmatch_nuslice_tree->Branch("flash_rr", &_flash_rr, "flash_rr/D");
438  _flashmatch_nuslice_tree->Branch("flash_pe", &_flash_pe, "flash_pe/D");
439  _flashmatch_nuslice_tree->Branch("flash_unpe", &_flash_unpe, "flash_unpe/D");
440  _flashmatch_nuslice_tree->Branch("flash_ratio", &_flash_ratio, "flash_ratio/D");
441  _flashmatch_nuslice_tree->Branch("charge_x_gl", &_charge_x_gl, "charge_x_gl/D");
442  _flashmatch_nuslice_tree->Branch("charge_x", &_charge_x, "charge_x/D");
443  _flashmatch_nuslice_tree->Branch("charge_y", &_charge_y, "charge_y/D");
444  _flashmatch_nuslice_tree->Branch("charge_z", &_charge_z, "charge_z/D");
445  _flashmatch_nuslice_tree->Branch("charge_q", &_charge_q, "charge_q/D");
446  _flashmatch_nuslice_tree->Branch("hypo_x", &_hypo_x, "hypo_x/D");
447  _flashmatch_nuslice_tree->Branch("hypo_x_err", &_hypo_x_err, "hypo_x_err/D");
448  _flashmatch_nuslice_tree->Branch("hypo_x_rr", &_hypo_x_rr, "hypo_x_rr/D");
449  _flashmatch_nuslice_tree->Branch("hypo_x_ratio", &_hypo_x_ratio, "hypo_x_ratio/D");
450  _flashmatch_nuslice_tree->Branch("score", &_score, "score/D");
451  _flashmatch_nuslice_tree->Branch("scr_y", &_scr_y, "scr_y/D");
452  _flashmatch_nuslice_tree->Branch("scr_z", &_scr_z, "scr_z/D");
453  _flashmatch_nuslice_tree->Branch("scr_rr", &_scr_rr, "scr_rr/D");
454  _flashmatch_nuslice_tree->Branch("scr_ratio", &_scr_ratio, "scr_ratio/D");
455  _flashmatch_nuslice_tree->Branch("y_skew", &_y_skew, "y_skew/D");
456  _flashmatch_nuslice_tree->Branch("z_skew", &_z_skew, "z_skew/D");
457  _flashmatch_nuslice_tree->Branch("true_nus", &_true_nus, "true_nus/I");
458  _flashmatch_nuslice_tree->Branch("mcT0", &_mcT0, "mcT0/D");
459 }
double _charge_q
double _charge_x_gl
double _flash_time
unsigned _sub
double _flash_zb
double _charge_y
double _charge_x
double _scr_ratio
double _hypo_x_err
double _flash_pe
double _flash_unpe
double _flash_x_gl
TTree * _flashmatch_nuslice_tree
unsigned _evt
double _flash_ratio
double _hypo_x_rr
unsigned _slices
art::ServiceHandle< art::TFileService > tfs
double _flash_rr
double _charge_z
double _hypo_x_ratio
unsigned _true_nus
unsigned _run
double _flash_yb
bool FlashPredict::isPDInCryo ( const int  pdChannel) const
private

Definition at line 1394 of file FlashPredict_module.cc.

1395 {
1396  if(fSBND) return true;
1397  else { // fICARUS
1398  // BUG: I believe this function is not working, every now and then
1399  // I get ophits from the other cryo
1400  auto& p = fGeometry->OpDetGeoFromOpChannel(pdChannel).GetCenter();
1401  // if the channel is in the Cryostat is relevant
1402  return fGeoCryo->ContainsPosition(p);
1403  }
1404  return false;
1405 }
const bool fSBND
const art::ServiceHandle< geo::Geometry > fGeometry
pdgs p
Definition: selectors.fcl:22
const std::unique_ptr< geo::CryostatGeo > fGeoCryo
bool FlashPredict::isSBNDPDRelevant ( const int  pdChannel,
const std::set< unsigned > &  tpcWithHits 
) const
private

Definition at line 1418 of file FlashPredict_module.cc.

1420 {
1421  // if there's hits on all TPCs all channels are relevant
1422  if(tpcWithHits.size() == fNTPC) return true;
1423  unsigned pdTPC = sbndPDinTPC(pdChannel);
1424  for(auto itpc: tpcWithHits) if(itpc == pdTPC) return true;
1425  return false;
1426 }
const size_t fNTPC
unsigned sbndPDinTPC(const int pdChannel) const
void FlashPredict::loadMetrics ( void  )
private

Definition at line 462 of file FlashPredict_module.cc.

463 {
464  // TODO: fill histos with less repetition and range for loops
465  // read histograms and fill vectors for match score calculation
466  std::string fname;
467  cet::search_path sp("FW_SEARCH_PATH");
468  if(!sp.find_file(fInputFilename, fname)) {
469  mf::LogError("FlashPredict")
470  << "Could not find the light-charge match root file '"
471  << fInputFilename << "' on FW_SEARCH_PATH\n";
472  throw cet::exception("FlashPredict")
473  << "Could not find the light-charge match root file '"
474  << fInputFilename << "' on FW_SEARCH_PATH\n";
475  }
476  mf::LogInfo("FlashPredict") << "Opening file with metrics: " << fname;
477  TFile *infile = new TFile(fname.c_str(), "READ");
478  auto metricsInFile = infile->GetListOfKeys();
479  if(!metricsInFile->Contains("dy_h1") ||
480  !metricsInFile->Contains("dz_h1") ||
481  !metricsInFile->Contains("rr_h1") ||
482  !metricsInFile->Contains("ratio_h1") ||
483  !metricsInFile->Contains("rr_fit_l") ||
484  !metricsInFile->Contains("rr_fit_m") ||
485  !metricsInFile->Contains("rr_fit_h") ||
486  !metricsInFile->Contains("ratio_fit_l") ||
487  !metricsInFile->Contains("ratio_fit_m") ||
488  !metricsInFile->Contains("ratio_fit_h"))
489  {
490  mf::LogError("FlashPredict")
491  << "The metrics file '" << fname << "' lacks at least one metric.";
492  throw cet::exception("FlashPredict")
493  << "The metrics file '" << fname << "'lacks at least one metric.";
494  }
495  //
496  TH1 *temphisto = (TH1*)infile->Get("dy_h1");
497  int bins = temphisto->GetNbinsX();
498  if (bins <= 0 || fNoAvailableMetrics) {
499  std::cout << " problem with input histos for dy " << bins << " bins " << std::endl;
500  bins = 1;
501  fdYMeans.push_back(0);
502  fdYSpreads.push_back(0.001);
503  }
504  else {
505  for (int ib = 1; ib <= bins; ++ib) {
506  fdYMeans.push_back(temphisto->GetBinContent(ib));
507  double tt = temphisto->GetBinError(ib);
508  if (tt <= 0) {
509  std::cout << "zero value for bin spread in dy" << std::endl;
510  std::cout << "ib:\t" << ib << "\n";
511  std::cout << "temphisto->GetBinContent(ib):\t" << temphisto->GetBinContent(ib) << "\n";
512  std::cout << "temphisto->GetBinError(ib):\t" << temphisto->GetBinError(ib) << "\n";
513  tt = 100.;
514  }
515  fdYSpreads.push_back(tt);
516  }
517  }
518  //
519  temphisto = (TH1*)infile->Get("dz_h1");
520  bins = temphisto->GetNbinsX();
521  if (bins <= 0 || fNoAvailableMetrics) {
522  std::cout << " problem with input histos for dz " << bins << " bins " << std::endl;
523  bins = 1;
524  fdZMeans.push_back(0);
525  fdZSpreads.push_back(0.001);
526  }
527  else {
528  for (int ib = 1; ib <= bins; ++ib) {
529  fdZMeans.push_back(temphisto->GetBinContent(ib));
530  double tt = temphisto->GetBinError(ib);
531  if (tt <= 0) {
532  std::cout << "zero value for bin spread in dz" << std::endl;
533  std::cout << "ib:\t" << ib << "\n";
534  std::cout << "temphisto->GetBinContent(ib):\t" << temphisto->GetBinContent(ib) << "\n";
535  std::cout << "temphisto->GetBinError(ib):\t" << temphisto->GetBinError(ib) << "\n";
536  tt = 100.;
537  }
538  fdZSpreads.push_back(tt);
539  }
540  }
541  //
542  temphisto = (TH1*)infile->Get("rr_h1");
543  bins = temphisto->GetNbinsX();
544  if (bins <= 0 || fNoAvailableMetrics) {
545  std::cout << " problem with input histos for rr " << bins << " bins " << std::endl;
546  bins = 1;
547  fRRMeans.push_back(0);
548  fRRSpreads.push_back(0.001);
549  }
550  else {
551  for (int ib = 1; ib <= bins; ++ib) {
552  double me = temphisto->GetBinContent(ib);
553  fRRMeans.push_back(me);
554  double tt = temphisto->GetBinError(ib);
555  if (tt <= 0) {
556  std::cout << "zero value for bin spread in rr" << std::endl;
557  std::cout << "ib:\t" << ib << "\n";
558  std::cout << "temphisto->GetBinContent(ib):\t" << temphisto->GetBinContent(ib) << "\n";
559  std::cout << "temphisto->GetBinError(ib):\t" << temphisto->GetBinError(ib) << "\n";
560  tt = 100.;
561  }
562  fRRSpreads.push_back(tt);
563  }
564  unsigned s = 0;
565  for(auto& rrF : fRRFits){
566  std::string nold = "rr_fit_" + kSuffixes[s];
567  std::string nnew = "rrFit_" + kSuffixes[s];
568  TF1* tempF1 = (TF1*)infile->Get(nold.c_str());
569  auto params = tempF1->GetParameters();
570  rrF.f = std::make_unique<TF1>(nnew.c_str(), fRR_TF1_fit.c_str(),
571  0., fDriftDistance);
572  rrF.f->SetParameters(params);
573  rrF.min = rrF.f->GetMinimum(0., fDriftDistance, kEps);
574  rrF.max = rrF.f->GetMaximum(0., fDriftDistance, kEps);
575  s++;
576  }
577  }
578  //
579  temphisto = (TH1*)infile->Get("ratio_h1");
580  bins = temphisto->GetNbinsX();
581  if (bins <= 0 || fNoAvailableMetrics) {
582  std::cout << " problem with input histos for ratio " << bins << " bins " << std::endl;
583  bins = 1;
584  fRatioMeans.push_back(0);
585  fRatioSpreads.push_back(0.001);
586  }
587  else {
588  for (int ib = 1; ib <= bins; ++ib) {
589  double me = temphisto->GetBinContent(ib);
590  fRatioMeans.push_back(me);
591  double tt = temphisto->GetBinError(ib);
592  if (tt <= 0) {
593  std::cout << "zero value for bin spread in ratio" << std::endl;
594  std::cout << "ib:\t" << ib << "\n";
595  std::cout << "temphisto->GetBinContent(ib):\t" << temphisto->GetBinContent(ib) << "\n";
596  std::cout << "temphisto->GetBinError(ib):\t" << temphisto->GetBinError(ib) << "\n";
597  tt = 100.;
598  }
599  fRatioSpreads.push_back(tt);
600  }
601  unsigned s = 0;
602  for(auto& ratioF : fRatioFits){
603  std::string nold = "ratio_fit_" + kSuffixes[s];
604  std::string nnew = "ratioFit_" + kSuffixes[s];
605  TF1* tempF1 = (TF1*)infile->Get(nold.c_str());
606  auto params = tempF1->GetParameters();
607  ratioF.f = std::make_unique<TF1>(nnew.c_str(), fRatio_TF1_fit.c_str(),
608  0., fDriftDistance);
609  ratioF.f->SetParameters(params);
610  ratioF.min = ratioF.f->GetMinimum(0., fDriftDistance, kEps);
611  ratioF.max = ratioF.f->GetMaximum(0., fDriftDistance, kEps);
612  s++;
613  }
614  }
615 
616  TH2* tmp1_h2 = (TH2*)infile->Get("rr_h2");
617  fRRH2 = (TH2D*)tmp1_h2->Clone("fRRH2");
618  TH2* tmp2_h2 = (TH2*)infile->Get("ratio_h2");
619  fRatioH2 = (TH2D*)tmp2_h2->Clone("fRatioH2");
620 
621  infile->Close();
622  delete infile;
623  mf::LogInfo("FlashPredict") << "Finish loading metrics";
624 }
std::vector< double > fdZMeans
std::vector< double > fRRSpreads
string fname
Definition: demo.py:5
const std::array< std::string, 3 > kSuffixes
std::vector< double > fRatioMeans
std::vector< double > fdYMeans
const std::string fInputFilename
std::array< Fits, 3 > fRRFits
const std::string fRR_TF1_fit
static constexpr double kEps
std::vector< double > fdZSpreads
const double fDriftDistance
std::array< Fits, 3 > fRatioFits
std::vector< double > fRRMeans
then echo File list $list not found else cat $list while read file do echo $file sed s
Definition: file_to_url.sh:60
std::vector< double > fdYSpreads
const std::string fRatio_TF1_fit
std::vector< double > fRatioSpreads
const bool fNoAvailableMetrics
BEGIN_PROLOG could also be cout
FlashPredict::ChargeDigestMap FlashPredict::makeChargeDigest ( const art::Event &  evt,
const art::ValidHandle< std::vector< recob::PFParticle >> &  pfps_h 
)
private

Definition at line 976 of file FlashPredict_module.cc.

979 {
980  // grab spacepoints associated with PFParticles
981  const art::FindManyP<recob::SpacePoint>
982  pfp_spacepoints_assns(pfps_h, evt, fPandoraProducer);
983  const auto& spacepoints_h =
984  evt.getValidHandle<std::vector<recob::SpacePoint>>(fSpacePointProducer);
985  const art::FindManyP<recob::Hit>
986  spacepoint_hits_assns(spacepoints_h, evt, fSpacePointProducer);
987  // grab tracks associated with PFParticles
988  // auto const& track_h = evt.getValidHandle<std::vector<recob::Track> >(fTrackProducer);
989  // art::FindManyP<recob::Track> pfp_track_assn_v(track_h, evt, fTrackProducer);
990  // grab calorimetry info for tracks
991  // auto const& calo_h =
992  // evt.getValidHandle<std::vector<anab::Calorimetry> >(fCaloProducer);
993  // art::FindManyP<anab::Calorimetry> track_calo_assn_v(calo_h,
994  // evt, fCaloProducer);
995 
996  std::vector<art::Ptr<simb::MCParticle>> mcParticles;
997  if(fStoreCheatMCT0){
999  mcParticles);
1000  }
1001 
1002  std::unordered_map<size_t, size_t> pfpMap;
1003  for(size_t pId=0; pId<pfps_h->size(); pId++) {
1004  pfpMap[pfps_h->at(pId).Self()] = pId;
1005  }
1006 
1007  // TODO: this block is meant to only load the charge related objects
1008  // of the slices that are potentially neutrinos, this to improve
1009  // performance. Loading these objects is currently what that takes
1010  // the most time. Unfortunately, the objects are not very robust and
1011  // frequently there's out of bounds requests.
1012  //
1013  // std::vector<art::Ptr<recob::PFParticle>> particles_in_slices;
1014  // for(size_t pId=0; pId<pfps_h->size(); pId++) {
1015  // if(!pfps_h->at(pId).IsPrimary()) continue;
1016  // const art::Ptr<recob::PFParticle> pfp_ptr(pfps_h, pId);
1017  // unsigned pfpPDGC = std::abs(pfp_ptr->PdgCode());
1018  // if(fSelectNeutrino &&
1019  // (pfpPDGC != 12) && (pfpPDGC != 14) && (pfpPDGC != 16)) continue;
1020  // std::vector<art::Ptr<recob::PFParticle>> particles;
1021  // addDaughters(pfpMap, pfp_ptr, pfps_h, particles);
1022  // particles_in_slices.insert(particles_in_slices.end(),
1023  // particles.begin(), particles.end());
1024  // }
1025  // const art::FindManyP<recob::SpacePoint>
1026  // pfp_spacepoints_assns(particles_in_slices, evt, fPandoraProducer);
1027  // std::vector<art::Ptr<recob::SpacePoint>> spacepoints_in_slices;
1028  // // for(auto& p: particles_in_slices) {
1029  // for(size_t p=0; p<particles_in_slices.size(); ++p) {
1030  // auto& sps = pfp_spacepoints_assns.at(p);
1031  // spacepoints_in_slices.insert(spacepoints_in_slices.end(), sps.begin(), sps.end());
1032  // }
1033  // const art::FindManyP<recob::Hit>
1034  // spacepoint_hits_assns(spacepoints_in_slices, evt, fSpacePointProducer);
1035  // TODO: this block is meant to only load the charge related objects
1036 
1037  ChargeDigestMap chargeDigestMap;
1038  // Loop over pandora pfp particles
1039  for(size_t pId=0; pId<pfps_h->size(); pId++) {
1040  if(!pfps_h->at(pId).IsPrimary()) continue;
1041  const art::Ptr<recob::PFParticle> pfp_ptr(pfps_h, pId);
1042  unsigned pfpPDGC = std::abs(pfp_ptr->PdgCode());
1043  if(fSelectNeutrino &&
1044  (pfpPDGC != 12) && (pfpPDGC != 14) && (pfpPDGC != 16)){
1045  chargeDigestMap[-10.-pId] =
1046  ChargeDigest(pId, pfpPDGC, pfp_ptr, flashmatch::QCluster_t{}, std::set<unsigned>{}, -9999.);
1047  continue;
1048  }
1049  std::set<unsigned> tpcWithHits;
1050 
1051  std::vector<art::Ptr<recob::PFParticle>> particles_in_slice;
1052  addDaughters(pfpMap, pfp_ptr, pfps_h, particles_in_slice);
1053 
1054  double sliceQ = 0.;
1055  std::vector<art::Ptr<recob::Hit>> hits_in_slice;
1056  flashmatch::QCluster_t particlesClusters;
1057  for(const auto& particle: particles_in_slice) {
1058  const auto particle_key = particle.key();
1059  const auto& particle_spacepoints = pfp_spacepoints_assns.at(particle_key);
1060  double particleQ = 0.;
1061  flashmatch::QCluster_t spsClusters;
1062  for(const auto& spacepoint : particle_spacepoints) {
1063  const auto spacepoint_key = spacepoint.key();
1064  const auto& hits = spacepoint_hits_assns.at(spacepoint_key);
1065  if(fStoreCheatMCT0){
1066  hits_in_slice.insert(hits_in_slice.end(),
1067  hits.begin(), hits.end());
1068  }
1069  const auto& pos = spacepoint->XYZ();
1070  double spacepointQ = 0.;
1071  flashmatch::QCluster_t hitsClusters;
1072  for(const auto& hit : hits) {
1073  geo::WireID wId = hit->WireID();
1074  if(fOnlyCollectionWires &&
1075  fGeometry->SignalType(wId) != geo::kCollection) continue;
1076  const double hitQ = hit->Integral();
1077  if(hitQ < fMinHitQ) continue;
1078  spacepointQ += hitQ;
1079  hitsClusters.emplace_back(pos[0], pos[1], pos[2], hitQ);
1080  const auto itpc = wId.TPC;
1081  tpcWithHits.insert(itpc);
1082  } // for hits associated to spacepoint
1083  if(spacepointQ < fMinSpacePointQ) continue;
1084  particleQ += spacepointQ;
1085  spsClusters.insert(spsClusters.end(),
1086  hitsClusters.begin(), hitsClusters.end());
1087  } // for spacepoints in particle
1088  double chargeToNPhots = lar_pandora::LArPandoraHelper::IsTrack(particle) ?
1090  particleQ *= chargeToNPhots;
1091  if(particleQ < fMinParticleQ) continue;
1092  sliceQ += particleQ;
1093  particlesClusters.insert(particlesClusters.end(),
1094  spsClusters.begin(), spsClusters.end());
1095  } // for particles in slice
1096  if(sliceQ < fMinSliceQ) continue;
1097  double mcT0 = (fStoreCheatMCT0) ? cheatMCT0(hits_in_slice, mcParticles)/1000.
1098  : -9999.;
1099  chargeDigestMap[sliceQ] =
1100  ChargeDigest(pId, pfpPDGC, pfp_ptr, particlesClusters, tpcWithHits, mcT0);
1101  } // over all slices
1102  return chargeDigestMap;
1103 }
double cheatMCT0(const std::vector< art::Ptr< recob::Hit >> &hits, const std::vector< art::Ptr< simb::MCParticle >> &mcParticles)
const double fMinSliceQ
const art::ServiceHandle< geo::Geometry > fGeometry
const art::InputTag fPandoraProducer
process_name hit
Definition: cheaterreco.fcl:51
T abs(T value)
const double fMinParticleQ
const double fMinHitQ
Collection of charge deposition 3D point (cluster)
const bool fSelectNeutrino
const double fMinSpacePointQ
const art::InputTag fSpacePointProducer
std::map< double, ChargeDigest, std::greater< double >> ChargeDigestMap
const double fChargeToNPhotonsTrack
static bool IsTrack(const art::Ptr< recob::PFParticle > particle)
Determine whether a particle has been reconstructed as track-like.
const bool fOnlyCollectionWires
const bool fStoreCheatMCT0
TCEvent evt
Definition: DataStructs.cxx:8
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
static void CollectMCParticles(const art::Event &evt, const std::string &label, MCParticleVector &particleVector)
Collect a vector of MCParticle objects from the ART event record.
void addDaughters(const std::unordered_map< size_t, size_t > &pfpMap, const art::Ptr< recob::PFParticle > &pfp_ptr, const art::ValidHandle< std::vector< recob::PFParticle >> &pfps_h, std::vector< art::Ptr< recob::PFParticle >> &pfp_v) const
const double fChargeToNPhotonsShower
Signal from collection planes.
Definition: geo_types.h:146
std::vector< FlashPredict::SimpleFlash > FlashPredict::makeSimpleFlashes ( std::vector< recob::OpHit > &  opHits,
std::vector< recob::OpHit > &  opHitsRght,
std::vector< recob::OpHit > &  opHitsLeft 
) const
private

Definition at line 1220 of file FlashPredict_module.cc.

1224 {
1225  std::unique_ptr<TH1D> opHitsTimeHist = std::make_unique<TH1D>(
1226  "opHitsTimeHist", "ophittime", fTimeBins, fBeamWindowStart, fBeamWindowEnd);
1227  opHitsTimeHist->SetOption("HIST");
1228  opHitsTimeHist->SetDirectory(0);//turn off ROOT's object ownership
1229  std::unique_ptr<TH1D> opHitsTimeHistRght = std::make_unique<TH1D>(
1230  "opHitsTimeHistRght", "ophittimer", fTimeBins, fBeamWindowStart, fBeamWindowEnd);
1231  opHitsTimeHistRght->SetOption("HIST");
1232  opHitsTimeHistRght->SetDirectory(0);//turn off ROOT's object ownership
1233  std::unique_ptr<TH1D> opHitsTimeHistLeft = std::make_unique<TH1D>(
1234  "opHitsTimeHistLeft", "ophittimel", fTimeBins, fBeamWindowStart, fBeamWindowEnd);
1235  opHitsTimeHistLeft->SetOption("HIST");
1236  opHitsTimeHistLeft->SetDirectory(0);//turn off ROOT's object ownership
1238  opHits, opHitsRght, opHitsLeft,
1239  opHitsTimeHist, opHitsTimeHistRght, opHitsTimeHistLeft)) return {};
1240 
1241  bool oph_in_rght = false, oph_in_left = false;
1242  std::vector<FlashPredict::SimpleFlash> simpleFlashes;
1243  if(opHitsRght.size() > 0 && opHitsTimeHistRght->GetEntries() > 0){
1244  oph_in_rght = true;
1245  findSimpleFlashes(simpleFlashes, opHitsRght,
1246  kActivityInRght, opHitsTimeHistRght);
1247  }
1248  if(opHitsLeft.size() > 0 && opHitsTimeHistLeft->GetEntries() > 0){
1249  oph_in_left = true;
1250  findSimpleFlashes(simpleFlashes, opHitsLeft,
1251  kActivityInLeft, opHitsTimeHistLeft);
1252  }
1253  if(oph_in_rght && oph_in_left ){
1254  findSimpleFlashes(simpleFlashes, opHits, kActivityInBoth, opHitsTimeHist);
1255  }
1256  else if(!oph_in_rght && !oph_in_left){
1257  return {};
1258  }
1259  return simpleFlashes;
1260 }
static constexpr unsigned kActivityInRght
const double fBeamWindowStart
static constexpr unsigned kActivityInBoth
const double fBeamWindowEnd
const unsigned fTimeBins
bool findSimpleFlashes(std::vector< SimpleFlash > &simpleFlashes, std::vector< recob::OpHit > &opHits, const unsigned ophsInVolume, std::unique_ptr< TH1D > &opHitsTimeHist) const
static constexpr unsigned kActivityInLeft
bool createOpHitsTimeHist(const std::vector< recob::OpHit > &opHits, std::vector< recob::OpHit > &opHitsRght, std::vector< recob::OpHit > &opHitsLeft, std::unique_ptr< TH1D > &opHitsTimeHist, std::unique_ptr< TH1D > &opHitsTimeHistRght, std::unique_ptr< TH1D > &opHitsTimeHistLeft) const
std::vector< FlashPredict::SimpleFlash > FlashPredict::makeSimpleFlashes ( std::vector< recob::OpHit > &  opHits) const
private

Definition at line 1264 of file FlashPredict_module.cc.

1266 {
1267  std::unique_ptr<TH1D> opHitsTimeHist = std::make_unique<TH1D>(
1268  "opHitsTimeHist", "ophittime", fTimeBins, fBeamWindowStart, fBeamWindowEnd);
1269  opHitsTimeHist->SetOption("HIST");
1270  opHitsTimeHist->SetDirectory(0);//turn off ROOT's object ownership
1271  unsigned ophsInVolume = createOpHitsTimeHist(opHits, opHitsTimeHist);
1272  if(ophsInVolume == 0) return {};
1273 
1274  std::vector<FlashPredict::SimpleFlash> simpleFlashes;
1275  if(!findSimpleFlashes(simpleFlashes, opHits, ophsInVolume, opHitsTimeHist))
1276  return {};
1277  return simpleFlashes;
1278 }
const double fBeamWindowStart
const double fBeamWindowEnd
const unsigned fTimeBins
bool findSimpleFlashes(std::vector< SimpleFlash > &simpleFlashes, std::vector< recob::OpHit > &opHits, const unsigned ophsInVolume, std::unique_ptr< TH1D > &opHitsTimeHist) const
bool createOpHitsTimeHist(const std::vector< recob::OpHit > &opHits, std::vector< recob::OpHit > &opHitsRght, std::vector< recob::OpHit > &opHitsLeft, std::unique_ptr< TH1D > &opHitsTimeHist, std::unique_ptr< TH1D > &opHitsTimeHistRght, std::unique_ptr< TH1D > &opHitsTimeHistLeft) const
FlashPredict& FlashPredict::operator= ( FlashPredict const &  )
delete
FlashPredict& FlashPredict::operator= ( FlashPredict &&  )
delete
bool FlashPredict::pfpNeutrinoOnEvent ( const art::ValidHandle< std::vector< recob::PFParticle >> &  pfps_h) const
inlineprivate

Definition at line 1184 of file FlashPredict_module.cc.

1186 {
1187  for (auto const& p : (*pfps_h)) {
1188  unsigned pfpPDGC = std::abs(p.PdgCode());
1189  if ((pfpPDGC == 12) ||
1190  (pfpPDGC == 14) ||
1191  (pfpPDGC == 16)) {
1192  return true;
1193  }
1194  }
1195  return false;
1196 }
pdgs p
Definition: selectors.fcl:22
T abs(T value)
template<typename Stream >
void FlashPredict::printBookKeeping ( Stream &&  out)
private

Definition at line 1537 of file FlashPredict_module.cc.

1538 {
1539  std::ostringstream m;
1540  m << "Bookkeeping\n";
1541  m << "----------------------------------------\n"
1542  << "Job Tally\n"
1543  << "\tEvents: \t " << bk.events << "\n";
1544  if(bk.nopfpneutrino) {
1545  m << "\tNo PFP Neutrino: \t -" << bk.nopfpneutrino
1546  << ", scored as: " << kNoPFPInEvt << "\n";
1547  }
1548  if(bk.noslice) {
1549  m << "\tNo Slice: \t -" << bk.noslice
1550  << ", scored as: " << kNoSlcInEvt << "\n";
1551  }
1552  if(bk.nonvalidophit) {
1553  m << "\tNon Valid OpHits: \t -" << bk.nonvalidophit
1554  << ", scored as: " << kNoOpHInEvt << "\n";
1555  }
1556  if(bk.nullophittime) {
1557  m << "\tNo OpHits in-time:\t -" << bk.nullophittime
1558  << ", scored as: " << kNoOpHInEvt << "\n";
1559  }
1560  m << "\t----------------------\n";
1562  m << "\tJob Bookkeeping: \t" << bk.job_bookkeeping << " ERROR!\n";
1563  m << "\tEvents Processed: \t" << bk.events_processed << "\n"
1564  << "----------------------------------------\n"
1565  << "PFP Tally\n"
1566  << "\tPFP to Score: \t " << bk.pfp_to_score << "\n";
1567  if(bk.no_oph_hits) {
1568  m << "\tHits w/o OpHits:\t -" << bk.no_oph_hits
1569  << ", scored as: " << kQNoOpHScr << "\n";
1570  }
1571  if(bk.no_charge) {
1572  m << "\tNo Charge: \t -" << bk.no_charge
1573  << ", scored as: " << kNoChrgScr << "\n";
1574  }
1575  if(bk.no_flash_pe) {
1576  m << "\tNo VUV PE: \t -" << bk.no_flash_pe
1577  << ", scored as: " << k0VUVPEScr << " ERROR!\n";
1578  }
1579  if(bk.no_nu_candidate) {
1580  m << "\tNo Nu Candidate:\t -" << bk.no_nu_candidate
1581  << ", scored as: " << kNotANuScr << "\n";
1582  }
1583  m << "\t----------------------\n";
1585  m << "\tPFP Bookkeeping: \t" << bk.pfp_bookkeeping << " ERROR!\n";
1586  m << "\tScored PFP: \t" << bk.scored_pfp << "\n"
1587  << "----------------------------------------";
1588  out << m.str();
1589 }
static constexpr int kNoPFPInEvt
static constexpr int kNotANuScr
static constexpr int kNoChrgScr
tuple m
now if test mode generate materials, CRT shell, world, gdml header else just generate CRT shell for u...
BookKeeping bk
static constexpr int kQNoOpHScr
static constexpr int kNoSlcInEvt
static constexpr int kNoOpHInEvt
static constexpr int k0VUVPEScr
template<typename Stream >
void FlashPredict::printMetrics ( const std::string  metric,
const ChargeMetrics charge,
const FlashMetrics flash,
const int  pdgc,
const std::set< unsigned > &  tpcWithHits,
const double  term,
Stream &&  out 
) const
private

Definition at line 1611 of file FlashPredict_module.cc.

1618 {
1619  int xbin = static_cast<int>(fXBins * (_charge_x / fDriftDistance));
1620  std::string tpcs;
1621  for(auto itpc: tpcWithHits) tpcs += std::to_string(itpc) + ' ';
1622  out
1623  << "Big term " << metric << ":\t" << term << "\n"
1624  << std::left << std::setw(12) << std::setfill(' ')
1625  << "xbin: \t" << xbin << "\n"
1626  << "pfp.PdgCode:\t" << pdgc << "\n"
1627  << "tpcWithHits:\t" << tpcs << "\n"
1628  << "_slices: \t" << std::setw(8) << _slices << "\n"
1629  << "_true_nus: \t" << std::setw(8) << _true_nus << "\n"
1630  << "_mcT0: \t" << std::setw(8) << _mcT0 << "\n"
1631  << "charge metrics:\n" << charge.dumpMetrics() << "\n"
1632  << "flash metrics:\n" << flash.dumpMetrics() << "\n";
1633 }
const int fXBins
double _charge_x
walls no left
Definition: selectors.fcl:105
const double fDriftDistance
std::string to_string(WindowPattern const &pattern)
unsigned _slices
unsigned _true_nus
void FlashPredict::produce ( art::Event &  evt)
override

Definition at line 139 of file FlashPredict_module.cc.

140 {
141  // sFM is an alias for sbn::SimpleFlashMatch
142  std::unique_ptr< std::vector<sFM> >
143  sFM_v(new std::vector<sFM>);
144  std::unique_ptr< art::Assns <recob::PFParticle, sFM> >
145  pfp_sFM_assn_v(new art::Assns<recob::PFParticle, sFM>);
146 
147  // reset TTree variables
148  if(fMakeTree){
149  _evt = evt.event();
150  _sub = evt.subRun();
151  _run = evt.run();
152  _slices = 0;
153  _true_nus = -9999;
154  _flash_time = -9999.;
155  _flash_pe = -9999.;
156  _flash_unpe = -9999.;
157  _flash_rr = -9999.;
158  _flash_ratio = -9999.;
159  _hypo_x = -9999.;
160  _hypo_x_err = -9999.;
161  _hypo_x_rr = -9999.;
162  _hypo_x_ratio = -9999.;
163  _score = -9999.;
164  _mcT0 = -9999.;
165  }
166  bk.events++;
167 
168  if(fMakeTree && fStoreTrueNus){
169  _true_nus = trueNus(evt);
170  }
171 
172  // grab PFParticles in event
173  const auto pfps_h =
174  evt.getValidHandle<std::vector<recob::PFParticle>>(fPandoraProducer);
175  if (fSelectNeutrino &&
176  !pfpNeutrinoOnEvent(pfps_h)) {
177  mf::LogInfo("FlashPredict")
178  << "No pfp neutrino on event. Skipping...";
179  bk.nopfpneutrino++;
181  for(size_t pId=0; pId<pfps_h->size(); pId++) {
182  if(!pfps_h->at(pId).IsPrimary()) continue;
183  const art::Ptr<recob::PFParticle> pfp_ptr(pfps_h, pId);
184  sFM_v->emplace_back(sFM(kNoScr, kNoScrTime, Charge(kNoScrQ),
186  util::CreateAssn(*this, evt, *sFM_v, pfp_ptr, *pfp_sFM_assn_v);
187  }
188  evt.put(std::move(sFM_v));
189  evt.put(std::move(pfp_sFM_assn_v));
190  return;
191  }
192 
193  if(fMakeTree){
194  auto const& slice_h = evt.getValidHandle<std::vector<recob::Slice>>(fPandoraProducer);
195  _slices = slice_h.product()->size();
196  if (_slices == 0) {
197  mf::LogWarning("FlashPredict")
198  << "No recob:Slice on event. Skipping...";
199  bk.noslice++;
201  for(size_t pId=0; pId<pfps_h->size(); pId++) {
202  if(!pfps_h->at(pId).IsPrimary()) continue;
203  const art::Ptr<recob::PFParticle> pfp_ptr(pfps_h, pId);
204  sFM_v->emplace_back(sFM(kNoScr, kNoScrTime, Charge(kNoScrQ),
206  util::CreateAssn(*this, evt, *sFM_v, pfp_ptr, *pfp_sFM_assn_v);
207  }
208  evt.put(std::move(sFM_v));
209  evt.put(std::move(pfp_sFM_assn_v));
210  return;
211  }
212  }
213 
214  // load OpHits previously created
215  art::Handle<std::vector<recob::OpHit>> ophits_h;
216  evt.getByLabel(fOpHitProducer, ophits_h);
217  if(!ophits_h.isValid()) {
218  mf::LogError("FlashPredict")
219  << "No optical hits from producer module "
220  << fOpHitProducer;
221  bk.nonvalidophit++;
223  for(size_t pId=0; pId<pfps_h->size(); pId++) {
224  if(!pfps_h->at(pId).IsPrimary()) continue;
225  const art::Ptr<recob::PFParticle> pfp_ptr(pfps_h, pId);
226  sFM_v->emplace_back(sFM(kNoScr, kNoScrTime, Charge(kNoScrQ),
228  util::CreateAssn(*this, evt, *sFM_v, pfp_ptr, *pfp_sFM_assn_v);
229  }
230  evt.put(std::move(sFM_v));
231  evt.put(std::move(pfp_sFM_assn_v));
232  return;
233  }
234 
235  std::vector<recob::OpHit> opHits(ophits_h->size());
236  copyOpHitsInBeamWindow(opHits, ophits_h);
237 
238  if(fUseARAPUCAS && !fOpHitARAProducer.empty()){
239  art::Handle<std::vector<recob::OpHit>> ophits_ara_h;
240  evt.getByLabel(fOpHitARAProducer, ophits_ara_h);
241  if(!ophits_ara_h.isValid()) {
242  mf::LogError("FlashPredict")
243  << "Non valid ophits from ARAPUCAS"
244  << "\nfUseARAPUCAS: " << std::boolalpha << fUseARAPUCAS
245  << "\nfOpHitARAProducer: " << fOpHitARAProducer;
246  }
247  else{
248  std::vector<recob::OpHit> opHitsARA(ophits_ara_h->size());
249  copyOpHitsInBeamWindow(opHitsARA, ophits_ara_h);
250  opHits.insert(opHits.end(),
251  opHitsARA.begin(), opHitsARA.end());
252  }
253  }
254 
255 
256  std::vector<recob::OpHit> opHitsRght, opHitsLeft;
257  const std::vector<SimpleFlash> simpleFlashes = (fSBND) ?
258  makeSimpleFlashes(opHits, opHitsRght, opHitsLeft) : makeSimpleFlashes(opHits);
259  if(simpleFlashes.empty()){
260  mf::LogWarning("FlashPredict")
261  << "\nNo OpHits in beam window,"
262  << "\nor the integral is less or equal to " << fMinFlashPE << " or 0."
263  << "\nSkipping...";
264  bk.nullophittime++;
266  for(size_t pId=0; pId<pfps_h->size(); pId++) {
267  if(!pfps_h->at(pId).IsPrimary()) continue;
268  const art::Ptr<recob::PFParticle> pfp_ptr(pfps_h, pId);
269  sFM_v->emplace_back(sFM(kNoScr, kNoScrTime, Charge(kNoScrQ),
271  util::CreateAssn(*this, evt, *sFM_v, pfp_ptr, *pfp_sFM_assn_v);
272  }
273  evt.put(std::move(sFM_v));
274  evt.put(std::move(pfp_sFM_assn_v));
275  return;
276  }
277 
278  ChargeDigestMap chargeDigestMap = makeChargeDigest(evt, pfps_h);
279 
280  std::map<unsigned, FlashMetrics> flashMetricsMap;
281  for(auto& chargeDigest : chargeDigestMap) {
282  //const size_t pId = chargeDigest.second.pId;
283  const int pfpPDGC = chargeDigest.second.pfpPDGC;
284  const auto& pfp_ptr = chargeDigest.second.pfp_ptr;
285  const auto& qClusters = chargeDigest.second.qClusters;
286  const auto& tpcWithHits = chargeDigest.second.tpcWithHits;
287  bk.pfp_to_score++;
288  if(chargeDigest.first < 0.){
289  mf::LogDebug("FlashPredict") << "Not a nu candidate slice. Skipping...";
291  mf::LogDebug("FlashPredict") << "Creating sFM and PFP-sFM association";
292  sFM_v->emplace_back(sFM(kNoScr, kNoScrTime, Charge(kNoScrQ),
294  util::CreateAssn(*this, evt, *sFM_v, pfp_ptr, *pfp_sFM_assn_v);
295  continue;
296  }
297 
298  unsigned hitsInVolume = 0;
299  bool in_right = false, in_left = false;
300  for(unsigned t : tpcWithHits){
301  if(t/fTPCPerDriftVolume == kRght) in_right = true;
302  else if(t/fTPCPerDriftVolume == kLeft) in_left = true;
303  }
304  if(in_right && in_left) hitsInVolume = kActivityInBoth;
305  else if(in_right && !in_left) hitsInVolume = kActivityInRght;
306  else if(!in_right && in_left) hitsInVolume = kActivityInLeft;
307  else {
308  mf::LogError("FlashPredict")
309  << "ERROR!!! tpcWithHits.size() " << tpcWithHits.size();
310  }
311 
312  ChargeMetrics charge = computeChargeMetrics(qClusters);
313  if(!charge.metric_ok){
314  mf::LogWarning("FlashPredict") << "Clusters with No Charge. Skipping...";
315  bk.no_charge++;
316  mf::LogDebug("FlashPredict") << "Creating sFM and PFP-sFM association";
317  sFM_v->emplace_back(sFM(kNoScr, kNoScrTime, Charge(kNoScrQ),
319  util::CreateAssn(*this, evt, *sFM_v, pfp_ptr, *pfp_sFM_assn_v);
320  continue;
321  }
322 
323  FlashMetrics flash = {};
324  Score score = {std::numeric_limits<double>::max()};
325  bool hits_ophits_concurrence = false;
326  for(auto& simpleFlash : simpleFlashes) {
327  unsigned ophsInVolume = simpleFlash.ophsInVolume;
328  if(hitsInVolume != ophsInVolume){
329  if(fSBND){
330  if(fForceConcurrence) continue;
331  else if((hitsInVolume < kActivityInBoth) &&
332  (ophsInVolume < kActivityInBoth)) {
333  continue;
334  }
335  }
336  else if(fICARUS){
337  if((hitsInVolume < kActivityInBoth) &&
338  (ophsInVolume < kActivityInBoth)) {
339  continue;
340  }
341  else if(fForceConcurrence && hitsInVolume == kActivityInBoth) continue;
342  }
343  }
344  hits_ophits_concurrence = true;
345 
346  unsigned flashUId = simpleFlash.ophsInVolume * 10 + simpleFlash.flashId;
347  bool mets_in_map = flashMetricsMap.find(flashUId) != flashMetricsMap.end();
348  FlashMetrics flash_tmp = (mets_in_map) ?
349  flashMetricsMap[flashUId] : computeFlashMetrics(simpleFlash);
350  if(mets_in_map){
351  mf::LogDebug("FlashPredict")
352  << "Reusing metrics previously computed, for flashUId " << flashUId;
353  }
354  else
355  flashMetricsMap[flashUId] = flash_tmp;
356 
357  Score score_tmp = computeScore(charge, flash_tmp, tpcWithHits, pfpPDGC);
358  if(score_tmp.total > 0. && score_tmp.total < score.total
359  && flash_tmp.metric_ok){
360  score = score_tmp;
361  flash = flash_tmp;
362  }
363  } // for simpleFlashes
364  if(!hits_ophits_concurrence) {
365  std::string extra_message = (!fForceConcurrence) ? "" :
366  "\nConsider setting ForceConcurrence to false to lower requirements";
367  mf::LogInfo("FlashPredict")
368  << "No OpHits where there's charge. Skipping..." << extra_message;
369  bk.no_oph_hits++;
370  mf::LogDebug("FlashPredict") << "Creating sFM and PFP-sFM association";
371  sFM_v->emplace_back(sFM(kNoScr, kNoScrTime, Charge(kNoScrQ),
373  util::CreateAssn(*this, evt, *sFM_v, pfp_ptr, *pfp_sFM_assn_v);
374  continue;
375  }
376  else if(!flash.metric_ok){
377  printMetrics("ERROR", charge, flash, pfpPDGC, tpcWithHits, 0, mf::LogError("FlashPredict"));
378  bk.no_flash_pe++;
379  mf::LogDebug("FlashPredict") << "Creating sFM and PFP-sFM association";
380  sFM_v->emplace_back(sFM(kNoScr, kNoScrTime, Charge(kNoScrQ),
382  util::CreateAssn(*this, evt, *sFM_v, pfp_ptr, *pfp_sFM_assn_v);
383  continue;
384  }
385 
386  if(score.total > 0. &&
387  score.total < std::numeric_limits<double>::max()){
388  if(fMakeTree) {
389  _mcT0 = chargeDigest.second.mcT0;
390  updateChargeMetrics(charge);
391  updateFlashMetrics(flash);
392  updateScore(score);
393  _flashmatch_nuslice_tree->Fill();
394  }
395  bk.scored_pfp++;
396  mf::LogDebug("FlashPredict") << "Creating sFM and PFP-sFM association";
397  Charge c{charge.q, TVector3(charge.x_gl, charge.y, charge.z)};
398  Flash f{flash.pe, TVector3(flash.x_gl, flash.y, flash.z)};
399  sFM_v->emplace_back(sFM(true, flash.time, c, f, score));
400  util::CreateAssn(*this, evt, *sFM_v, pfp_ptr, *pfp_sFM_assn_v);
401  }
402  else{
403  mf::LogError("FlashPredict") << "ERROR: score <= 0. Dumping info."
404  << "\n_score: " << _score
405  << "\n_scr_y: " << _scr_y
406  << "\n_scr_z: " << _scr_z
407  << "\n_scr_rr: " << _scr_rr
408  << "\n_scr_ratio: " << _scr_ratio;
409  printMetrics("ERROR", charge, flash, pfpPDGC, tpcWithHits, 0, mf::LogError("FlashPredict"));
410  }
411 
412  } // chargeDigestMap: PFparticles that pass criteria
415 
416  evt.put(std::move(sFM_v));
417  evt.put(std::move(pfp_sFM_assn_v));
418 
419 }// end of producer module
const bool fForceConcurrence
const bool fSBND
static constexpr int kNoPFPInEvt
unsigned int event
Definition: DataStructs.h:634
static constexpr unsigned kActivityInRght
static constexpr unsigned kRght
static constexpr int kNotANuScr
void printMetrics(const std::string metric, const ChargeMetrics &charge, const FlashMetrics &flash, const int pdgc, const std::set< unsigned > &tpcWithHits, const double term, Stream &&out) const
double _flash_time
unsigned int run
Definition: DataStructs.h:635
unsigned _sub
void updateChargeMetrics(const ChargeMetrics &chargeMetrics)
static constexpr unsigned kLeft
BEGIN_PROLOG or score(default)}sbnd_crttrackmatchingalg_crID
static constexpr int kNoChrgScr
void copyOpHitsInBeamWindow(std::vector< recob::OpHit > &opHits, const art::Handle< std::vector< recob::OpHit >> &ophit_h) const
const art::InputTag fOpHitARAProducer
const art::InputTag fPandoraProducer
const bool fICARUS
static constexpr unsigned kActivityInBoth
Score computeScore(const ChargeMetrics &charge, const FlashMetrics &flash, const std::set< unsigned > &tpcWithHits, const int pdgc) const
bool pfpNeutrinoOnEvent(const art::ValidHandle< std::vector< recob::PFParticle >> &pfps_h) const
static constexpr double kNoScrQ
double _scr_ratio
const bool fMakeTree
BookKeeping bk
double _hypo_x_err
static constexpr int kQNoOpHScr
sbn::SimpleFlashMatch::Score Score
const bool fUseARAPUCAS
sbn::SimpleFlashMatch::Charge Charge
unsigned trueNus(art::Event &evt) const
double _flash_pe
double _flash_unpe
static constexpr int kNoSlcInEvt
TTree * _flashmatch_nuslice_tree
unsigned _evt
unsigned fTPCPerDriftVolume
static constexpr double kNoScrTime
ChargeDigestMap makeChargeDigest(const art::Event &evt, const art::ValidHandle< std::vector< recob::PFParticle >> &pfps_h)
std::vector< SimpleFlash > makeSimpleFlashes(std::vector< recob::OpHit > &opHits, std::vector< recob::OpHit > &opHitsRght, std::vector< recob::OpHit > &opHitsLeft) const
const bool fSelectNeutrino
static constexpr unsigned kActivityInLeft
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.
sbn::SimpleFlashMatch::Flash Flash
double _flash_ratio
ChargeMetrics computeChargeMetrics(const flashmatch::QCluster_t &qClusters) const
void updateFlashMetrics(const FlashMetrics &flashMetrics)
void updateScore(const Score &score)
const double fMinFlashPE
sbn::SimpleFlashMatch sFM
static constexpr bool kNoScr
double _hypo_x_rr
std::map< double, ChargeDigest, std::greater< double >> ChargeDigestMap
static constexpr double kNoScrPE
const art::InputTag fOpHitProducer
unsigned _slices
unsigned int subRun
Definition: DataStructs.h:636
FlashMetrics computeFlashMetrics(const SimpleFlash &simpleFlash) const
const bool fStoreTrueNus
double pe
photo-electrons on flash
static constexpr int kNoOpHInEvt
TCEvent evt
Definition: DataStructs.cxx:8
static constexpr int k0VUVPEScr
double _flash_rr
double _hypo_x_ratio
unsigned _true_nus
unsigned _run
unsigned FlashPredict::sbndPDinTPC ( const int  pdChannel) const
private

Definition at line 1430 of file FlashPredict_module.cc.

1431 {
1432  auto p = fGeometry->OpDetGeoFromOpChannel(pdChannel).GetCenter();
1433  p.SetX(p.X()/2.);//OpDets are outside the TPCs
1434  return (fGeometry->PositionToTPCID(p)).TPC;
1435 }
const art::ServiceHandle< geo::Geometry > fGeometry
pdgs p
Definition: selectors.fcl:22
BEGIN_PROLOG TPC
double FlashPredict::scoreTerm ( const double  m,
const double  n,
const double  mean,
const double  spread 
) const
inlineprivate

Definition at line 1168 of file FlashPredict_module.cc.

1170 {
1171  return std::abs(std::abs(m - n) - mean) / spread;
1172 }
tuple m
now if test mode generate materials, CRT shell, world, gdml header else just generate CRT shell for u...
T abs(T value)
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:13
double FlashPredict::scoreTerm ( const double  m,
const double  mean,
const double  spread 
) const
inlineprivate

Definition at line 1176 of file FlashPredict_module.cc.

1178 {
1179  return std::abs(m - mean) / spread;
1180 }
tuple m
now if test mode generate materials, CRT shell, world, gdml header else just generate CRT shell for u...
T abs(T value)
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:13
unsigned FlashPredict::trueNus ( art::Event &  evt) const
private

Definition at line 1122 of file FlashPredict_module.cc.

1123 {
1124  art::Handle<std::vector<simb::MCTruth> > mctruthList_h;
1125  std::vector<art::Ptr<simb::MCTruth> > mclist;
1126  if(evt.getByLabel("generator", mctruthList_h))
1127  art::fill_ptr_vector(mclist, mctruthList_h);
1128  unsigned true_nus = 0;
1129  for(auto const& mc: mclist){
1130  if(mc->Origin() == simb::kBeamNeutrino) ++true_nus;
1131  }
1132  return true_nus;
1133 }
TCEvent evt
Definition: DataStructs.cxx:8
void FlashPredict::updateBookKeeping ( )
private

Definition at line 1592 of file FlashPredict_module.cc.

1593 {
1594  // account for the reasons that an event could lack
1598 
1599  // account for the reasons that a particle might lack
1603 
1606  printBookKeeping(mf::LogWarning("FlashPredict"));
1607 }
BookKeeping bk
T abs(T value)
void printBookKeeping(Stream &&out)
void FlashPredict::updateChargeMetrics ( const ChargeMetrics chargeMetrics)
inlineprivate

Definition at line 1137 of file FlashPredict_module.cc.

1138 {
1139  const auto& c = chargeMetrics;
1140  _charge_x = c.x; _charge_x_gl = c.x_gl; _charge_y = c.y;
1141  _charge_z = c.z; _charge_q = c.q;
1142 }
double _charge_q
double _charge_x_gl
double _charge_y
double _charge_x
double _charge_z
void FlashPredict::updateFlashMetrics ( const FlashMetrics flashMetrics)
inlineprivate

Definition at line 1146 of file FlashPredict_module.cc.

1147 {
1148  const auto& f = flashMetrics;
1149  _flash_x = f.x; _flash_x_gl = f.x_gl;
1150  _flash_y = f.y; _flash_yb = f.yb; _flash_z = f.z; _flash_zb = f.zb;
1151  _flash_rr = f.rr; _flash_pe = f.pe; _flash_unpe = f.unpe;
1152  _flash_ratio = f.ratio; _flash_time = f.time;
1153  _hypo_x = f.h_x; _hypo_x_err = f.h_xerr; _hypo_x_rr = f.h_xrr;
1154  _hypo_x_ratio = f.h_xratio;
1155  _y_skew = f.y_skew; _z_skew = f.z_skew;
1156 }
double _flash_time
double _flash_zb
double _hypo_x_err
double _flash_pe
double _flash_unpe
double _flash_x_gl
double _flash_ratio
double _hypo_x_rr
double _flash_rr
double _hypo_x_ratio
double _flash_yb
void FlashPredict::updateScore ( const Score score)
inlineprivate

Definition at line 1160 of file FlashPredict_module.cc.

1161 {
1162  _score = score.total, _scr_y = score.y, _scr_z = score.z,
1163  _scr_rr = score.rr, _scr_ratio = score.ratio;
1164 }
BEGIN_PROLOG or score(default)}sbnd_crttrackmatchingalg_crID
double _scr_ratio
double FlashPredict::wallXWithMaxPE ( const OpHitIt  opH_beg,
const OpHitIt  opH_end 
) const
private

Definition at line 1438 of file FlashPredict_module.cc.

1440 {
1441  std::map<double, double> opdetX_PE {{-99999., 0.}};
1442  for(auto oph=opH_beg; oph!=opH_end; ++oph){
1443  double ophPE = oph->PE();
1444  double ophPE2 = ophPE*ophPE;
1445  double opdetX = fGeometry->OpDetGeoFromOpChannel(
1446  oph->OpChannel()).GetCenter().X();
1447  bool stored = false;
1448  for(auto& m : opdetX_PE){
1449  if(std::abs(m.first - opdetX) < 5.) {
1450  m.second += ophPE2;
1451  stored = true;
1452  break;
1453  }
1454  }
1455  if(!stored) opdetX_PE[opdetX] = ophPE2;
1456  }
1457  auto maxIt = std::max_element(
1458  opdetX_PE.begin(), opdetX_PE.end(),
1459  [] (const auto& a, const auto& b) ->bool{return a.second < b.second;});
1460  return maxIt->first;
1461 }
const art::ServiceHandle< geo::Geometry > fGeometry
tuple m
now if test mode generate materials, CRT shell, world, gdml header else just generate CRT shell for u...
process_name gaushit a
T abs(T value)
std::list< double > FlashPredict::wiresXGl ( ) const
private

Definition at line 1464 of file FlashPredict_module.cc.

1465 {
1466  std::list<double> wiresX_gl;
1467  for (size_t t = 0; t < fNTPC; t++) {
1468  const geo::TPCGeo& tpcg = fGeoCryo->TPC(t);
1469  wiresX_gl.push_back(tpcg.LastPlane().GetCenter().X());
1470  }
1471  wiresX_gl.unique([](double l, double r) { return std::abs(l - r) < 0.00001;});
1472  return wiresX_gl;
1473 }
Geometry information for a single TPC.
Definition: TPCGeo.h:38
T abs(T value)
Point GetCenter() const
Returns the centre of the wire plane in world coordinates [cm].
Definition: PlaneGeo.h:479
geo::PlaneGeo const & LastPlane() const
Returns the last wire plane (the farther from TPC center).
Definition: TPCGeo.h:248
const size_t fNTPC
const std::unique_ptr< geo::CryostatGeo > fGeoCryo
esac echo uname r
std::tuple< double, double > FlashPredict::xEstimateAndRMS ( double  metric_value,
const TH2D *  metric_h2 
) const
private

Definition at line 941 of file FlashPredict_module.cc.

943 {
944  int bin = metric_h2->GetYaxis()->FindBin(metric_value);
945  int bins = metric_h2->GetNbinsY();
946  double metric_hypoX = -1.;
947  double metric_hypoXWgt = 0.;
948  int bin_buff = 0;
949  while(0 < bin-bin_buff || bin+bin_buff <= bins){
950  int low_bin = (0 < bin-bin_buff) ? bin-bin_buff : 0;
951  int high_bin = (bin+bin_buff <= bins) ? bin+bin_buff : -1;
952  TH1D* metric_px = metric_h2->ProjectionX("metric_px", low_bin, high_bin);
953  if(metric_px->GetEntries() > kMinEntriesInProjection){
954  metric_hypoX = metric_px->GetRandom();
955  // metric_hypoX = metric_px->GetMean(); TODO: which one is more justified?
956  double metric_rmsX = metric_px->GetRMS();
957  if(metric_rmsX < fXBinWidth){//something went wrong
958  mf::LogDebug("FlashPredict")
959  << "metric_h2 projected on metric_value: "<< metric_value
960  << ", bin: " << bin
961  << ", bin_buff: " << bin_buff
962  << "; has " << metric_px->GetEntries() << " entries."
963  << "\nmetric_hypoX: " << metric_hypoX
964  << ", metric_rmsX: " << metric_rmsX;
965  return {-1., 0.}; // no estimate
966  }
967  metric_hypoXWgt = 1/(metric_rmsX*metric_rmsX);
968  return {metric_hypoX, metric_hypoXWgt};
969  }
970  bin_buff += 1;
971  }
972  return {-1., 0.}; // no estimate
973 }
const double fXBinWidth
constexpr details::BinObj< T > bin(T value)
Returns a wrapper to print the specified data in binary format.
static constexpr unsigned kMinEntriesInProjection

Member Data Documentation

double FlashPredict::_charge_q
private

Definition at line 351 of file FlashPredict.hh.

double FlashPredict::_charge_x
private

Definition at line 351 of file FlashPredict.hh.

double FlashPredict::_charge_x_gl
private

Definition at line 351 of file FlashPredict.hh.

double FlashPredict::_charge_y
private

Definition at line 351 of file FlashPredict.hh.

double FlashPredict::_charge_z
private

Definition at line 351 of file FlashPredict.hh.

unsigned FlashPredict::_evt
private

Definition at line 358 of file FlashPredict.hh.

double FlashPredict::_flash_pe
private

Definition at line 353 of file FlashPredict.hh.

double FlashPredict::_flash_ratio
private

Definition at line 353 of file FlashPredict.hh.

double FlashPredict::_flash_rr
private

Definition at line 353 of file FlashPredict.hh.

double FlashPredict::_flash_time
private

Definition at line 353 of file FlashPredict.hh.

double FlashPredict::_flash_unpe
private

Definition at line 353 of file FlashPredict.hh.

double FlashPredict::_flash_x
private

Definition at line 353 of file FlashPredict.hh.

double FlashPredict::_flash_x_gl
private

Definition at line 353 of file FlashPredict.hh.

double FlashPredict::_flash_y
private

Definition at line 353 of file FlashPredict.hh.

double FlashPredict::_flash_yb
private

Definition at line 353 of file FlashPredict.hh.

double FlashPredict::_flash_z
private

Definition at line 353 of file FlashPredict.hh.

double FlashPredict::_flash_zb
private

Definition at line 353 of file FlashPredict.hh.

TTree* FlashPredict::_flashmatch_nuslice_tree
private

Definition at line 338 of file FlashPredict.hh.

double FlashPredict::_hypo_x
private

Definition at line 353 of file FlashPredict.hh.

double FlashPredict::_hypo_x_err
private

Definition at line 353 of file FlashPredict.hh.

double FlashPredict::_hypo_x_ratio
private

Definition at line 353 of file FlashPredict.hh.

double FlashPredict::_hypo_x_rr
private

Definition at line 353 of file FlashPredict.hh.

double FlashPredict::_mcT0 = -9999.
private

Definition at line 360 of file FlashPredict.hh.

unsigned FlashPredict::_run
private

Definition at line 358 of file FlashPredict.hh.

double FlashPredict::_score
private

Definition at line 357 of file FlashPredict.hh.

double FlashPredict::_scr_ratio
private

Definition at line 357 of file FlashPredict.hh.

double FlashPredict::_scr_rr
private

Definition at line 357 of file FlashPredict.hh.

double FlashPredict::_scr_y
private

Definition at line 357 of file FlashPredict.hh.

double FlashPredict::_scr_z
private

Definition at line 357 of file FlashPredict.hh.

unsigned FlashPredict::_slices = -1
private

Definition at line 359 of file FlashPredict.hh.

unsigned FlashPredict::_sub
private

Definition at line 358 of file FlashPredict.hh.

unsigned FlashPredict::_true_nus = -1
private

Definition at line 359 of file FlashPredict.hh.

double FlashPredict::_y_skew
private

Definition at line 353 of file FlashPredict.hh.

double FlashPredict::_z_skew
private

Definition at line 353 of file FlashPredict.hh.

BookKeeping FlashPredict::bk
private

Definition at line 385 of file FlashPredict.hh.

const double FlashPredict::fBeamWindowEnd
private

Definition at line 293 of file FlashPredict.hh.

const double FlashPredict::fBeamWindowStart
private

Definition at line 293 of file FlashPredict.hh.

const double FlashPredict::fChargeToNPhotonsShower
private

Definition at line 305 of file FlashPredict.hh.

const double FlashPredict::fChargeToNPhotonsTrack
private

Definition at line 305 of file FlashPredict.hh.

detinfo::DetectorClocksData const FlashPredict::fClockData
private

Definition at line 291 of file FlashPredict.hh.

const int FlashPredict::fCryostat
private

Definition at line 314 of file FlashPredict.hh.

const std::string FlashPredict::fDetector
private

Definition at line 310 of file FlashPredict.hh.

const double FlashPredict::fDriftDistance
private

Definition at line 318 of file FlashPredict.hh.

unsigned FlashPredict::fDriftVolumes
private

Definition at line 325 of file FlashPredict.hh.

std::vector<double> FlashPredict::fdYMeans
private

Definition at line 362 of file FlashPredict.hh.

std::vector<double> FlashPredict::fdYSpreads
private

Definition at line 363 of file FlashPredict.hh.

std::vector<double> FlashPredict::fdZMeans
private

Definition at line 362 of file FlashPredict.hh.

std::vector<double> FlashPredict::fdZSpreads
private

Definition at line 363 of file FlashPredict.hh.

const double FlashPredict::fFlashEnd
private

Definition at line 294 of file FlashPredict.hh.

const double FlashPredict::fFlashStart
private

Definition at line 294 of file FlashPredict.hh.

const bool FlashPredict::fForceConcurrence
private

Definition at line 298 of file FlashPredict.hh.

const std::unique_ptr<geo::CryostatGeo> FlashPredict::fGeoCryo
private

Definition at line 316 of file FlashPredict.hh.

const art::ServiceHandle<geo::Geometry> FlashPredict::fGeometry
private

Definition at line 309 of file FlashPredict.hh.

const bool FlashPredict::fICARUS
private

Definition at line 311 of file FlashPredict.hh.

const std::string FlashPredict::fInputFilename
private

Definition at line 303 of file FlashPredict.hh.

const bool FlashPredict::fMakeTree
private

Definition at line 304 of file FlashPredict.hh.

const unsigned FlashPredict::fMaxFlashes
private

Definition at line 307 of file FlashPredict.hh.

const double FlashPredict::fMinFlashPE
private

Definition at line 308 of file FlashPredict.hh.

const double FlashPredict::fMinHitQ
private

Definition at line 306 of file FlashPredict.hh.

const double FlashPredict::fMinOpHPE
private

Definition at line 308 of file FlashPredict.hh.

const double FlashPredict::fMinParticleQ
private

Definition at line 306 of file FlashPredict.hh.

const double FlashPredict::fMinSliceQ
private

Definition at line 306 of file FlashPredict.hh.

const double FlashPredict::fMinSpacePointQ
private

Definition at line 306 of file FlashPredict.hh.

const bool FlashPredict::fNoAvailableMetrics
private

Definition at line 304 of file FlashPredict.hh.

const size_t FlashPredict::fNTPC
private

Definition at line 313 of file FlashPredict.hh.

const bool FlashPredict::fOnlyCollectionWires
private

Definition at line 297 of file FlashPredict.hh.

const unsigned FlashPredict::fOpDetNormalizer
private

Definition at line 327 of file FlashPredict.hh.

const art::InputTag FlashPredict::fOpHitARAProducer
private

Definition at line 289 of file FlashPredict.hh.

const art::InputTag FlashPredict::fOpHitProducer
private

Definition at line 289 of file FlashPredict.hh.

const art::InputTag FlashPredict::fPandoraProducer
private

Definition at line 289 of file FlashPredict.hh.

const std::unique_ptr<opdet::PDMapAlg> FlashPredict::fPDMapAlgPtr
private

Definition at line 312 of file FlashPredict.hh.

const std::string FlashPredict::fRatio_TF1_fit
private

Definition at line 321 of file FlashPredict.hh.

std::array<Fits, 3> FlashPredict::fRatioFits
private

Definition at line 346 of file FlashPredict.hh.

TH2D* FlashPredict::fRatioH2
private

Definition at line 343 of file FlashPredict.hh.

std::vector<double> FlashPredict::fRatioMeans
private

Definition at line 362 of file FlashPredict.hh.

std::vector<double> FlashPredict::fRatioSpreads
private

Definition at line 363 of file FlashPredict.hh.

const std::string FlashPredict::fRR_TF1_fit
private

Definition at line 321 of file FlashPredict.hh.

std::array<Fits, 3> FlashPredict::fRRFits
private

Definition at line 345 of file FlashPredict.hh.

TH2D* FlashPredict::fRRH2
private

Definition at line 343 of file FlashPredict.hh.

std::vector<double> FlashPredict::fRRMeans
private

Definition at line 362 of file FlashPredict.hh.

std::vector<double> FlashPredict::fRRSpreads
private

Definition at line 363 of file FlashPredict.hh.

const bool FlashPredict::fSBND
private

Definition at line 311 of file FlashPredict.hh.

const bool FlashPredict::fSelectNeutrino
private

Definition at line 296 of file FlashPredict.hh.

const art::InputTag FlashPredict::fSpacePointProducer
private

Definition at line 289 of file FlashPredict.hh.

const bool FlashPredict::fStoreCheatMCT0
private

Definition at line 302 of file FlashPredict.hh.

const bool FlashPredict::fStoreTrueNus
private

Definition at line 301 of file FlashPredict.hh.

const double FlashPredict::fTermThreshold
private

Definition at line 328 of file FlashPredict.hh.

const double FlashPredict::fTickPeriod
private

Definition at line 292 of file FlashPredict.hh.

const unsigned FlashPredict::fTimeBins
private

Definition at line 295 of file FlashPredict.hh.

unsigned FlashPredict::fTPCPerDriftVolume
private

Definition at line 326 of file FlashPredict.hh.

const bool FlashPredict::fUseARAPUCAS
private

Definition at line 300 of file FlashPredict.hh.

const bool FlashPredict::fUseOppVolMetric
private

Definition at line 299 of file FlashPredict.hh.

const bool FlashPredict::fUseUncoatedPMT
private

Definition at line 299 of file FlashPredict.hh.

const std::list<double> FlashPredict::fWiresX_gl
private

Definition at line 317 of file FlashPredict.hh.

const int FlashPredict::fXBins
private

Definition at line 319 of file FlashPredict.hh.

const double FlashPredict::fXBinWidth
private

Definition at line 320 of file FlashPredict.hh.

const double FlashPredict::fYBiasSlope
private

Definition at line 324 of file FlashPredict.hh.

const unsigned FlashPredict::fYBins
private

Definition at line 322 of file FlashPredict.hh.

const double FlashPredict::fYHigh
private

Definition at line 323 of file FlashPredict.hh.

const double FlashPredict::fYLow
private

Definition at line 323 of file FlashPredict.hh.

const double FlashPredict::fZBiasSlope
private

Definition at line 324 of file FlashPredict.hh.

const unsigned FlashPredict::fZBins
private

Definition at line 322 of file FlashPredict.hh.

const double FlashPredict::fZHigh
private

Definition at line 323 of file FlashPredict.hh.

const double FlashPredict::fZLow
private

Definition at line 323 of file FlashPredict.hh.

constexpr int FlashPredict::k0VUVPEScr = -3
staticprivate

Definition at line 371 of file FlashPredict.hh.

constexpr unsigned FlashPredict::kActivityInBoth = 300
staticprivate

Definition at line 335 of file FlashPredict.hh.

constexpr unsigned FlashPredict::kActivityInLeft = 200
staticprivate

Definition at line 334 of file FlashPredict.hh.

constexpr unsigned FlashPredict::kActivityInRght = 100
staticprivate

Definition at line 333 of file FlashPredict.hh.

constexpr double FlashPredict::kEps = 1e-4
staticprivate

Definition at line 348 of file FlashPredict.hh.

constexpr unsigned FlashPredict::kLeft = 1
staticprivate

Definition at line 331 of file FlashPredict.hh.

constexpr unsigned FlashPredict::kMinEntriesInProjection = 100
staticprivate

Definition at line 344 of file FlashPredict.hh.

constexpr int FlashPredict::kNoChrgScr = -2
staticprivate

Definition at line 370 of file FlashPredict.hh.

constexpr int FlashPredict::kNoOpHInEvt = -11
staticprivate

Definition at line 373 of file FlashPredict.hh.

constexpr int FlashPredict::kNoPFPInEvt = -12
staticprivate

Definition at line 374 of file FlashPredict.hh.

constexpr bool FlashPredict::kNoScr = false
staticprivate

Definition at line 365 of file FlashPredict.hh.

constexpr double FlashPredict::kNoScrPE = -9999.
staticprivate

Definition at line 368 of file FlashPredict.hh.

constexpr double FlashPredict::kNoScrQ = -9999.
staticprivate

Definition at line 367 of file FlashPredict.hh.

constexpr double FlashPredict::kNoScrTime = -9999.
staticprivate

Definition at line 366 of file FlashPredict.hh.

constexpr int FlashPredict::kNoSlcInEvt = -13
staticprivate

Definition at line 375 of file FlashPredict.hh.

constexpr int FlashPredict::kNotANuScr = -4
staticprivate

Definition at line 372 of file FlashPredict.hh.

constexpr int FlashPredict::kQNoOpHScr = -1
staticprivate

Definition at line 369 of file FlashPredict.hh.

constexpr unsigned FlashPredict::kRght = 0
staticprivate

Definition at line 330 of file FlashPredict.hh.

const std::array<std::string, 3> FlashPredict::kSuffixes {"l", "h", "m"}
private

Definition at line 347 of file FlashPredict.hh.


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