EnergyCont_help.h

Go to the documentation of this file.

#pragma once

#include "CAFAna/Core/Var.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "StandardRecord/Proxy/SRProxy.h"

namespace ana {
  
  static const float MassNeut = 0.939565;  
  static const float MassProt = 0.938272;

  static const bool DEBUGGING = false;

  // --- Return the true energy or true uncontained energy of the primaries, summed by particle type.
  double CalcPrimE( const caf::SRProxy* sr, int WhPDG, bool GetUncontE) {
    // --- Define my variable to return.
    double MyEnergy = 0.;
    // --- Loop through all primary particles.
    for (unsigned int PrimL=0; PrimL < sr->mc.nu[0].prim.size(); ++PrimL) {
      // --- Continue if not the particle type we want.
      if ( abs(sr->mc.nu[0].prim[PrimL].pdg) != WhPDG ) continue;
      // --- Do I want to get the primary energy, or the uncontained energy?
      if (GetUncontE) {
        // --- Check that uncontained energy is set.
        if (sr->mc.nu[0].prim[PrimL].totEscE>0.0) MyEnergy += sr->mc.nu[0].prim[PrimL].totEscE;
      } else {
        // --- Check that true energy is set.
        if(sr->mc.nu[0].prim[PrimL].p.E<0.0) continue;
        // --- Grab kinetic energy for neutrons and protons, and total energy for all other particle types. 
        double E = 0.;
        if (sr->mc.nu[0].prim[PrimL].pdg == 2112) {
          E = sr->mc.nu[0].prim[PrimL].p.E - MassNeut;
        } else if (sr->mc.nu[0].prim[PrimL].pdg == 2212) {
          E = sr->mc.nu[0].prim[PrimL].p.E - MassProt;
        } else if (abs(sr->mc.nu[0].prim[PrimL].pdg) == 211) {
          E = sr->mc.nu[0].prim[PrimL].p.E;
        } else if (sr->mc.nu[0].prim[PrimL].pdg == 111) {
          E = sr->mc.nu[0].prim[PrimL].p.E;
        } else if (abs(sr->mc.nu[0].prim[PrimL].pdg) == 13) {
          E = sr->mc.nu[0].prim[PrimL].p.E;
        } else if (sr->mc.nu[0].prim[PrimL].pdg == 22 ) {
          E = sr->mc.nu[0].prim[PrimL].p.E;
        } else if (abs(sr->mc.nu[0].prim[PrimL].pdg) == 321) {
          E = sr->mc.nu[0].prim[PrimL].p.E;
        } else if (abs(sr->mc.nu[0].prim[PrimL].pdg) == 11) {
          E = sr->mc.nu[0].prim[PrimL].p.E;
        }
        // --- Sum energy for the particle type
        MyEnergy += E;
      }
    } // Loop through prim particles.
    if (DEBUGGING) std::cout << "\tThe total " << (GetUncontE==true?"uncontained":"kinetic") << " energy from " << WhPDG << " was " << MyEnergy << std::endl;
    // --- Return the value we have just worked out.
    return MyEnergy;
  }

  // ----- Neutrino energy
  const Var kNeutrinoEn([](const caf::SRProxy *sr) {
      float ThisEn = 0.;
      ThisEn = sr->mc.nu[0].E;
      return ThisEn;
    });

  // ---- Find out the total uncontained energy
  const Var kTotUncontEn([](const caf::SRProxy *sr) {
      float TotEn = 0.0;
      for (unsigned int PrimL=0; PrimL < sr->mc.nu[0].prim.size(); ++PrimL) {
        TotEn += sr->mc.nu[0].prim[PrimL].totEscE;
      }
      // --- want to know the total uncontained energy and run / subrun / slice / event number 
      // so we can go look at the events in these "bumps" in the evd
      std::cout << "Energy " << TotEn << " Run " << sr->hdr.run << ", SubRun " << sr->hdr.subrun << ", Cycle " << sr->hdr.cycle << ", Batch " << sr->hdr.batch << ", Event " << sr->hdr.evt << ", Slice " << sr->hdr.subevt << std::endl;
      // --- Now return our final value...
      return TotEn;
    });

