10 #include "tensorflow/core/public/session.h" 22 #include "larevt/CalibrationDBI/Interface/ChannelStatusService.h" 23 #include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h" 27 #include "CLHEP/Random/RandGauss.h" 37 if ((samples == 0) || inps.empty() || inps.front().empty() || inps.front().front().empty())
40 if ((samples == -1) || (samples > (int)inps.size())) { samples = inps.size(); }
42 std::vector< std::vector<float> > results;
43 for (
int i = 0; i < samples; ++i)
45 results.push_back(
Run(inps[i]));
53 std::string fname_out;
54 cet::search_path sp(
"FW_SEARCH_PATH");
55 if (!sp.find_file(fileName, fname_out))
58 if (stat(fileName, &buffer) == 0) { fname_out = fileName; }
62 <<
"Could not find the model file " << fileName;
76 mf::LogInfo(
"KerasModelInterface") <<
"Keras model loaded.";
82 std::vector< std::vector< std::vector<float> > > inp3d;
83 inp3d.push_back(inp2d);
102 mf::LogInfo(
"TfModelInterface") <<
"TF model loaded.";
108 if ((samples == 0) || inps.empty() || inps.front().empty() || inps.front().front().empty())
111 if ((samples == -1) || (samples > (
long long int)inps.size())) { samples = inps.size(); }
113 long long int rows = inps.front().size(), cols = inps.front().front().size();
115 tensorflow::Tensor _x(tensorflow::DT_FLOAT, tensorflow::TensorShape({ samples, rows, cols, 1 }));
116 auto input_map = _x.tensor<float, 4>();
117 for (
long long int s = 0;
s < samples; ++
s) {
118 const auto & sample = inps[
s];
119 for (
long long int r = 0; r < rows; ++r) {
120 const auto & row = sample[r];
121 for (
long long int c = 0; c < cols; ++c) {
122 input_map(
s, r, c, 0) = row[c];
133 long long int rows = inp2d.size(), cols = inp2d.front().size();
135 tensorflow::Tensor _x(tensorflow::DT_FLOAT, tensorflow::TensorShape({ 1, rows, cols, 1 }));
136 auto input_map = _x.tensor<float, 4>();
137 for (
long long int r = 0; r < rows; ++r) {
138 const auto & row = inp2d[r];
139 for (
long long int c = 0; c < cols; ++c) {
140 input_map(0, r, c, 0) = row[c];
144 auto out = g->run(_x);
145 if (!out.empty())
return out.front();
146 else return std::vector<float>();
156 fPatchSizeW(config.PatchSizeW()), fPatchSizeD(config.PatchSizeD()),
157 fCurrentWireIdx(99999), fCurrentScaledDrift(99999)
176 mf::LogError(
"PointIdAlg") <<
"File name extension not supported.";
179 if (!
fNNet) {
throw cet::exception(
"nnet::PointIdAlg") <<
"Loading model from file failed."; }
204 mf::LogError(
"PointIdAlg") <<
"Patch buffering failed.";
211 if (!out.empty()) { result = out[outIdx]; }
214 mf::LogError(
"PointIdAlg") <<
"Problem with applying model to input.";
224 std::vector<float> result;
228 mf::LogError(
"PointIdAlg") <<
"Patch buffering failed.";
237 mf::LogError(
"PointIdAlg") <<
"Problem with applying model to input.";
247 if (points.empty() || !
fNNet) {
return std::vector< std::vector<float> >(); }
249 std::vector< std::vector< std::vector<float> > > inps(
250 points.size(), std::vector< std::vector<float> >(
252 for (
size_t i = 0; i < points.size(); ++i)
254 unsigned int wire = points[i].first;
255 float drift = points[i].second;
258 throw cet::exception(
"PointIdAlg") <<
"Patch buffering failed" << std::endl;
273 if ((wire1 == wire2) && (sd1 == sd2))
278 if ((wire1 == wire2) && ((size_t)drift1 == (
size_t)drift2))
305 std::vector<float> flat;
306 if (patch.empty() || patch.front().empty())
312 flat.resize(patch.size() * patch.front().size());
314 for (
size_t w = 0, i = 0;
w < patch.size(); ++
w)
316 auto const & wire = patch[
w];
317 for (
size_t d = 0;
d < wire.size(); ++
d, ++i)
333 if ((wire >= marginW) && (wire <
fNWires - marginW) &&
334 (scaledDrift >= marginD) && (scaledDrift <
fNScaledDrifts - marginD))
return true;
345 fWireProducerLabel(config.WireLabel()),
346 fHitProducerLabel(config.HitLabel()),
347 fTrackModuleLabel(config.TrackLabel()),
348 fSimulationProducerLabel(config.SimulationLabel()),
349 fSaveVtxFlags(config.SaveVtxFlags()),
350 fAdcDelay(config.AdcDelayTicks()),
351 fEventsPerBin(100, 0)
383 std::vector<float>
const & edeps, std::vector<int>
const & pdgs,
size_t wireIdx)
385 if ((wireIdx >=
fWireDriftEdep.size()) || (edeps.size() != pdgs.size())) {
return false; }
397 size_t i0 = i * dstep;
398 size_t i1 = (i + 1) * dstep;
403 float max_edep = edeps[i0];
404 for (
size_t k = i0 + 1; k < i1; ++k)
418 best_pdg |= type_flags;
429 auto const* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
431 wd.
Wire = 0; wd.Drift = 0; wd.TPC = -1; wd.Cryo = -1;
435 double vtx[3] = {tvec.X(), tvec.Y(), tvec.Z()};
439 unsigned int tpc = tpcid.
TPC, cryo = tpcid.
Cryostat;
442 float dx = tvec.T() * 1.e-3 * detprop->DriftVelocity();
444 if (driftDir == 1) { dx *= -1; }
445 else if (driftDir != -1)
447 throw cet::exception(
"nnet::TrainingDataAlg") <<
"drift direction is not X." << std::endl;
449 vtx[0] = tvec.X() + dx;
453 wd.TPC = tpc; wd.Cryo = cryo;
458 mf::LogWarning(
"TrainingDataAlg") <<
"Vertex projection out of wire planes, just skipping this vertex.";
462 mf::LogWarning(
"TrainingDataAlg") <<
"Vertex projection out of wire planes, skip MC vertex.";
469 const std::unordered_map< int, const simb::MCParticle* > & particleMap)
const 471 const float minElectronLength2 = 2.5*2.5;
472 const float maxDeltaLength2 = 0.15*0.15;
474 int pdg = abs(particle.
PdgCode());
475 if (pdg != 11)
return false;
478 for (
size_t d = 0;
d < nSec; ++
d)
480 auto d_search = particleMap.find(particle.
Daughter(
d));
481 if (d_search != particleMap.end())
483 auto const & daughter = *((*d_search).second);
484 int d_pdg = abs(daughter.PdgCode());
485 if (d_pdg != 22) {
return false; }
489 float trkLength2 = 0;
490 auto const * p = &particle;
491 bool branching =
false;
495 auto m_search = particleMap.find(p->Mother());
496 if (m_search != particleMap.end())
498 p = (*m_search).second;
499 int m_pdg = abs(p->PdgCode());
502 nSec = p->NumberDaughters();
504 for (
size_t d = 0;
d < nSec; ++
d)
506 auto d_search = particleMap.find(p->Daughter(
d));
507 if (d_search != particleMap.end())
509 auto const & daughter = *((*d_search).second);
510 int d_pdg = abs(daughter.PdgCode());
517 if (ne > 1) { branching =
true; }
524 return (trkLength2 > minElectronLength2);
528 const std::unordered_map< int, const simb::MCParticle* > & particleMap)
const 530 bool hasElectron =
false, hasNuMu =
false, hasNuE =
false;
532 int pdg = abs(particle.
PdgCode());
533 if ((pdg == 13) && (particle.
EndProcess() ==
"FastScintillation"))
536 for (
size_t d = 0;
d < nSec; ++
d)
538 auto d_search = particleMap.find(particle.
Daughter(
d));
539 if (d_search != particleMap.end())
541 auto const & daughter = *((*d_search).second);
542 int d_pdg = abs(daughter.PdgCode());
543 if (d_pdg == 11) hasElectron =
true;
544 else if (d_pdg == 14) hasNuMu =
true;
545 else if (d_pdg == 12) hasNuE =
true;
550 return (hasElectron && hasNuMu && hasNuE);
554 std::unordered_map<
size_t, std::unordered_map< int, int > > & wireToDriftToVtxFlags,
555 const std::unordered_map< int, const simb::MCParticle* > & particleMap,
556 unsigned int plane)
const 558 for (
auto const & p : particleMap)
560 auto const & particle = *p.second;
562 double ekStart = 1000. * (particle.E() - particle.Mass());
563 double ekEnd = 1000. * (particle.EndE() - particle.Mass());
565 int pdg = abs(particle.PdgCode());
572 if ((particle.EndProcess() ==
"conv") &&
606 if (particle.Mother() != 0)
608 auto search = particleMap.find(particle.Mother());
609 if (search != particleMap.end())
611 auto const & mother = *((*search).second);
612 int m_pdg = abs(mother.PdgCode());
613 unsigned int nSec = mother.NumberDaughters();
614 unsigned int nVisible = 0;
617 for (
size_t d = 0;
d < nSec; ++
d)
619 auto d_search = particleMap.find(mother.Daughter(
d));
620 if (d_search != particleMap.end())
622 auto const & daughter = *((*d_search).second);
623 int d_pdg = abs(daughter.PdgCode());
624 if (((d_pdg == 2212) || (d_pdg == 211) || (d_pdg == 321)) &&
625 (1000. * (daughter.E() - daughter.Mass()) > 50.0))
635 if (((m_pdg != pdg) && (m_pdg != 2112)) || ((m_pdg != 2112) && (nVisible > 0)) || ((m_pdg == 2112) && (nVisible > 1)))
646 if (particle.EndProcess() ==
"FastScintillation")
648 unsigned int nSec = particle.NumberDaughters();
649 for (
size_t d = 0;
d < nSec; ++
d)
651 auto d_search = particleMap.find(particle.Daughter(
d));
652 if (d_search != particleMap.end())
654 auto const & daughter = *((*d_search).second);
655 int d_pdg = abs(daughter.PdgCode());
656 if ((pdg == 321) && (d_pdg == 13))
662 if ((pdg == 211) && (d_pdg == 13))
672 if ((particle.EndProcess() ==
"Decay") && (ekEnd > 200.0))
674 unsigned int nSec = particle.NumberDaughters();
675 for (
size_t d = 0;
d < nSec; ++
d)
677 auto d_search = particleMap.find(particle.Daughter(
d));
678 if (d_search != particleMap.end())
680 auto const & daughter = *((*d_search).second);
681 int d_pdg = abs(daughter.PdgCode());
682 if ((pdg == 321) && (d_pdg == 13))
688 if ((pdg == 211) && (d_pdg == 13))
703 if (particle.Process() ==
"primary")
713 if ((wd.TPC == (
int)
fTPC) && (wd.Cryo == (int)
fCryo))
715 wireToDriftToVtxFlags[wd.Wire][wd.Drift] |= flagsStart;
723 if ((wd.TPC == (
int)
fTPC) && (wd.Cryo == (int)
fCryo))
726 wireToDriftToVtxFlags[wd.Wire][wd.Drift] |= flagsEnd;
744 unsigned int plane,
unsigned int tpc,
unsigned int cryo)
748 std::vector< art::Ptr<recob::Wire> > Wirelist;
754 mf::LogError(
"TrainingDataAlg") <<
"Wire data not set.";
760 std::vector< art::Ptr<recob::Hit> > Hitlist;
767 std::vector< art::Ptr<recob::Track> > Tracklist;
775 for (
size_t widx = 0; widx < 240; ++widx) {
777 std::vector< float > labels_deposit(
fNDrifts, 0);
778 std::vector< int > labels_pdg(
fNDrifts, 0);
781 for(
size_t subwidx = 0; subwidx < Wirelist.size(); ++subwidx) {
782 if(widx+240 == Wirelist[subwidx]->Channel()) {
783 labels_deposit = Wirelist[subwidx]->Signal();
799 for(
size_t iHit = 0; iHit < Hitlist.size(); ++iHit) {
801 if(Hitlist[iHit]->Channel() != widx+240) {
continue; }
802 if(Hitlist[iHit]->View() != 1) {
continue; }
805 if(ass_trk_hits.at(iHit).size() == 0) {
continue; }
810 if(ass_trk_hits.at(iHit)[0]->Length() < 5) {
continue; }
816 for(
size_t jHit = 0; jHit < Hitlist.size(); ++jHit) {
817 if(jHit == iHit) {
continue; }
818 if(Hitlist[jHit]->View() != 1) {
continue; }
820 if(ass_trk_hits.at(jHit).size() == 0) {
continue; }
821 if(ass_trk_hits.at(jHit)[0]->ID() !=
822 ass_trk_hits.at(iHit)[0]->ID()) {
continue; }
824 double dist = sqrt((Hitlist[iHit]->Channel()-Hitlist[jHit]->Channel()) *
825 (Hitlist[iHit]->Channel()-Hitlist[jHit]->Channel()) +
826 (Hitlist[iHit]->PeakTime()-Hitlist[jHit]->PeakTime()) *
827 (Hitlist[iHit]->PeakTime()-Hitlist[jHit]->PeakTime()));
839 for(
size_t jHit = 0; jHit < Hitlist.size(); ++jHit) {
840 if(jHit == iHit or
int(jHit) == far_index) {
continue; }
841 if(Hitlist[jHit]->View() != 1) {
continue; }
843 if(ass_trk_hits.at(jHit).size() == 0) {
continue; }
844 if(ass_trk_hits.at(jHit)[0]->ID() !=
845 ass_trk_hits.at(iHit)[0]->ID()) {
continue; }
847 double dist = sqrt((Hitlist[far_index]->Channel()-Hitlist[jHit]->Channel()) *
848 (Hitlist[far_index]->Channel()-Hitlist[jHit]->Channel()) +
849 (Hitlist[far_index]->PeakTime()-Hitlist[jHit]->PeakTime()) *
850 (Hitlist[far_index]->PeakTime()-Hitlist[jHit]->PeakTime()));
852 if(other_dist < dist){
859 double del_wire = double(Hitlist[other_end]->Channel() - Hitlist[far_index]->Channel());
860 double del_time = double(Hitlist[other_end]->PeakTime() - Hitlist[far_index]->PeakTime());
861 double hypo = sqrt(del_wire * del_wire + del_time * del_time);
863 if(hypo == 0) {
continue; }
865 double cosser = TMath::Abs(del_wire / hypo);
866 double norm_ang = TMath::ACos(cosser) * 2 / TMath::Pi();
878 labels_pdg[Hitlist[iHit]->PeakTime()] = 211;
897 unsigned int plane,
unsigned int tpc,
unsigned int cryo)
904 mf::LogError(
"TrainingDataAlg") <<
"Wire data not set.";
910 mf::LogInfo(
"TrainingDataAlg") <<
"Skip MC simulation info.";
921 std::unordered_map< int, const simb::MCParticle* > particleMap;
922 for (
auto const & particle : *particleHandle)
924 particleMap[particle.TrackId()] = &particle;
927 std::unordered_map< size_t, std::unordered_map< int, int > > wireToDriftToVtxFlags;
932 std::map< int, int > trackToPDG;
933 for (
size_t widx = 0; widx <
fNWires; ++widx)
938 std::vector< float > labels_deposit(
fNDrifts, 0);
939 std::vector< int > labels_pdg(labels_deposit.size(), 0);
940 int labels_size = labels_deposit.size();
942 std::map< int, std::map< int, double > > timeToTrackToCharge;
943 for (
auto const & channel : *simChannelHandle)
945 if (channel.Channel() != wireChannelNumber)
continue;
947 auto const & timeSlices = channel.TDCIDEMap();
948 for (
auto const & timeSlice : timeSlices)
950 int time = timeSlice.first;
952 auto const & energyDeposits = timeSlice.second;
953 for (
auto const & energyDeposit : energyDeposits)
956 int tid = energyDeposit.trackID;
959 pdg = 11; tid = -tid;
961 auto search = particleMap.find(tid);
962 if (search == particleMap.end())
967 auto const & mother = *((*search).second);
968 int mPdg = abs(mother.PdgCode());
969 if ((mPdg == 13) || (mPdg == 211) || (mPdg == 2212))
976 auto search = particleMap.find(tid);
977 if (search == particleMap.end())
982 auto const & particle = *((*search).second);
983 pdg = abs(particle.PdgCode());
985 if (particle.Process() ==
"primary")
997 auto msearch = particleMap.find(particle.Mother());
998 if (msearch != particleMap.end())
1000 auto const & mother = *((*msearch).second);
1011 trackToPDG[energyDeposit.trackID] = pdg;
1013 double energy = energyDeposit.numElectrons * electronsToGeV;
1014 timeToTrackToCharge[time][energyDeposit.trackID] +=
energy;
1022 for (
auto const & ttc : timeToTrackToCharge)
1024 float max_deposit = 0.0;
1026 for (
auto const & tc : ttc.second) {
1028 if( tc.second > max_deposit )
1030 max_deposit = tc.second;
1031 max_pdg = trackToPDG[tc.first];
1035 if (ttc.first < labels_size)
1038 if (tick_idx < labels_size)
1040 labels_deposit[tick_idx] = max_deposit;
1041 labels_pdg[tick_idx] = max_pdg & type_pdg_mask;
1046 for (
auto const & drift_flags : wireToDriftToVtxFlags[widx])
1048 int drift = drift_flags.first, flags = drift_flags.second;
1049 if ((drift >= 0) && (drift < labels_size))
1051 labels_pdg[drift] |= flags;
1067 float max_cut = 0.25 * max_e_cut;
1074 if (cut < max_cut) w0++;
1082 if (cut < max_cut) w1--;
1092 if (cut < max_cut) d0++;
1100 if (cut < max_cut) d1--;
1105 unsigned int margin = 20;
1106 if ((w1 - w0 > 8) && (d1 - d0 > 8))
1108 if (w0 < margin) w0 = 0;
1114 if (d0 < margin) d0 = 0;
bool isInsideFiducialRegion(unsigned int wire, float drift) const
virtual void resizeView(size_t wires, size_t drifts)
Store parameters for running LArG4.
geo::GeometryCore const * fGeometry
bool bufferPatch(size_t wire, float drift, std::vector< std::vector< float > > &patch) const
void collectVtxFlags(std::unordered_map< size_t, std::unordered_map< int, int > > &wireToDriftToVtxFlags, const std::unordered_map< int, const simb::MCParticle * > &particleMap, unsigned int plane) const
bool setWireEdepsAndLabels(std::vector< float > const &edeps, std::vector< int > const &pdgs, size_t wireIdx)
std::vector< raw::ChannelID_t > fWireChannels
Utilities related to art service access.
bool isMuonDecaying(const simb::MCParticle &particle, const std::unordered_map< int, const simb::MCParticle * > &particleMap) const
std::vector< std::vector< float > > Run(std::vector< std::vector< std::vector< float > > > const &inps, int samples=-1) override
std::vector< std::vector< int > > fWireDriftPdg
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
art::InputTag fTrackModuleLabel
bool isValid
Whether this ID points to a valid element.
static std::vector< float > flattenData2D(std::vector< std::vector< float > > const &patch)
virtual std::vector< float > Run(std::vector< std::vector< float > > const &inp2d)=0
virtual void set_data(std::vector< std::vector< std::vector< float > > > const &)
CryostatID_t Cryostat
Index of cryostat.
DataProviderAlg(const fhicl::ParameterSet &pset)
bool setEventData(const art::Event &event, unsigned int plane, unsigned int tpc, unsigned int cryo)
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
std::vector< float > compute_output(keras::DataChunk *dc)
bool isCurrentPatch(unsigned int wire, float drift) const
int NumberDaughters() const
~TrainingDataAlg(void) override
std::vector< std::string > fNNetOutputs
bool setWireDriftData(const std::vector< recob::Wire > &wires, unsigned int plane, unsigned int tpc, unsigned int cryo)
static std::unique_ptr< Graph > create(const char *graph_file_name, const std::vector< std::string > &outputs={})
int Daughter(const int i) const
virtual double ConvertXToTicks(double X, int p, int t, int c) const =0
art::InputTag fWireProducerLabel
geo::TPCID FindTPCAtPosition(double const worldLoc[3]) const
Returns the ID of the TPC at specified location.
~PointIdAlg(void) override
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
bool isElectronEnd(const simb::MCParticle &particle, const std::unordered_map< int, const simb::MCParticle * > &particleMap) const
std::vector< float > predictIdVector(unsigned int wire, float drift) const
calculate multi-class probabilities for [wire, drift] point
constexpr ChannelID_t InvalidChannelID
ID of an invalid channel.
std::vector< std::vector< float > > fWireDriftPatch
std::string EndProcess() const
detinfo::DetectorProperties const * fDetProp
fhicl::Sequence< std::string > NNetOutputs
unsigned int fNCachedDrifts
geo::WireID::WireID_t NearestWire(geo::Point_t const &point, geo::PlaneID const &planeid) const
Returns the index of wire closest to position in the specified TPC.
std::vector< size_t > fEventsPerBin
KerasModelInterface(const char *modelFileName)
void fill(const art::PtrVector< recob::Hit > &hits, int only_plane)
bool findCrop(float max_e_cut, unsigned int &w0, unsigned int &w1, unsigned int &d0, unsigned int &d1) const
WireDrift getProjection(const TLorentzVector &tvec, unsigned int plane) const
The data type to uniquely identify a TPC.
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
fhicl::Atom< std::string > NNetModelFile
TfModelInterface(const char *modelFileName)
std::vector< float > Run(std::vector< std::vector< float > > const &inp2d) override
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
std::string findFile(const char *fileName) const
short int DetectDriftDirection() const
Returns the expected drift direction based on geometry.
std::vector< std::vector< float > > predictIdVectors(std::vector< std::pair< unsigned int, float > > points) const
unsigned int fNScaledDrifts
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
static float particleRange2(const simb::MCParticle &particle)
object containing MC truth information necessary for making RawDigits and doing back tracking ...
size_t fCurrentScaledDrift
nnet::ModelInterface * fNNet
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
Exception thrown on invalid wire number (e.g. NearestWireID())
art::InputTag fSimulationProducerLabel
Interface to algorithm class for a specific detector channel mapping.
std::string fNNetModelFilePath
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
TPCID_t TPC
Index of the TPC within its cryostat.
bool isSamePatch(unsigned int wire1, float drift1, unsigned int wire2, float drift2) const
test if two wire/drift coordinates point to the same patch
bool setDataEventData(const art::Event &event, unsigned int plane, unsigned int tpc, unsigned int cryo)
std::vector< std::vector< float > > fWireDriftEdep
void resizeView(size_t wires, size_t drifts) override
art::InputTag fHitProducerLabel
double GeVToElectrons() const
cet::coded_exception< error, detail::translate > exception
Event finding and building.
TrainingDataAlg(const fhicl::ParameterSet &pset)
float predictIdValue(unsigned int wire, float drift, size_t outIdx=0) const
calculate single-value prediction (2-class probability) for [wire, drift] point
PointIdAlg(const fhicl::ParameterSet &pset)