35 #include "art_root_io/TFileService.h" 39 #include "cetlib_except/exception.h" 53 #include <unordered_map> 69 Name(
"CalorimetryAlg"),
70 Comment(
"Used to calculate electron lifetime correction.")};
73 Comment(
"Simulation producer")};
76 Name(
"PfpModuleLabel"),
77 Comment(
"PFP producer tag, to compare with NNet results")};
103 void countTruthDep(
const std::vector<sim::SimChannel>& channels,
105 float& trackLike)
const;
109 const std::vector<recob::PFParticle>& pfparticles,
113 float& trackLike)
const;
116 const std::unordered_map<int, const simb::MCParticle*>& particleMap)
const;
120 const std::vector<sim::SimChannel>& channels,
122 const std::array<float, MVA_LENGTH>& cnn_out,
197 fEventTree = tfs->make<TTree>(
"event",
"event info");
225 fClusterTree = tfs->make<TTree>(
"cluster",
"clusters info");
231 fHitTree = tfs->make<TTree>(
"hits",
"hits info");
247 for (
size_t i = 0; i < 100; ++i) {
261 TTree* thrTree = tfs->make<TTree>(
"threshold",
"error rate vs threshold");
263 float thr, shErr, trkErr;
264 thrTree->Branch(
"thr", &thr,
"thr/F");
265 thrTree->Branch(
"shErr", &shErr,
"shErr/F");
266 thrTree->Branch(
"trkErr", &trkErr,
"trkErr/F");
268 for (
size_t i = 0; i < 100; ++i) {
294 for (
size_t i = 0; i < 100; ++i) {
322 for (
auto const& particle : *particleHandle) {
354 <<
"No em/track labeled columns in MVA data products." << std::endl;
361 cluResults->dataHandle(),
e, cluResults->dataTag());
363 for (
size_t c = 0; c < cluResults->size(); ++c) {
368 const std::vector<art::Ptr<recob::Hit>>&
hits = hitsFromClusters.at(c);
369 std::array<float, MVA_LENGTH> cnn_out = cluResults->getOutput(c);
397 for (
size_t i = 0; i < 100; ++i) {
403 if (fHitTrack_0p5 > 0)
427 mf::LogWarning(
"TrainingDataAlg") <<
"MVA FOR CLUSTERS NOT FOUND";
436 float& trackLike)
const 440 for (
auto const& channel : channels) {
442 auto const& timeSlices = channel.TDCIDEMap();
443 for (
auto const& timeSlice : timeSlices) {
445 auto const& energyDeposits = timeSlice.second;
446 for (
auto const& energyDeposit : energyDeposits) {
447 int trackID = energyDeposit.trackID;
451 if (trackID < 0) { emLike +=
energy; }
452 else if (trackID > 0) {
463 if (!pdg) pdg = particle.
PdgCode();
466 if ((pdg == 11) || (pdg == -11) || (pdg == 22))
479 const std::vector<recob::PFParticle>& pfparticles,
483 float& trackLike)
const 487 for (
size_t i = 0; i < pfparticles.size(); ++i) {
488 const auto& pfp = pfparticles[i];
489 const auto& clus = pfpclus.at(i);
492 for (
const auto& c : clus) {
493 const auto&
hits = cluhits.at(c.key());
494 for (
const auto& h :
hits) {
495 if (h->View() ==
fView) {
502 if ((pfp.PdgCode() == 11) || pfp.PdgCode() == -11) { emLike += hitdep; }
512 const std::unordered_map<int, const simb::MCParticle*>& particleMap)
const 514 bool hasElectron =
false, hasNuMu =
false, hasNuE =
false;
517 if ((pdg == 13) && (particle.
EndProcess() ==
"FastScintillation"))
520 for (
size_t d = 0;
d < nSec; ++
d) {
521 auto d_search = particleMap.find(particle.
Daughter(
d));
522 if (d_search != particleMap.end()) {
523 auto const& daughter = *((*d_search).second);
524 int d_pdg =
abs(daughter.PdgCode());
527 else if (d_pdg == 14)
529 else if (d_pdg == 12)
535 return (hasElectron && hasNuMu && hasNuE);
541 const std::vector<sim::SimChannel>& channels,
543 const std::array<float, MVA_LENGTH>& cnn_out,
549 std::unordered_map<int, int> mcHitPid;
558 double totEnSh = 0, totEnTrk = 0;
561 auto hitChannelNumber =
hit->Channel();
563 double hitEn = 0, hitEnSh = 0, hitEnTrk = 0, hitEnMichel = 0;
565 auto const& vout = hit_outs[
hit.key()];
571 for (
auto const& channel : channels) {
572 if (channel.Channel() != hitChannelNumber)
continue;
575 auto const& timeSlices = channel.TDCIDEMap();
576 for (
auto const& timeSlice : timeSlices) {
577 int time = timeSlice.first;
578 if (
std::abs(
hit->TimeDistanceAsRMS(time)) < 1.0) {
580 auto const& energyDeposits = timeSlice.second;
582 for (
auto const& energyDeposit : energyDeposits) {
583 int trackID = energyDeposit.trackID;
588 if (trackID < 0) { hitEnSh +=
energy; }
589 else if (trackID > 0) {
595 if ((pdg == 11) || (pdg == -11) || (pdg == 22))
604 auto const& mother = *((*msearch).second);
618 totEnTrk += hitEnTrk;
624 int hitPidMc_0p5 = -1;
625 if (hitEnSh > hitEnTrk) {
633 mcHitPid[
hit.key()] = hitPidMc_0p5;
634 auto const& hout = hit_outs[
hit.key()];
638 fHitEM_mc = hitEnSh / (hitEnSh + hitEnTrk);
640 int hitPidMc_0p85 = -1;
641 double hitDep = hitEnSh + hitEnTrk;
642 if (hitEnSh > 0.85 * hitDep) {
647 else if (hitEnTrk > 0.85 * hitDep) {
653 bool mcMichel =
false;
663 for (
size_t i = 0; i < 100; ++i) {
664 double thr = 0.01 * i;
666 if (p_michel > thr) {
704 if (totEnSh > 1.5 * totEnTrk)
708 else if (totEnTrk > 1.5 * totEnSh) {
712 for (
size_t i = 0; i < 100; ++i) {
713 double thr = 0.01 * i;
743 for (
auto const& h : hits) {
744 auto const& vout = hit_outs[h.key()];
745 double hitPidValue = 0;
747 if (h_trk_or_sh > 0) hitPidValue = vout[
fTrkLikeIdx] / h_trk_or_sh;
750 <<
" " << h->PeakTime() <<
" " 753 <<
" " << mcHitPid[h.key()] <<
" " <<
fPidValue <<
" " << hitPidValue;
float fHitsTrack_OK_0p85[100]
art::InputTag fSimulationProducerLabel
Store parameters for running LArG4.
std::ofstream fHitsOutFile
float fHitsMichel_OK_0p5[100]
Declaration of signal hit object.
art::InputTag fNNetModuleLabel
constexpr auto abs(T v)
Returns the absolute value of the argument.
int getIndex(const std::string &name) const
Index of column with given name, or -1 if name not found.
void countPfpDep(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const std::vector< recob::PFParticle > &pfparticles, const art::FindManyP< recob::Cluster > &pfpclus, const art::FindManyP< recob::Hit > &cluhits, float &emLike, float &trackLike) const
Set of hits with a 2D structure.
void countTruthDep(const std::vector< sim::SimChannel > &channels, float &emLike, float &trackLike) const
geo::PlaneID Plane() const
Returns the plane ID this cluster lies on.
PointIdEffTest & operator=(PointIdEffTest const &)=delete
EDAnalyzer(fhicl::ParameterSet const &pset)
fhicl::Atom< bool > SaveHitsFile
int NumberDaughters() const
int Daughter(const int i) const
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
void analyze(art::Event const &e) override
std::unordered_map< int, const simb::MCParticle * > fParticleMap
PointIdEffTest(Parameters const &config)
#define DEFINE_ART_MODULE(klass)
void beginRun(const art::Run &run) override
std::string EndProcess() const
int testCNN(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const std::vector< sim::SimChannel > &channels, const std::vector< art::Ptr< recob::Hit >> &hits, const std::array< float, MVA_LENGTH > &cnn_out, const std::vector< anab::FeatureVector< MVA_LENGTH >> &hit_outs, size_t cidx)
bool isMuonDecaying(const simb::MCParticle &particle, const std::unordered_map< int, const simb::MCParticle * > &particleMap) const
float fHitsEM_OK_0p5[100]
std::vector< FeatureVector< N > > const & outputs() const
Access the vector of the feature vectors.
Provides recob::Track data product.
float fHitsTrack_OK_0p5[100]
fhicl::Atom< art::InputTag > PfpModuleLabel
PlaneID_t Plane
Index of the plane within its TPC.
fhicl::Table< calo::CalorimetryAlg::Config > CalorimetryAlg
Declaration of cluster object.
float fHitRecoFractionEM[100]
Detector simulation of raw signals on wires.
fhicl::Atom< unsigned int > View
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
fhicl::Atom< art::InputTag > SimModuleLabel
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Contains all timing reference information for the detector.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
object containing MC truth information necessary for making RawDigits and doing back tracking ...
fhicl::Atom< art::InputTag > NNetModuleLabel
EventNumber_t event() const
Declaration of basic channel signal object.
double LifetimeCorrection(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, double time, double T0=0) const
Collection of Physical constants used in LArSoft.
float fHitsMichel_False_0p5[100]
calo::CalorimetryAlg fCalorimetryAlg
static std::unique_ptr< MVAReader > create(const art::Event &evt, const art::InputTag &tag)
art::InputTag fPfpModuleLabel
double GeVToElectrons() const
cet::coded_exception< error, detail::translate > exception
float fHitsEM_OK_0p85[100]