  // ----- Muon energy
  // Total primary energy
  const Var kMuonEn([](const caf::SRProxy* sr) {
      float ThisEn = 0.;
      ThisEn = CalcPrimE( sr, 13, false );
      return ThisEn;
    });
  // Uncontained energy
  const Var kUncontMuonEn([](const caf::SRProxy* sr) {
      float ThisEn = 0.;
      ThisEn = CalcPrimE( sr, 13, true );
      return ThisEn;
    });

  // ----- Neutron energy
  // Total primary energy
  const Var kNeutronEn([](const caf::SRProxy* sr) {
      float ThisEn = 0.;
      ThisEn = CalcPrimE( sr, 2112, false );
      return ThisEn;
    });
  // Uncontained energy
  const Var kUncontNeutronEn([](const caf::SRProxy* sr) {
      float ThisEn = 0.;
      ThisEn = CalcPrimE( sr, 2112, true );
      return ThisEn;
    });

  // ----- Proton energy
  // Total primary energy
  const Var kProtonEn([](const caf::SRProxy* sr) {
      float ThisEn = 0.;
      ThisEn = CalcPrimE( sr, 2212, false );
      return ThisEn;
    });
  // Uncontained energy
  const Var kUncontProtonEn([](const caf::SRProxy* sr) {
      float ThisEn = 0.;
      ThisEn = CalcPrimE( sr, 2212, true );
      return ThisEn;
    });


  // ----- Pion energy
  // Total primary energy
  const Var kPionEn([](const caf::SRProxy* sr) {
      float ThisEn = 0.;
      ThisEn = CalcPrimE( sr, 211, false );
      return ThisEn;
    });
  // Uncontained energy
  const Var kUncontPionEn([](const caf::SRProxy* sr) {
      float ThisEn = 0.;
      ThisEn = CalcPrimE( sr, 211, true );
      return ThisEn;
    });

  // ----- Pi0 energy
  // Total primary energy
  const Var kPi0En([](const caf::SRProxy* sr) {
      float ThisEn = 0.;
      ThisEn = CalcPrimE( sr, 111, false );
      return ThisEn;
    });
  // Uncontained energy
  const Var kUncontPi0En([](const caf::SRProxy* sr) {
      float ThisEn = 0.;
      ThisEn = CalcPrimE( sr, 111, true );
      return ThisEn;
    });

  // ----- Photon energy
  // Total primary energy
  const Var kPhotonEn([](const caf::SRProxy* sr) {
      float ThisEn = 0.;
      ThisEn = CalcPrimE( sr, 22, false );
      return ThisEn;
    });
  // Uncontained energy
  const Var kUncontPhotonEn([](const caf::SRProxy* sr) {
      float ThisEn = 0.;
      ThisEn = CalcPrimE( sr, 22, true );
      return ThisEn;
    });

  // ----- Kaon energy
  // Total primary energy
  const Var kKaonEn([](const caf::SRProxy* sr) {
      float ThisEn = 0.;
      ThisEn = CalcPrimE( sr, 321, false );
      return ThisEn;
    });
  // Uncontained energy
  const Var kUncontKaonEn([](const caf::SRProxy* sr) {
      float ThisEn = 0.;
      ThisEn = CalcPrimE( sr, 321, true );
      return ThisEn;
    });

  // ----- Electron energy
  // Total primary energy
  const Var kElectronEn([](const caf::SRProxy* sr) {
      float ThisEn = 0.;
      ThisEn = CalcPrimE( sr, 11, false ); 
      return ThisEn;
    });
  // Uncontained energy
  const Var kUncontElectronEn([](const caf::SRProxy* sr) {
      float ThisEn = 0.;
      ThisEn = CalcPrimE( sr, 11, true );
      return ThisEn;
    });


// ------------------------------------------------------------
// A few extra Vars that make use of the previous ones...

    // nu E - energy from each selected particle type
    const Var kNuMinusAll([](const caf::SRProxy* sr) {
      return kNeutrinoEn(sr) - kMuonEn(sr) - kNeutronEn(sr) - kProtonEn(sr) - kPionEn(sr) - kPi0En(sr) - kPhotonEn(sr) - kKaonEn(sr) - kElectronEn(sr);
    });

    // total uncontained E - uncontianed E from each selected particle type
    const Var kUncontMiscEn([](const caf::SRProxy* sr) {
      return kTotUncontEn(sr) - kUncontMuonEn(sr) - kUncontNeutronEn(sr) - kUncontProtonEn(sr) - kUncontPionEn(sr) - kUncontPi0En(sr) - kUncontPhotonEn(sr) - kUncontKaonEn(sr) - kUncontElectronEn(sr);
    });





// ------------------------------------------------------------


