54 const std::vector<const sim::Particle*>& particles);
74 std::map<unsigned int, MRCCParticle*>
Particles()
const;
99 float DepositedCharge()
const;
100 float DepositedMRCCCharge()
const;
110 std::vector<art::Ptr<rb::CellHit> >
fHits;
117 const std::vector<const sim::Particle*>& particles) {
120 for (std::vector<const sim::Particle*>::const_iterator particleIt = particles.begin();
121 particleIt != particles.end(); ++particleIt) {
124 fMuon =
new unsigned int((*particleIt)->TrackId());
175 for (std::map<unsigned int, MRCCParticle*>::const_iterator particleIt =
fParticles.begin();
177 charge += particleIt->second->DepositedCharge();
184 for (std::map<unsigned int, MRCCParticle*>::const_iterator particleIt =
fParticles.begin();
186 charge += particleIt->second->DepositedMRCCCharge();
193 for (std::map<unsigned int, MRCCParticle*>::const_iterator particleIt =
fParticles.begin();
194 particleIt !=
fParticles.end(); ++particleIt) {
197 charge += particleIt->second->DepositedCharge();
205 for (std::map<unsigned int, MRCCParticle*>::const_iterator particleIt =
fParticles.begin();
206 particleIt !=
fParticles.end(); ++particleIt) {
209 charge += particleIt->second->DepositedMRCCCharge();
218 charge =
fParticles.at(*fMuon)->DepositedCharge();
226 charge =
fParticles.at(*fMuon)->DepositedMRCCCharge();
249 fTotalMRCCCharge(0.),
250 fTrueParticle(particle) {
255 fHits.push_back(hit);
349 fTruthModuleLabel(pset.
get<
std::
string>("TruthModuleLabel")),
350 fHitsModuleLabel(pset.
get<
std::
string>("HitsModuleLabel")),
351 fClusterModuleLabel(pset.
get<
std::
string>("ClusterModuleLabel")),
352 fElasticModuleLabel(pset.
get<
std::
string>("ElasticModuleLabel")),
353 fFuzzy3DModuleLabel(pset.
get<
std::
string>("Fuzzy3DModuleLabel")),
354 fBPFTrackModuleLabel(pset.
get<
std::
string>("BPFTrackModuleLabel")),
355 fMRCCHitsModuleLabel(pset.
get<
std::
string>("MRCCHitsModuleLabel")),
356 fMRCCClusterModuleLabel(pset.
get<
std::
string>("MRCCClusterModuleLabel")) {
367 = tfs->make<TH1I>(
"NumNuClusters",
";Number of Clusters Associated with True Neutrino;", 5,0,5);
368 hNumNuClusters->GetXaxis()->CenterLabels();
370 = tfs->make<TH1I>(
"NumNuClustersMRCC",
";Number of Clusters Associated with True Neutrino;", 5,0,5);
371 hNumNuClustersMRCC->GetXaxis()->CenterLabels();
373 = tfs->make<TProfile>(
"SliceHadronicFrac",
";Neutrino Energy (GeV);Hadronic Fraction of Total Slice Charge;", 100,0,40);
374 hSliceHadronicFrac->SetMinimum(0);
375 hSliceHadronicFrac->SetMaximum(1);
376 hSliceHadronicFracMRCC
377 = tfs->make<TProfile>(
"SliceHadronicFracMRCC",
";Neutrino Energy (GeV);Hadronic Fraction of Total Slice Charge;", 100,0,40);
378 hSliceHadronicFracMRCC->SetMinimum(0);
379 hSliceHadronicFracMRCC->SetMaximum(1);
380 hAnyMRCCNeutrinoEnergy
381 = tfs->make<TProfile>(
"AnyMRCCNeutrinoEnergy",
";Neutrino Energy (GeV);Fraction of Neutrino Interactions with Charge Removed;", 100,0,40);
382 hAnyMRCCNeutrinoEnergy->SetMinimum(0);
383 hAnyMRCCNeutrinoEnergy->SetMaximum(1);
385 = tfs->make<TProfile>(
"AnyMRCCMuonEnergy",
";Muon Energy (GeV);Fraction of CC Muons with Charge Removed;", 100,0,40);
386 hAnyMRCCMuonEnergy->SetMinimum(0);
387 hAnyMRCCMuonEnergy->SetMaximum(1);
389 = tfs->make<TProfile>(
"SliceEFracRemoved",
";Neutrino Energy (GeV);Fraction of Slice Charge Removed;", 100,0,40);
390 hSliceEFracRemoved->SetMinimum(0);
391 hSliceEFracRemoved->SetMaximum(1);
392 hSliceHadEFracRemoved
393 = tfs->make<TProfile>(
"SliceHadEFracRemoved",
";Neutrino Energy (GeV);Fraction of Slice Hadronic Charge Removed;", 100,0,40);
394 hSliceHadEFracRemoved->SetMinimum(0);
395 hSliceHadEFracRemoved->SetMaximum(1);
396 hSliceLepEFracRemoved
397 = tfs->make<TProfile>(
"SliceLepEFracRemoved",
";Neutrino Energy (GeV);Fraction of Slice Leptonic Charge Removed;", 100,0,40);
398 hSliceLepEFracRemoved->SetMinimum(0);
399 hSliceLepEFracRemoved->SetMaximum(1);
401 = tfs->make<TH1F>(
"SliceDeltaE",
";Fraction of Slice Charge Removed (GeV);", 100,0,1);
403 = tfs->make<TH1F>(
"SliceDeltaLepE",
";Fraction of Slice Leptonic Charge Removed (GeV);", 100,0,1);
405 = tfs->make<TH1F>(
"SliceDeltaHadE",
";Fraction of Slice Hadronic Charge Removed (GeV);", 100,0,1);
407 = tfs->make<TProfile>(
"LeptonEFracRemoved",
";Lepton Energy (GeV);Fraction of Charge Removed;", 100,0,40);
408 hLeptonEFracRemoved->SetMinimum(0);
409 hLeptonEFracRemoved->SetMaximum(1);
411 = tfs->make<TProfile>(
"HadronEFracRemoved",
";Hadron Energy (GeV);Fraction of Charge Removed;", 100,0,40);
412 hHadronEFracRemoved->SetMinimum(0);
413 hHadronEFracRemoved->SetMaximum(1);
414 hMRCCSliceEFracRemoved
415 = tfs->make<TProfile>(
"MRCCSliceEFracRemoved",
";Neutrino Energy (GeV);Fraction of Slice Charge Removed;", 100,0,40);
416 hMRCCSliceEFracRemoved->SetMinimum(0);
417 hMRCCSliceEFracRemoved->SetMaximum(1);
418 hMRCCSliceHadEFracRemoved
419 = tfs->make<TProfile>(
"MRCCSliceHadEFracRemoved",
";Neutrino Energy (GeV);Fraction of Slice Hadronic Charge Removed;", 100,0,40);
420 hMRCCSliceHadEFracRemoved->SetMinimum(0);
421 hMRCCSliceHadEFracRemoved->SetMaximum(1);
422 hMRCCSliceLepEFracRemoved
423 = tfs->make<TProfile>(
"MRCCSliceLepEFracRemoved",
";Neutrino Energy (GeV);Fraction of Slice Leptonic Charge Removed;", 100,0,40);
424 hMRCCSliceLepEFracRemoved->SetMinimum(0);
425 hMRCCSliceLepEFracRemoved->SetMaximum(1);
426 hMRCCLeptonEFracRemoved
427 = tfs->make<TProfile>(
"MRCCLeptonEFracRemoved",
";Lepton Energy (GeV);Fraction of Charge Removed;", 100,0,40);
428 hMRCCLeptonEFracRemoved->SetMinimum(0);
429 hMRCCLeptonEFracRemoved->SetMaximum(1);
430 hMRCCHadronEFracRemoved
431 = tfs->make<TProfile>(
"MRCCHadronEFracRemoved",
";Hadron Energy (GeV);Fraction of Charge Removed;", 100,0,40);
432 hMRCCHadronEFracRemoved->SetMinimum(0);
433 hMRCCHadronEFracRemoved->SetMaximum(1);
434 hHadYMRCC = tfs->make<TH1D>(
"HadYMRCC",
";Hadronic Y;", 100, 0, 1);
435 hHadYNoMRCC = tfs->make<TH1D>(
"HadYNoMRCC",
";Hadronic Y;", 100, 0, 1);
437 = tfs->make<TProfile>(
"SliceEFracRemovedQE",
";Fraction of Slice Charge Removed (QE);", 100,0,40);
438 hSliceEFracRemovedDIS
439 = tfs->make<TProfile>(
"SliceEFracRemovedDIS",
";Fraction of Slice Charge Removed (DIS);", 100,0,40);
440 hSliceEFracRemovedRES
441 = tfs->make<TProfile>(
"SliceEFracRemovedRES",
";Fraction of Slice Charge Removed (RES);", 100,0,40);
452 std::vector<art::Ptr<simb::MCTruth> > truths;
453 if (evt.
getByLabel(fTruthModuleLabel, truthHandle))
458 std::vector<art::Ptr<rb::CellHit> >
hits;
459 if (evt.
getByLabel(fHitsModuleLabel, hitsHandle))
464 std::vector<art::Ptr<rb::Cluster> > clusters;
465 if (evt.
getByLabel(fClusterModuleLabel, clusterHandle))
470 std::vector<art::Ptr<rb::CellHit> > mrccHits;
471 if (evt.
getByLabel(fMRCCHitsModuleLabel, mrccHitsHandle))
476 std::vector<art::Ptr<rb::Cluster> > mrccClusters;
477 if (evt.
getByLabel(fMRCCClusterModuleLabel, mrccClusterHandle))
481 std::map<unsigned int, MRCCNeutrino*> true_neutrinos;
484 for (std::vector<
art::Ptr<simb::MCTruth> >::const_iterator truthIt = truths.begin(); truthIt != truths.end(); ++truthIt) {
485 std::vector<const sim::Particle*> sim_particles = bt->MCTruthToParticles(*truthIt);
486 true_neutrinos[truthIt->key()] =
new MRCCNeutrino(*truthIt, sim_particles);
491 clusterIt != clusters.end(); ++clusterIt) {
492 if ((*clusterIt)->IsNoise() or !SliceSelection(evt, *clusterIt))
494 std::vector<cheat::NeutrinoEffPur> nus = bt->SliceToNeutrinoInteractions(*clusterIt, hits);
497 true_neutrinos[nus[0].neutrinoInt.key()]->AddAssociatedCluster(clusterIt->key());
500 hitIt != cluster_hits.
end(); ++hitIt) {
501 const sim::Particle* true_particle = bt->HitToParticle(*(*hitIt));
502 if (true_particle ==
nullptr)
504 true_neutrinos[nus[0].neutrinoInt.key()]->AddParticleHit(true_particle->
TrackId(), *hitIt);
510 clusterIt != mrccClusters.end(); ++clusterIt) {
511 if ((*clusterIt)->IsNoise())
513 std::vector<cheat::NeutrinoEffPur> nus = bt->SliceToNeutrinoInteractions(*clusterIt, mrccHits);
516 true_neutrinos[nus[0].neutrinoInt.key()]->AddAssociatedMRCCCluster(clusterIt->key());
519 hitIt != cluster_hits.
end(); ++hitIt) {
520 const sim::Particle* true_particle = bt->HitToParticle(*(*hitIt));
521 if (true_particle ==
nullptr)
523 true_neutrinos[nus[0].neutrinoInt.key()]->AddParticleMRCCHit(true_particle->
TrackId(), *hitIt);
528 for (std::map<unsigned int, MRCCNeutrino*>::const_iterator neutrinoIt = true_neutrinos.begin();
529 neutrinoIt != true_neutrinos.end(); ++neutrinoIt) {
546 bool mrcc = charge_before != charge_after;
549 hAnyMRCCNeutrinoEnergy->Fill(true_nu->
Energy(),
mrcc);
552 hSliceEFracRemoved->Fill(true_nu->
Energy(), (charge_before-charge_after)/charge_before);
553 hSliceDeltaE->Fill((charge_before-charge_after)/charge_before);
555 hSliceEFracRemovedQE->Fill(true_nu->
Energy(), (charge_before-charge_after)/charge_before);
557 hSliceEFracRemovedRES->Fill(true_nu->
Energy(), (charge_before-charge_after)/charge_before);
559 hSliceEFracRemovedDIS->Fill(true_nu->
Energy(), (charge_before-charge_after)/charge_before);
561 hMRCCSliceEFracRemoved->Fill(true_nu->
Energy(), (charge_before-charge_after)/charge_before);
565 if (had_charge_before) {
566 hSliceHadEFracRemoved->Fill(true_nu->
Energy(), (had_charge_before-had_charge_after)/had_charge_before);
567 hSliceDeltaHadE->Fill((had_charge_before-had_charge_after)/had_charge_before);
569 hMRCCSliceHadEFracRemoved->Fill(true_nu->
Energy(), (had_charge_before-had_charge_after)/had_charge_before);
573 if (lep_charge_before) {
574 hSliceLepEFracRemoved->Fill(true_nu->
Energy(), (lep_charge_before-lep_charge_after)/lep_charge_before);
575 hSliceDeltaLepE->Fill((lep_charge_before-lep_charge_after)/lep_charge_before);
577 hMRCCSliceLepEFracRemoved->Fill(true_nu->
Energy(), (lep_charge_before-lep_charge_after)/lep_charge_before);
586 std::map<unsigned int, MRCCParticle*> particles = true_nu->
Particles();
587 int muon = true_nu->
Muon();
588 for (std::map<unsigned int, MRCCParticle*>::const_iterator particleIt = particles.begin();
589 particleIt != particles.end(); ++particleIt) {
590 float p_charge_before = particleIt->second->DepositedCharge();
591 float p_charge_after = particleIt->second->DepositedMRCCCharge();
592 if (!p_charge_before)
594 if ((
int)particleIt->first == muon) {
595 hLeptonEFracRemoved->Fill(particleIt->second->Energy(), (p_charge_before-p_charge_after)/p_charge_before);
596 hAnyMRCCMuonEnergy->Fill(particleIt->second->Energy(),
mrcc);
598 hMRCCLeptonEFracRemoved->Fill(particleIt->second->Energy(), (p_charge_before-p_charge_after)/p_charge_before);
600 hHadronEFracRemoved->Fill(particleIt->second->Energy(), (p_charge_before-p_charge_after)/p_charge_before);
602 hMRCCHadronEFracRemoved->Fill(particleIt->second->Energy(), (p_charge_before-p_charge_after)/p_charge_before);
623 if (!fmElastic.isValid())
626 std::vector<art::Ptr<rb::Vertex> > elastics = fmElastic.at(0);
627 if (elastics.empty())
632 std::vector<art::Ptr<rb::Prong> > prongs = fmFuzzyProng3D.at(0);
637 double muonEmax = -1;
639 std::vector<art::Ptr<rb::Track> > tracks =
fmt.at(0);
640 for (
unsigned int trackIt = 0; trackIt < tracks.size(); ++trackIt) {
642 if (fmBPFFitSums.isValid()) {
643 std::vector<art::Ptr<rb::FitSum> > bpfFitSums = fmBPFFitSums.at(0);
644 if (
abs(bpfFitSums[0]->PDG()) == 13) {
645 double muonE = bpfFitSums[0]->FourMom().E();
646 if (muonE > muonEmax)
double E(const int i=0) const
float ClusterCharge() const
back track the reconstruction to the simulation
float ClusterMRCCCharge() const
MRCCNeutrino(art::Ptr< simb::MCTruth > neutrino_int, const std::vector< const sim::Particle * > &particles)
const simb::MCNeutrino & GetNeutrino() const
TProfile * hSliceEFracRemovedRES
TProfile * hSliceHadEFracRemoved
art::Ptr< simb::MCTruth > fTruth
rb::Cluster cluster_hits(std::vector< art::Ptr< rb::CellHit > > const &hits, int32_t t_start, int32_t t_end, int16_t adc_cut=0)
std::vector< art::Ptr< rb::CellHit > > fMRCCHits
const simb::MCParticle & Nu() const
std::string fTruthModuleLabel
void AddHit(art::Ptr< rb::CellHit > hit)
float ClusterMRCCLeptonicCharge() const
TProfile * hSliceEFracRemoved
unsigned int NumAssociatedClusters() const
DEFINE_ART_MODULE(TestTMapFile)
TProfile * hLeptonEFracRemoved
const simb::MCParticle * TrueParticle() const
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
std::string fFuzzy3DModuleLabel
TProfile * hMRCCSliceHadEFracRemoved
TProfile * hSliceLepEFracRemoved
int InteractionType() const
TProfile * hAnyMRCCNeutrinoEnergy
void AddAssociatedCluster(unsigned int clusterID)
TH1I * hNumNuClustersMRCC
TProfile * hMRCCSliceEFracRemoved
std::string fBPFTrackModuleLabel
bool SliceSelection(const art::Event &evt, const art::Ptr< rb::Cluster > &slice)
float DepositedMRCCCharge() const
std::string fHitsModuleLabel
float DepositedCharge() const
std::string fElasticModuleLabel
art::Ptr< simb::MCTruth > Truth() const
TProfile * hMRCCSliceLepEFracRemoved
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
float ClusterMRCCHadronicCharge() const
TProfile * hSliceHadronicFracMRCC
TProfile * hHadronEFracRemoved
std::string fMRCCClusterModuleLabel
Vertex location in position and time.
TProfile * hMRCCLeptonEFracRemoved
TProfile * hAnyMRCCMuonEnergy
TProfile * hMRCCHadronEFracRemoved
art::ServiceHandle< cheat::BackTracker > bt
std::vector< art::Ptr< rb::CellHit > > fHits
data_t::const_iterator const_iterator
std::string fClusterModuleLabel
float ClusterHadronicCharge() const
TProfile * hSliceHadronicFrac
void AddAssociatedMRCCCluster(unsigned int clusterID)
std::string fMRCCHitsModuleLabel
void AddMRCCHit(art::Ptr< rb::CellHit > hit)
std::vector< unsigned int > fAssociatedClusters
int16_t ADC(uint32_t i) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
std::map< unsigned int, MRCCParticle * > Particles() const
std::map< unsigned int, MRCCParticle * > fParticles
void AddParticleMRCCHit(unsigned int trackID, art::Ptr< rb::CellHit > hit)
TProfile * hSliceEFracRemovedQE
MRCCParticle(const simb::MCParticle *particle)
std::vector< unsigned int > fAssociatedMRCCClusters
void analyze(const art::Event &evt)
art::ServiceHandle< art::TFileService > tfs
TProfile * hSliceEFracRemovedDIS
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
void AddParticleHit(unsigned int trackID, art::Ptr< rb::CellHit > hit)
unsigned int NumAssociatedMRCCClusters() const
const simb::MCParticle * fTrueParticle
float ClusterLeptonicCharge() const