439 art::Handle<std::vector<simb::MCTruth>> MCtruthHandle;
441 std::vector<art::Ptr<simb::MCTruth>> MCtruthlist;
442 art::fill_ptr_vector(MCtruthlist, MCtruthHandle);
443 art::Ptr<simb::MCTruth> MCtruth;
446 MCtruth = MCtruthlist[0];
447 if (MCtruth->NeutrinoSet()) {
448 simb::MCNeutrino
nu = MCtruth->GetNeutrino();
451 else if (nu.CCNC() == 1)
453 simb::MCParticle neutrino = nu.Nu();
455 const TLorentzVector& nu_momentum = nu.Nu().Momentum(0);
457 const TLorentzVector&
vertex = neutrino.Position(0);
463 double tmp_leadingPionPlusE = 0.0;
464 double tmp_leadingPionMinusE = 0.0;
465 double tmp_leadingProtonE = 0.0;
467 simb::MCParticle* MClepton =
nullptr;
468 simb::MCParticle* MCproton =
nullptr;
469 simb::MCParticle* MCpion_plus =
nullptr;
470 simb::MCParticle* MCpion_minus =
nullptr;
471 simb::MCParticle* MCkaon =
nullptr;
472 simb::MCParticle* MCmichel =
nullptr;
474 art::ServiceHandle<cheat::ParticleInventoryService const> pi_serv;
475 const sim::ParticleList& plist = pi_serv->ParticleList();
476 simb::MCParticle* particle = 0;
479 auto const clockData =
480 art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
482 art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(event);
484 for (sim::ParticleList::const_iterator ipar = plist.begin(); ipar != plist.end(); ++ipar) {
485 particle = ipar->second;
486 if (particle->PdgCode() ==
fLeptonPDGcode && particle->Mother() == 0) {
487 const TLorentzVector& lepton_momentum = particle->Momentum(0);
490 MC_leptonP = sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
491 pow(particle->Momentum().Pz(), 2));
494 if (particle->Mother() == 0) {
496 if (particle->PdgCode() == 2212) {
497 if (particle->Momentum().E() > tmp_leadingProtonE) {
498 tmp_leadingProtonE = particle->Momentum().E();
501 sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
502 pow(particle->Momentum().Pz(), 2));
506 else if (particle->PdgCode() == 211) {
507 if (particle->Momentum().E() > tmp_leadingPionPlusE) {
508 tmp_leadingPionPlusE = particle->Momentum().E();
511 sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
512 pow(particle->Momentum().Pz(), 2));
513 MCpion_plus = particle;
516 else if (particle->PdgCode() == -211) {
517 if (particle->Momentum().E() > tmp_leadingPionMinusE) {
518 tmp_leadingPionMinusE = particle->Momentum().E();
521 sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
522 pow(particle->Momentum().Pz(), 2));
523 MCpion_minus = particle;
533 if (particle->Mother() == 0) {
534 const TLorentzVector& positionStart = particle->Position(0);
537 if (particle->PdgCode() == 321) {
539 MC_kaonP = sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
540 pow(particle->Momentum().Pz(), 2));
545 MC_leptonP = sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
546 pow(particle->Momentum().Pz(), 2));
549 else if (particle->PdgCode() == 2212) {
550 if (particle->Momentum().E() > tmp_leadingProtonE) {
551 tmp_leadingProtonE = particle->Momentum().E();
554 sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
555 pow(particle->Momentum().Pz(), 2));
559 else if (particle->PdgCode() == 211) {
560 if (particle->Momentum().E() > tmp_leadingPionPlusE) {
561 tmp_leadingPionPlusE = particle->Momentum().E();
564 sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
565 pow(particle->Momentum().Pz(), 2));
566 MCpion_plus = particle;
569 else if (particle->PdgCode() == -211) {
570 if (particle->Momentum().E() > tmp_leadingPionMinusE) {
571 tmp_leadingPionMinusE = particle->Momentum().E();
574 sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
575 pow(particle->Momentum().Pz(), 2));
576 MCpion_minus = particle;
579 else if (particle->Process() ==
"Decay" &&
580 particle->PdgCode() == -11) {
582 MC_michelP = sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
583 pow(particle->Momentum().Pz(), 2));
586 else if (TMath::Abs(particle->PdgCode() == 321)) {
588 MC_kaonP = sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
589 pow(particle->Momentum().Pz(), 2));
604 theta_mu *= (180.0 / 3.14159);
670 art::Handle<std::vector<recob::Track>> trackListHandle;
672 std::vector<art::Ptr<recob::Track>> tracklist;
673 art::fill_ptr_vector(tracklist, trackListHandle);
674 int n_recoTrack = tracklist.size();
676 art::FindManyP<recob::Hit> track_hits(trackListHandle, event,
fTrackModuleLabel);
677 if (n_recoTrack == 0) {
678 MF_LOG_DEBUG(
"NeutrinoTrackingEff") <<
"There are no reco tracks... bye";
681 MF_LOG_DEBUG(
"NeutrinoTrackingEff") <<
"Found this many reco tracks " << n_recoTrack;
683 double Efrac_lepton = 0.0;
684 double Ecomplet_lepton = 0.0;
685 double Efrac_proton = 0.0;
686 double Ecomplet_proton = 0.0;
687 double Efrac_pionplus = 0.0;
688 double Ecomplet_pionplus = 0.0;
689 double Efrac_pionminus = 0.0;
690 double Ecomplet_pionminus = 0.0;
691 double Efrac_kaon = 0.0;
692 double Ecomplet_kaon = 0.0;
693 double Efrac_michel = 0.0;
694 double Ecomplet_michel = 0.0;
695 double trackLength_lepton = 0.0;
696 double trackLength_proton = 0.0;
697 double trackLength_pion_plus = 0.0;
698 double trackLength_pion_minus = 0.0;
699 double trackLength_kaon = 0.0;
700 double trackLength_michel = 0.0;
701 const simb::MCParticle* MClepton_reco =
nullptr;
702 const simb::MCParticle* MCproton_reco =
nullptr;
703 const simb::MCParticle* MCpion_plus_reco =
nullptr;
704 const simb::MCParticle* MCpion_minus_reco =
nullptr;
705 const simb::MCParticle* MCkaon_reco =
nullptr;
706 const simb::MCParticle* MCmichel_reco =
nullptr;
708 std::vector<art::Ptr<recob::Hit>> tmp_all_trackHits = track_hits.at(0);
709 std::vector<art::Ptr<recob::Hit>> all_hits;
710 art::Handle<std::vector<recob::Hit>> hithandle;
711 auto const pd =
event.getProductDescription(tmp_all_trackHits[0].
id());
712 if (pd && event.getByLabel(pd->inputTag(), hithandle)) {
713 art::fill_ptr_vector(all_hits, hithandle);
716 for (
int i = 0; i < n_recoTrack; i++) {
717 art::Ptr<recob::Track>
track = tracklist[i];
718 std::vector<art::Ptr<recob::Hit>> all_trackHits = track_hits.at(i);
720 double tmpEcomplet = 0;
721 const simb::MCParticle* particle;
722 truthMatcher(clockData, all_hits, all_trackHits, particle, tmpEfrac, tmpEcomplet);
723 if (!particle)
continue;
727 if (tmpEcomplet > Ecomplet_lepton) {
728 Ecomplet_lepton = tmpEcomplet;
729 Efrac_lepton = tmpEfrac;
730 trackLength_lepton = track->Length();
731 MClepton_reco = particle;
736 if (tmpEcomplet > Ecomplet_proton) {
737 Ecomplet_proton = tmpEcomplet;
738 Efrac_proton = tmpEfrac;
739 trackLength_proton = track->Length();
740 MCproton_reco = particle;
745 if (tmpEcomplet > Ecomplet_pionplus) {
746 Ecomplet_pionplus = tmpEcomplet;
747 Efrac_pionplus = tmpEfrac;
748 trackLength_pion_plus = track->Length();
749 MCpion_plus_reco = particle;
754 if (tmpEcomplet > Ecomplet_pionminus) {
755 Ecomplet_pionminus = tmpEcomplet;
756 Efrac_pionminus = tmpEfrac;
757 trackLength_pion_minus = track->Length();
758 MCpion_minus_reco = particle;
762 else if ((TMath::Abs(particle->PdgCode()) == 321) && (particle->TrackId() ==
MC_kaonID)) {
764 if (tmpEcomplet > Ecomplet_kaon) {
765 Ecomplet_kaon = tmpEcomplet;
766 Efrac_kaon = tmpEfrac;
767 trackLength_kaon = track->Length();
768 MCkaon_reco = particle;
772 else if ((particle->PdgCode() == -11) && (particle->TrackId() ==
MC_michelID)) {
774 if (tmpEcomplet > Ecomplet_michel) {
775 Ecomplet_michel = tmpEcomplet;
776 Efrac_michel = tmpEfrac;
777 trackLength_michel = track->Length();
778 MCmichel_reco = particle;
783 double Reco_LengthRes = truth_lengthLepton - trackLength_lepton;
784 double Reco_LengthResProton = proton_length - trackLength_proton;
785 double Reco_LengthResPionPlus = pion_plus_length - trackLength_pion_plus;
786 double Reco_LengthResPionMinus = pion_minus_length - trackLength_pion_minus;
788 if (MClepton_reco && MClepton) {
799 if (MCproton_reco && MCproton) {
808 if (MCpion_plus_reco && MCpion_plus) {
817 if (MCpion_minus_reco && MCpion_minus) {
826 if (MCkaon_reco && MCkaon) {
838 if (MClepton_reco && MClepton) {
845 if (MCkaon_reco && MCkaon) {
852 if (MCproton_reco && MCproton) {
859 if (MCpion_plus_reco && MCpion_plus) {
866 if (MCpion_minus_reco && MCpion_minus) {
873 if (MCmichel_reco && MCmichel) {
std::string fTrackModuleLabel
std::string fMCTruthModuleLabel
TH1D * h_michelwtrk_length
process_name use argoneut_mc_hitfinder track
TH1D * h_protonwtrk_length
TH1D * h_trackRes_pion_plus
void truthMatcher(detinfo::DetectorClocksData const &clockData, std::vector< art::Ptr< recob::Hit >> all_hits, std::vector< art::Ptr< recob::Hit >> track_hits, const simb::MCParticle *&MCparticle, double &Efrac, double &Ecomplet)
double MC_lepton_startMomentum[4]
TH1D * h_Efrac_pion_minus
double truthLength(const detinfo::DetectorClocksData &clockData, detinfo::DetectorPropertiesData const &detProp, const simb::MCParticle *MCparticle)
double MC_leading_PionPlusP
TH1D * h_Ecomplet_pion_plus
int MC_leading_PionMinusID
TH1D * h_trackRes_pion_minus
int MC_leading_PionPlusID
TH1D * h_pionmwtrk_length
double MC_leading_ProtonP
TH1D * h_pionpwtrk_length
double MC_leading_PionMinusP
TH1D * h_Ecomplet_pion_minus
bool insideFV(double vertex[4])