  // ----- Energy fraction
  // calculate particle type E / true nu E
  // or
  // uncontained E per particle type / total uncontained E
  // or
  // uncontained E per particle type / true nu E
  double CalcFracE( const caf::SRProxy* sr, int WhPDG, bool GetUncontE, bool UncontOverTotal=false) {
    // calculate the total energy of the particle type
    double E = CalcPrimE(sr, WhPDG, GetUncontE); 
    // is it uncontained?
    if(GetUncontE) {
      // uncont. E / true nu E
      if(UncontOverTotal){
        if(sr->mc.nu[0].E<=0.0) return -5.0; // avoid nan
        else return (E/sr->mc.nu[0].E);
      }
      // uncont. E / total uncont. E
      else{
        if(kTotUncontEn(sr)<=0.0) return -5.0; // avoid nan
        else return (E/kTotUncontEn(sr));
      }
    }

    // particle type E / true nu E
    else {
      if(sr->mc.nu[0].E<=0.0) return -5.0; // avoid nan
      else return (E/sr->mc.nu[0].E);
    }
  }

  //  total uncont. E / true nu E
  const Var kTotUncontEnFrac([](const caf::SRProxy* sr) {
    if(sr->mc.nu[0].E<=0.0) return -5.0;
    else return (kTotUncontEn(sr)/sr->mc.nu[0].E);
  });

  // ----- Muon
  // energy fraction
  const Var kMuonEnFrac([](const caf::SRProxy* sr) {
    return CalcFracE(sr, 13, false);
  }); 
  // uncont./tot. uncont.
  const Var kUncontMuonEnFrac([](const caf::SRProxy* sr) {
    return CalcFracE(sr, 13, true);
  }); 
  // uncont./true nu E
  const Var kUncontMuonEnFracTot([](const caf::SRProxy* sr) {
    return CalcFracE(sr, 13, true, true);
  });

  // ----- Neutron
  // energy fraction
  const Var kNeutronEnFrac([](const caf::SRProxy* sr) {
    return CalcFracE(sr, 2112, false);  
  });
  // uncont./tot. uncont.
  const Var kUncontNeutronEnFrac([](const caf::SRProxy* sr) {
    return CalcFracE(sr, 2112, true);  
  });
  // uncont./true nu E
  const Var kUncontNeutronEnFracTot([](const caf::SRProxy* sr) {
    return CalcFracE(sr, 2112, true, true);
  });

  // ----- Proton
  // energy fraction
  const Var kProtonEnFrac([](const caf::SRProxy* sr) {
    return CalcFracE(sr, 2212, false);
  });
  // uncont./tot. uncont.
  const Var kUncontProtonEnFrac([](const caf::SRProxy* sr) {
    return CalcFracE(sr, 2212, true);
  });
  // uncont./true nu E
  const Var kUncontProtonEnFracTot([](const caf::SRProxy* sr) {
    return CalcFracE(sr, 2212, true, true);
  });

  // ----- Charged Pions
  // energy fraction
  const Var kPionEnFrac([](const caf::SRProxy *sr) {
    return CalcFracE(sr, 211, false);
  });
  // uncont./tot. uncont.
  const Var kUncontPionEnFrac([](const caf::SRProxy *sr) {
    return CalcFracE(sr, 211, true);
  });
  // uncont./true nu E
  const Var kUncontPionEnFracTot([](const caf::SRProxy *sr) {
    return CalcFracE(sr, 211, true, true);
  });

  // ----- Neutral Pions
  // energy fraction
  const Var kPi0EnFrac([](const caf::SRProxy *sr) {
    return CalcFracE(sr, 111, false);
  });
  // uncont./tot. uncont.
  const Var kUncontPi0EnFrac([](const caf::SRProxy *sr) {
    return CalcFracE(sr, 111, true);
  });
  // uncont./true nu E
  const Var kUncontPi0EnFracTot([](const caf::SRProxy *sr) {
    return CalcFracE(sr, 111, true, true);
  });

