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])