  // ----- Gammas
  // energy fraction
  const Var kPhotonEnFrac([](const caf::SRProxy *sr) {
    return CalcFracE(sr, 22, false);
  });
  // uncont./tot. uncont.
  const Var kUncontPhotonEnFrac([](const caf::SRProxy *sr) {
    return CalcFracE(sr, 22, true);
  });
  // uncont./true nu E
  const Var kUncontPhotonEnFracTot([](const caf::SRProxy *sr) {
    return CalcFracE(sr, 22, true, true);
  });

  // ----- Kaons
  // energy fraction
  const Var kKaonEnFrac([](const caf::SRProxy *sr) {
    return CalcFracE(sr, 321, false);
  });
  // uncont./tot. uncont.
  const Var kUncontKaonEnFrac([](const caf::SRProxy *sr) {
    return CalcFracE(sr, 321, true);
  });
  // uncont./true nu E
  const Var kUncontKaonEnFracTot([](const caf::SRProxy *sr) {
    return CalcFracE(sr, 321, true, true);
  });

  // ----- Electrons
  // energy fraction
  const Var kElectronEnFrac([](const caf::SRProxy *sr) {
    return CalcFracE(sr, 11, false);
  });
  // uncont./tot. uncont.
  const Var kUncontElectronEnFrac([](const caf::SRProxy *sr) {
    return CalcFracE(sr, 11, true);
  });
  // uncont./true nu E
  const Var kUncontElectronEnFracTot([](const caf::SRProxy *sr) {
    return CalcFracE(sr, 11, true, true);
  });

  // ----- Find out the number of neutrons
  const Var kNNeutrons([](const caf::SRProxy* sr) {
      float NPar = 0.;
      // --- Write out some basic event information.
      if (DEBUGGING) {
        std::cout << "\nLooking at event " << sr->hdr.evt << ", there are " << sr->mc.nu.size() << ", the 0th neutrino is a " << sr->mc.nu[0].pdg << ", with Energy " << sr->mc.nu[0].E 
                  << ". There were a total of " << sr->mc.nu[0].prim.size() << " other primaries associated with that neutrino, " << sr->mc.nu[0].nneutron << " were neutrons."
                  << std::endl;
        if ( kIsDytmanMEC(sr) ) std::cout << "  The neutrino was a MEC" << std::endl;
        else if ( kIsDIS(sr) )  std::cout << "  The neutrino was a DIS" << std::endl;
        else if ( kIsRes(sr) )  std::cout << "  The neutrino was a Resonance" << std::endl;
        else if ( kIsCoh(sr) )  std::cout << "  The neutrino was a Coherent" << std::endl;
        float TotEn = 0.0;
        for (unsigned int PrimL=0; PrimL < sr->mc.nu[0].prim.size(); ++PrimL) {
          std::cout << "\t\tLooking at Prim " << PrimL << " of " << sr->mc.nu[0].prim.size() << ", it was a " << sr->mc.nu[0].prim[PrimL].pdg 
                    << ", VisE " << sr->mc.nu[0].prim[PrimL].visE << ", VisEInslc " << sr->mc.nu[0].prim[PrimL].visEinslc << ", E0 " << sr->mc.nu[0].prim[PrimL].p.E
                    << ", enteringE " << sr->mc.nu[0].prim[PrimL].enteringE << ", totEscE " << sr->mc.nu[0].prim[PrimL].totEscE
                    << std::endl;
          TotEn += sr->mc.nu[0].prim[PrimL].totEscE;
        }
        std::cout << "\tThe total uncontained energy was " << TotEn << std::endl;
        // --- Return the number of neutrons
        NPar = sr->mc.nu[0].nneutron;
      }
      return NPar;
    });

  // ----- Find out the number of protons
  const Var kNProtons([](const caf::SRProxy* sr) {
      float NPar = 0.;
      NPar = sr->mc.nu[0].nproton;
      // --- Return the number of protons
      return NPar;
    });

  // ----- Find out the number of pions
  const Var kNPions([](const caf::SRProxy* sr) {
      float NPar = 0.;
      NPar = sr->mc.nu[0].npiplus + sr->mc.nu[0].npiminus;
      // --- Return the number of pions
      return NPar;
    });

  // ----- Find out the number of pi0s
  const Var kNPi0s([](const caf::SRProxy* sr) {
      float NPar = 0.;
      NPar = sr->mc.nu[0].npizero;
      // --- Return the number of pi0s
      return NPar;
    });

} // namespace