NuShowerEff Class Reference

Inheritance diagram for NuShowerEff:

Public Member Functions
	NuShowerEff (fhicl::ParameterSet const &p)
	NuShowerEff (NuShowerEff const &)=delete
	NuShowerEff (NuShowerEff &&)=delete
NuShowerEff &	operator= (NuShowerEff const &)=delete
NuShowerEff &	operator= (NuShowerEff &&)=delete
Public Member Functions inherited from art::EDAnalyzer
	EDAnalyzer (fhicl::ParameterSet const &pset)
template<typename Config >
	EDAnalyzer (Table< Config > const &config)
std::string	workerType () const
Public Member Functions inherited from art::detail::Analyzer
virtual	~Analyzer () noexcept
	Analyzer (fhicl::ParameterSet const &pset)
template<typename Config >
	Analyzer (Table< Config > const &config)
void	doBeginJob (SharedResources const &resources)
void	doEndJob ()
void	doRespondToOpenInputFile (FileBlock const &fb)
void	doRespondToCloseInputFile (FileBlock const &fb)
void	doRespondToOpenOutputFiles (FileBlock const &fb)
void	doRespondToCloseOutputFiles (FileBlock const &fb)
bool	doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
bool	doEndRun (RunPrincipal &rp, ModuleContext const &mc)
bool	doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
bool	doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
bool	doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
Public Member Functions inherited from art::Observer
	~Observer () noexcept
	Observer (Observer const &)=delete
	Observer (Observer &&)=delete
Observer &	operator= (Observer const &)=delete
Observer &	operator= (Observer &&)=delete
void	registerProducts (ProductDescriptions &, ModuleDescription const &)
void	fillDescriptions (ModuleDescription const &)
fhicl::ParameterSetID	selectorConfig () const
Public Member Functions inherited from art::ModuleBase
virtual	~ModuleBase () noexcept
	ModuleBase ()
ModuleDescription const &	moduleDescription () const
void	setModuleDescription (ModuleDescription const &)
std::array< std::vector< ProductInfo >, NumBranchTypes > const &	getConsumables () const
void	sortConsumables (std::string const &current_process_name)
template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)
template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

Private Member Functions
void	analyze (art::Event const &e) override
void	beginJob () override
void	endJob () override
void	beginRun (art::Run const &run) override
void	doEfficiencies ()
bool	insideFV (double vertex[4])
void	reset ()

Private Attributes
std::string	fMCTruthModuleLabel
std::string	fHitModuleLabel
std::string	fShowerModuleLabel
std::string	fTruthMatchDataModuleLabel
int	fNeutrinoPDGcode
int	fLeptonPDGcode
bool	fSaveMCTree
bool	fHitShowerThroughPFParticle
double	fMinPurity
double	fMinCompleteness
float	fFidVolCutX
float	fFidVolCutY
float	fFidVolCutZ
float	fFidVolXmin
float	fFidVolXmax
float	fFidVolYmin
float	fFidVolYmax
float	fFidVolZmin
float	fFidVolZmax
int	Event
int	Run
int	SubRun
int	MC_incoming_PDG
int	MC_lepton_PDG
int	MC_isCC
int	MC_channel
int	MC_target
double	MC_Q2
double	MC_W
double	MC_vertex [4]
double	MC_incoming_P [4]
double	MC_lepton_startMomentum [4]
double	MC_lepton_endMomentum [40]
double	MC_lepton_startXYZT [4]
double	MC_lepton_endXYZT [4]
double	MC_lepton_theta
int	MC_leptonID
int	MC_LeptonTrack
double	MC_incoming_E
double	MC_lepton_startEnergy
int	n_recoShowers
double	sh_direction_X [MAX_SHOWERS]
double	sh_direction_Y [MAX_SHOWERS]
double	sh_direction_Z [MAX_SHOWERS]
double	sh_start_X [MAX_SHOWERS]
double	sh_start_Y [MAX_SHOWERS]
double	sh_start_Z [MAX_SHOWERS]
double	sh_length [MAX_SHOWERS]
double	sh_ehit_Q [MAX_SHOWERS]
int	sh_TrackId [MAX_SHOWERS]
int	sh_hasPrimary_e [MAX_SHOWERS]
double	sh_allhit_Q [MAX_SHOWERS]
double	sh_purity [MAX_SHOWERS]
double	sh_completeness [MAX_SHOWERS]
double	esh_1_purity
double	esh_1_completeness
double	esh_each_purity [MAX_SHOWERS]
double	esh_each_completeness [MAX_SHOWERS]
double	esh_purity
double	esh_completeness
int	count_primary_e_in_Event
TTree *	fEventTree
TH1D *	h_Ev_den
TH1D *	h_Ev_num
TH1D *	h_Ee_den
TH1D *	h_Ee_num
TEfficiency *	h_Eff_Ev = 0
TEfficiency *	h_Eff_Ee = 0

Additional Inherited Members
Public Types inherited from art::EDAnalyzer
using	WorkerType = WorkerT< EDAnalyzer >
using	ModuleType = EDAnalyzer
Protected Member Functions inherited from art::Observer
std::string const &	processName () const
bool	wantAllEvents () const noexcept
bool	wantEvent (ScheduleID id, Event const &e) const
Handle< TriggerResults >	getTriggerResults (Event const &e) const
	Observer (fhicl::ParameterSet const &config)
	Observer (std::vector< std::string > const &select_paths, std::vector< std::string > const &reject_paths, fhicl::ParameterSet const &config)
Protected Member Functions inherited from art::ModuleBase
ConsumesCollector &	consumesCollector ()
template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)
template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)
template<typename T , BranchType = InEvent>
void	consumesMany ()
template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)
template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)
template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

Detailed Description

Definition at line 48 of file NuShowerEff_module.cc.

Constructor & Destructor Documentation

NuShowerEff::NuShowerEff ( fhicl::ParameterSet const &  p )

explicit

Definition at line 165 of file NuShowerEff_module.cc.

  :
  EDAnalyzer(p), // ,
 // More initializers here.
  //fMCTruthModuleLabel (p.get< std::string >("MCTruthModuleLabel", "generator")),
  fMCTruthModuleLabel (p.get< std::string >("MCTruthModuleLabel")),// get parameter from fcl file
  fHitModuleLabel     (p.get< std::string >("HitModuleLabel")),
  fShowerModuleLabel  (p.get< std::string >("ShowerModuleLabel")),
  fTruthMatchDataModuleLabel (p.get< std::string >("TruthMatchDataModuleLabel")),
  fNeutrinoPDGcode    (p.get<int>("NeutrinoPDGcode")),
  fLeptonPDGcode      (p.get<int>("LeptonPDGcode")),
  fSaveMCTree         (p.get<bool>("SaveMCTree")),
  fHitShowerThroughPFParticle (p.get<bool>("HitShowerThroughPFParticle")),
  fMinPurity          (p.get<double>("MinPurity")),
  fMinCompleteness    (p.get<double>("MinCompleteness")),
  fFidVolCutX         (p.get<float>("FidVolCutX")),
  fFidVolCutY         (p.get<float>("FidVolCutY")),
  fFidVolCutZ         (p.get<float>("FidVolCutZ"))
{
  //cout << "\n===== Please refer the fchicl for the values of preset parameters ====\n" << endl;
}

NuShowerEff::NuShowerEff ( NuShowerEff const &  )

delete

NuShowerEff::NuShowerEff ( NuShowerEff &&  )

delete

Member Function Documentation

void NuShowerEff::analyze ( art::Event const &  e )

overrideprivatevirtual

Implements art::EDAnalyzer.

Definition at line 284 of file NuShowerEff_module.cc.

{
  // Implementation of required member function here.
  //cout << "\n===== function: analyze() =====\n" << endl;

  reset(); // for some variables

  Event = e.id().event();
  Run = e.run();
  SubRun = e.subRun();
  //cout << "Event information:" << endl;
  //cout << "\tEvent: " << Event << endl;
  //cout << "\tRun: " << Run << endl;
  //cout << "\tSubRun: " << SubRun << endl;


  // -------- find Geant4 TrackId that corresponds to e+/e- from neutrino interaction ----------
  // MCTruth: Generator
  art::Handle<std::vector<simb::MCTruth>> MCtruthHandle;
  e.getByLabel(fMCTruthModuleLabel, MCtruthHandle);
  std::vector<art::Ptr<simb::MCTruth>> MCtruthlist;
  art::fill_ptr_vector(MCtruthlist, MCtruthHandle);
  art::Ptr<simb::MCTruth> MCtruth;
  // MC (neutrino) interaction
  int MCinteractions = MCtruthlist.size();
  //cout << "\nMCinteractions: " << MCinteractions << endl;
  for (int i=0; i<MCinteractions; i++){
    MCtruth = MCtruthlist[i];
    if ( MCtruth->NeutrinoSet() ) {  // NeutrinoSet(): whether the neutrino information has been set
      simb::MCNeutrino nu = MCtruth->GetNeutrino();// GetNeutrino(): reference to neutrino info
      if( nu.CCNC()==0 ) MC_isCC = 1;
      else if (nu.CCNC()==1) MC_isCC = 0;
      simb::MCParticle neutrino = nu.Nu(); // Nu(): the incoming neutrino
      MC_target = nu.Target(); // Target(): nuclear target, as PDG code
      MC_incoming_PDG = std::abs(nu.Nu().PdgCode());// here not use std::abs()
      MC_Q2 = nu.QSqr();// QSqr(): momentum transfer Q^2, in GeV^2
      MC_channel = nu.InteractionType();// 1001: CCQE
      MC_W = nu.W(); // W(): hadronic invariant mass
      const TLorentzVector& nu_momentum = nu.Nu().Momentum(0);
      nu_momentum.GetXYZT(MC_incoming_P);
      const TLorentzVector& vertex =neutrino.Position(0);
      vertex.GetXYZT(MC_vertex);
      simb::MCParticle lepton = nu.Lepton();// Lepton(): the outgoing lepton
      MC_lepton_PDG = lepton.PdgCode();

      MC_incoming_E = MC_incoming_P[3];

      //cout << "\tMCinteraction: " << i              << "\n\t"
      //     << "neutrino PDG: "  << MC_incoming_PDG  << "\n\t"
      //     << "MC_lepton_PDG: " << MC_lepton_PDG    << "\n\t"
      //     << "MC_channel: "    << MC_channel       << "\n\t"
      //     << "incoming E: "    << MC_incoming_P[3] << "\n\t"
      //     << "MC_vertex: " << MC_vertex[0] << " , " << MC_vertex[1] << " , " <<MC_vertex[2] << " , " <<MC_vertex[3] << endl;
    }
    // MCTruth Generator Particles
    int nParticles =  MCtruthlist[0]->NParticles();
    //cout << "\n\tNparticles: " << MCtruth->NParticles() <<endl;
    for (int i=0; i<nParticles; i++){
      simb::MCParticle pi = MCtruthlist[0]->GetParticle(i);
      //cout << "\tparticle: " << i << "\n\t\t"
           //<< "Pdgcode: " << pi.PdgCode() << "; Mother: " << pi.Mother() << "; TrackId: " << pi.TrackId() << endl;
    // Mother(): mother = -1 means that this particle has no mother
    // TrackId(): same as the index in the MCParticleList
    }
  }

  // Geant4: MCParticle -> lepton (e)
  // Note: generator level MCPartilceList is different from Geant4 level MCParticleList.  MCParticleList(Geant4) contains all particles in MCParticleList(generator) but their the indexes (TrackIds) are different.
  simb::MCParticle *MClepton = NULL; //Geant4 level
  art::ServiceHandle<cheat::ParticleInventoryService const> pi_serv;
  const sim::ParticleList& plist = pi_serv->ParticleList();
  simb::MCParticle *particle=0;

  for (sim::ParticleList::const_iterator ipar = plist.begin(); ipar!=plist.end();++ipar){
    particle = ipar->second; // first is index(TrackId), second is value (point address)

    auto & truth = pi_serv->ParticleToMCTruth_P(particle);// beam neutrino only
    if ( truth->Origin()==simb::kBeamNeutrino && std::abs(particle->PdgCode()) == fLeptonPDGcode &&  particle->Mother()==0){ // primary lepton; Mother() = 0 means e^{-} for v+n=e^{-}+p
      const TLorentzVector& lepton_momentum =particle->Momentum(0);
      const TLorentzVector& lepton_position =particle->Position(0);
      const TLorentzVector& lepton_positionEnd   = particle->EndPosition();
      const TLorentzVector& lepton_momentumEnd   = particle->EndMomentum();
      lepton_momentum.GetXYZT(MC_lepton_startMomentum);
      lepton_position.GetXYZT(MC_lepton_startXYZT);
      lepton_positionEnd.GetXYZT(MC_lepton_endXYZT);
      lepton_momentumEnd.GetXYZT(MC_lepton_endMomentum);
      MC_leptonID = particle->TrackId();

      MC_lepton_startEnergy = MC_lepton_startMomentum[3];
      //cout << "\nGeant Track ID for electron/positron: " << endl;
      //cout << "\tMClepton PDG: " << particle->PdgCode() <<  " ; MC_leptonID: " << MC_leptonID << endl;
      MClepton = particle;
      //cout << "\tMClepton PDG:" << MClepton->PdgCode() <<endl;
      //cout << "\tipar->first (TrackId): " << ipar->first << endl;
    }
  }

  // check if the interaction is inside the Fiducial Volume
  bool isFiducial = false;
  isFiducial = insideFV( MC_vertex);
  if (isFiducial) {
    //cout <<"\nInfo: Interaction is inside the Fiducial Volume.\n" << endl;
  }
  else {
    //cout << "\n********Interaction is NOT inside the Fiducial Volume. RETURN**********" << endl;
    return;
  }

  if (MC_isCC && (fNeutrinoPDGcode == std::abs(MC_incoming_PDG))) {
    if (MClepton){
      h_Ev_den->Fill(MC_incoming_P[3]);
      h_Ee_den->Fill(MC_lepton_startMomentum[3]);
    }
  }

  //if (MC_isCC && (fNeutrinoPDGcode == std::abs(MC_incoming_PDG)) && MC_incoming_E < 0.5) {
    // cout << "\n------output CC event info for neutrino energy less than 0.5 GeV ------\n" << endl;
   // cout << "\tEvent   : " << Event  << "\n"
   //   << "\tRun     : " << Run    << "\n"
   //   << "\tSubRun  : " << SubRun << "\n"
    //  << "\tMC_incoming_E: " << MC_incoming_E << endl;
 // }

  // recob::Hit
  // ---- build a map for total hit charges corresponding to MC_leptonID (Geant4) ----

  art::Handle<std::vector<recob::Hit>> hitHandle;
  std::vector<art::Ptr<recob::Hit>> all_hits;
  if(e.getByLabel(fHitModuleLabel,hitHandle)){
    art::fill_ptr_vector(all_hits, hitHandle);
  }
  //cout << "\nTotal hits:" << endl;
  //cout << "\tall_hits.size(): " << all_hits.size() << endl;

  double ehit_total =0.0;

  art::FindManyP<simb::MCParticle,anab::BackTrackerHitMatchingData> fmhitmc(hitHandle, e, fTruthMatchDataModuleLabel);// need #include "lardataobj/AnalysisBase/BackTrackerMatchingData.h"

  std::map<int,double> all_hits_trk_Q;//Q for charge: Integral()
  for (size_t i=0; i <  all_hits.size(); ++i) {
    art::Ptr<recob::Hit> hit = all_hits[i];
    auto particles = fmhitmc.at(hit.key());// particles here is a pointer. A hit may come from multiple particles.
    auto hitmatch = fmhitmc.data(hit.key());// metadata
    double maxenergy = -1e9;
    int hit_TrackId = 0;
    for (size_t j = 0; j < particles.size(); ++j){
      if (!particles[j]) continue;
      if (!pi_serv->TrackIdToMotherParticle_P(particles[j]->TrackId())) continue;
      size_t trkid = (pi_serv->TrackIdToMotherParticle_P(particles[j]->TrackId()))->TrackId();

      if ( (hitmatch[j]->energy) > maxenergy ){
        maxenergy = hitmatch[j]->energy;
        hit_TrackId = trkid;
      }
   }
   if (hit_TrackId == MC_leptonID) ehit_total += hit->Integral();
   all_hits_trk_Q[hit_TrackId] += hit->Integral(); // Integral(): the integral under the calibrated signal waveform of the hit, in tick x ADC units
  }
  //cout << "....ehit_total: " << ehit_total << endl;
  //cout << "\tall_hits_trk_Q.size(): " << all_hits_trk_Q.size() << endl;


  //--------- Loop over all showers: find the associated hits for each shower ------------
  const simb::MCParticle *MClepton_reco = NULL;

  double temp_sh_ehit_Q = 0.0;
  double temp_sh_allHit_Q = 0.0;
  int temp_sh_TrackId = -999;
  count_primary_e_in_Event = 0;

  art::Handle<std::vector<recob::Shower>> showerHandle;
  std::vector<art::Ptr<recob::Shower>> showerlist;
  if(e.getByLabel(fShowerModuleLabel,showerHandle)){
    art::fill_ptr_vector(showerlist, showerHandle);
  }
  n_recoShowers= showerlist.size();
  //cout << "\nRecon Showers: " << endl;
  //cout << "\tn_recoShowers: " << n_recoShowers << endl;
  art::FindManyP<recob::Hit> sh_hitsAll(showerHandle, e, fShowerModuleLabel);

  //std::vector<std::map<int,double>> showers_hits_trk_Q;
  std::map<int,double> showers_trk_Q;
  for (int i=0; i<n_recoShowers;i++){ // loop over showers
    //const simb::MCParticle *particle;
    sh_hasPrimary_e[i] = 0;

    std::map<int,double> sh_hits_trk_Q;//Q for charge: Integral()
    sh_hits_trk_Q.clear();

    art::Ptr<recob::Shower> shower = showerlist[i];
    sh_direction_X[i] = shower->Direction().X(); // Direction(): direction cosines at start of the shower
    sh_direction_Y[i] = shower->Direction().Y();
    sh_direction_Z[i] = shower->Direction().Z();
    sh_start_X[i] = shower->ShowerStart().X();
    sh_start_Y[i] = shower->ShowerStart().Y();
    sh_start_Z[i] = shower->ShowerStart().Z();
    sh_length[i] = shower->Length(); // shower length in [cm]
    //cout << "\tInfo of shower " << i << "\n\t\t"
    //<< "direction (cosines): " << sh_direction_X[i] << ", " << sh_direction_Y[i] << ", " << sh_start_Z[i] << "\n\t\t"
    //<< "start position: " << sh_start_X[i] << ", " << sh_start_Y[i] << ", " << sh_start_Z[i] << "\n\t\t"
    //<< "shower length: " << sh_length[i] << endl;


    std::vector<art::Ptr<recob::Hit>> sh_hits;// associated hits for ith shower
    //In mcc8, if we get hits associated with the shower through shower->hits association directly for pandora, the hit list is incomplete. The recommended way of getting hits is through association with pfparticle:
    //shower->pfparticle->clusters->hits
    //----------getting hits through PFParticle (an option here)-------------------
    if (fHitShowerThroughPFParticle) {
      //cout << "\n==== Getting Hits associated with shower THROUGH PFParticle ====\n" << endl;
      //cout << "\nHits in a shower through PFParticle:\n" << endl;

      // PFParticle
      art::Handle<std::vector<recob::PFParticle> > pfpHandle;
      std::vector<art::Ptr<recob::PFParticle> > pfps;
      if (e.getByLabel(fShowerModuleLabel, pfpHandle)) art::fill_ptr_vector(pfps, pfpHandle);
      // Clusters
      art::Handle<std::vector<recob::Cluster> > clusterHandle;
      std::vector<art::Ptr<recob::Cluster> > clusters;
      if (e.getByLabel(fShowerModuleLabel, clusterHandle)) art::fill_ptr_vector(clusters, clusterHandle);
      art::FindManyP<recob::PFParticle> fmps(showerHandle, e, fShowerModuleLabel);// PFParticle in Shower
      art::FindManyP<recob::Cluster> fmcp(pfpHandle, e, fShowerModuleLabel); // Cluster vs. PFParticle
      art::FindManyP<recob::Hit> fmhc(clusterHandle, e, fShowerModuleLabel); // Hit in Shower
      if (fmps.isValid()){
        std::vector<art::Ptr<recob::PFParticle>> pfs = fmps.at(i);
        for (size_t ipf = 0; ipf<pfs.size(); ++ipf){
          if (fmcp.isValid()){
            std::vector<art::Ptr<recob::Cluster>> clus = fmcp.at(pfs[ipf].key());
            for (size_t iclu = 0; iclu<clus.size(); ++iclu){
              if (fmhc.isValid()){
                std::vector<art::Ptr<recob::Hit>> hits = fmhc.at(clus[iclu].key());
                for (size_t ihit = 0; ihit<hits.size(); ++ihit){
                  sh_hits.push_back(hits[ihit]);
                }
              }
            }
          }
        }
      }

     // cout << "\tsh_hits.size(): " << sh_hits.size() << endl;

     // for (size_t k=0;k<sh_hits.size();k++){
     //   art::Ptr<recob::Hit> hit = sh_hits[k];
     //   cout << k << "\thit.key(): " << hit.key() << endl;
     //   cout << k << "\thit->Integral(): " << hit->Integral() << endl;
     // }
    } else {

      // ----- using shower->hit association for getting hits associated with shower-----
      sh_hits = sh_hitsAll.at(i);// associated hits for ith shower (using association of hits and shower)
      //cout << "\t\tsh_hits.size(): " << sh_hits.size() << endl;
    } // we only choose one method for hits associated with shower here.


    for (size_t k=0; k <  sh_hits.size(); ++k) {
      art::Ptr<recob::Hit> hit = sh_hits[k];
      auto particles = fmhitmc.at(hit.key());
      auto hitmatch = fmhitmc.data(hit.key());
      double maxenergy = -1e9;
      int hit_TrackId = 0;
      for (size_t j = 0; j < particles.size(); ++j){
        if (!particles[j]) continue;
        if (!pi_serv->TrackIdToMotherParticle_P(particles[j]->TrackId())) continue;
        size_t trkid = (pi_serv->TrackIdToMotherParticle_P(particles[j]->TrackId()))->TrackId();

        if ( (hitmatch[j]->energy) > maxenergy ){
          maxenergy = hitmatch[j]->energy;
          hit_TrackId = trkid;
        }
      }
      if (hit_TrackId == MC_leptonID) {
        sh_ehit_Q[i] += hit->Integral();
      }
      sh_allhit_Q[i] += hit->Integral();
      sh_hits_trk_Q[hit_TrackId] += hit->Integral();// Q from the same hit_TrackId
    }
    //cout << "\tsh_hits_trk_Q.size(): " << sh_hits_trk_Q.size() << endl;
    //showers_hits_trk_Q.push_back(sh_hits_trk_Q);

    // get TrackId for a shower
    double maxshowerQ = -1.0e9;
    //sh_TrackId[i] = 0;
    for(std::map<int,double>::iterator k = sh_hits_trk_Q.begin(); k != sh_hits_trk_Q.end(); k++) {
      //cout << k->first << "\t;\t" << k->second << endl;
      if (k->second > maxshowerQ) {
        maxshowerQ = k->second;
       sh_TrackId[i] = k->first;//assign a sh_TrackId with TrackId from hit(particle) that contributing largest to the shower.
      }
    }

    //---------------------------------------------------------------------------------
    //cout << "\nRecon Shower: " << i << endl;
    //cout << "\t*****shower primary TrackId: " << sh_TrackId[i] << endl;

    if (sh_TrackId[i] == MC_leptonID && sh_ehit_Q[i] >0) {
      temp_sh_TrackId = sh_TrackId[i];
      sh_hasPrimary_e[i] = 1;
      count_primary_e_in_Event += 1;
      MClepton_reco = pi_serv->TrackIdToParticle_P(sh_TrackId[i]);

      if (sh_allhit_Q[i] >0 && sh_ehit_Q[i] <= sh_allhit_Q[i]){
        sh_purity[i] = sh_ehit_Q[i]/sh_allhit_Q[i];
        //cout << "\t*****shower purity: " << sh_purity[i] << endl;
      } else {
        sh_purity[i] = 0;
      }
      if(ehit_total >0 && sh_ehit_Q[i] <= sh_allhit_Q[i]){
        sh_completeness[i] = sh_ehit_Q[i] / ehit_total;
        //cout << "\t*****shower completeness: " << sh_completeness[i] << endl;
      } else {
        sh_completeness[i] = 0;
      }
      temp_sh_ehit_Q += sh_ehit_Q[i];
      temp_sh_allHit_Q += sh_allhit_Q[i];
    }

    showers_trk_Q[sh_TrackId[i]] += maxshowerQ;
    //cout << "\tsh_TrackId: " << sh_TrackId[i] <<" ; maxshowerQ: " << maxshowerQ << endl;

  } // end: for loop over showers

  // ---------------------------------------------------------------
  if (MClepton_reco && MClepton) {
    if ((temp_sh_TrackId == MC_leptonID) && MC_isCC && (fNeutrinoPDGcode == std::abs(MC_incoming_PDG))) {
      if ((temp_sh_allHit_Q >= temp_sh_ehit_Q) && (temp_sh_ehit_Q > 0.0)) {
        esh_purity = temp_sh_ehit_Q/temp_sh_allHit_Q;
        esh_completeness = temp_sh_ehit_Q/ehit_total;

        if (esh_purity > fMinPurity &&
          esh_completeness > fMinCompleteness) {
          //cout << "\nInfo: fill h_Ev_num ........\n" << endl;
          h_Ev_num->Fill(MC_incoming_P[3]);
          h_Ee_num->Fill(MC_lepton_startMomentum[3]);
        }
      }
    }
  }
  // --------------------------------------------------------------

  if ( (MClepton_reco && MClepton) && MC_isCC && (fNeutrinoPDGcode == std::abs(MC_incoming_PDG))){
    //cout << "\n count_primary_e_in_Event: " << count_primary_e_in_Event << endl;
    for (int j=0; j<count_primary_e_in_Event; j++){
      esh_each_purity[j] = 0.0;
    }

    double temp_esh_1_allhit = -1.0e9;
    int temp_shower_index = -999;
    int temp_esh_index = 0;
    for (int i=0; i<n_recoShowers; i++) {
      if (sh_TrackId[i] == MC_leptonID) {

        // for each electron shower
        if (sh_ehit_Q[i] >0){
          esh_each_purity[temp_esh_index] = sh_purity[i];
          esh_each_completeness[temp_esh_index] = sh_completeness[i];
          temp_esh_index += 1;
        }

        // find largest shower
        if ((sh_allhit_Q[i] > temp_esh_1_allhit) && (sh_ehit_Q[i] > 0.0) ) {
          temp_esh_1_allhit = sh_allhit_Q[i];
          temp_shower_index = i;
        }
      }
    }
    //if (temp_esh_index != count_primary_e_in_Event){
    //  cout << "wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww" << endl;
    //}
    // largest shower with electron contribution
    esh_1_purity = sh_purity[temp_shower_index];
    esh_1_completeness = sh_completeness[temp_shower_index];
  }

  if (count_primary_e_in_Event>0 && MC_isCC && fSaveMCTree) { fEventTree->Fill();}// so far, only save CC event
}

void NuShowerEff::beginJob ( )

overrideprivatevirtual

Reimplemented from art::EDAnalyzer.

Definition at line 188 of file NuShowerEff_module.cc.

                          {
  //cout << "\n===== function: beginJob() ====\n" << endl;

  // Get geometry: Fiducial Volume
  auto const* geo = lar::providerFrom<geo::Geometry>();
  double minx = 1e9; // [cm]
  double maxx = -1e9;
  double miny = 1e9;
  double maxy = -1e9;
  double minz = 1e9;
  double maxz = -1e9;
  //cout << "\nGeometry:\n\tgeo->NTPC(): " << geo->NTPC() << endl;
  for (size_t i = 0; i<geo->NTPC(); ++i){
    double local[3] = {0.,0.,0.};
    double world[3] = {0.,0.,0.};
    const geo::TPCGeo &tpc = geo->TPC(i);
    tpc.LocalToWorld(local,world);
    //cout << "\tlocal: " << local[0] << " ; " << local[1] << " ; " << local[2] << endl;
    //cout << "\tworld: " << world[0] << " ; " << world[1] << " ; " << world[2] << endl;
    //cout << "\tgeo->DetHalfWidth(" << i << "): " << geo->DetHalfWidth(i) << endl;
    //cout << "\tgeo->DetHalfHeight(" << i << "): " << geo->DetHalfHeight(i) << endl;
    //cout << "\tgeo->DetLength(" << i << "): " << geo->DetLength(i) << endl;

    if (minx > world[0] - geo->DetHalfWidth(i))  minx = world[0]-geo->DetHalfWidth(i);
    if (maxx < world[0] + geo->DetHalfWidth(i))  maxx = world[0]+geo->DetHalfWidth(i);
    if (miny > world[1] - geo->DetHalfHeight(i)) miny = world[1]-geo->DetHalfHeight(i);
    if (maxy < world[1] + geo->DetHalfHeight(i)) maxy = world[1]+geo->DetHalfHeight(i);
    if (minz > world[2] - geo->DetLength(i)/2.)  minz = world[2]-geo->DetLength(i)/2.;
    if (maxz < world[2] + geo->DetLength(i)/2.)  maxz = world[2]+geo->DetLength(i)/2.;
  }

  fFidVolXmin = minx + fFidVolCutX;
  fFidVolXmax = maxx - fFidVolCutX;
  fFidVolYmin = miny + fFidVolCutY;
  fFidVolYmax = maxy - fFidVolCutY;
  fFidVolZmin = minz + fFidVolCutZ;
  fFidVolZmax = maxz - fFidVolCutZ;

  //cout << "\nFiducial Volume (length unit: cm):\n"
       //<< "\t" << fFidVolXmin<<"\t< x <\t"<<fFidVolXmax<<"\n"
       //<< "\t" << fFidVolYmin<<"\t< y <\t"<<fFidVolYmax<<"\n"
       //<< "\t" << fFidVolZmin<<"\t< z <\t"<<fFidVolZmax<<"\n";

  art::ServiceHandle<art::TFileService const> tfs;

  double E_bins[21] = {0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4,4.5,5.0,5.5,6.0,7.0,8.0,10.0,12.0,14.0,17.0,20.0,25.0};
//  double theta_bins[44]= { 0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.,16.,17.,18.,19.,20.,22.,24.,26.,28.,30.,32.,34.,36.,38.,40.,42.,44.,46.,48.,50.,55.,60.,65.,70.,75.,80.,85.,90.};
  //  double Pbins[18] ={0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.5,3.0};

  h_Ev_den = tfs->make<TH1D>("h_Ev_den", "Neutrino Energy; Neutrino Energy (GeV); Shower Reconstruction Efficiency", 20, E_bins);
  h_Ev_den->Sumw2();
  h_Ev_num = tfs->make<TH1D>("h_Ev_num","Neutrino Energy; Neutrino Energy (GeV); Shower Reconstruction Efficiency",20,E_bins);
  h_Ev_num->Sumw2();

  h_Ee_den = tfs->make<TH1D>("h_Ee_den","Electron Energy; Electron Energy (GeV); Shower Reconstruction Efficiency",20,E_bins);
  h_Ee_den->Sumw2();
  h_Ee_num = tfs->make<TH1D>("h_Ee_num","Electron Energy; Electron Energy (GeV); Shower Reconstruction Efficiency",20,E_bins);
  h_Ee_num->Sumw2();

  if (fSaveMCTree) {
    fEventTree = new TTree("Event", "Event Tree from Sim & Reco");
    fEventTree->Branch("eventNo", &Event); // only select event in FV
    fEventTree->Branch("runNo", &Run);
    fEventTree->Branch("subRunNo", &SubRun);

    fEventTree->Branch("MC_incoming_E", &MC_incoming_E);
    fEventTree->Branch("MC_lepton_startEnergy", &MC_lepton_startEnergy);

    fEventTree->Branch("n_showers", &n_recoShowers);
    fEventTree->Branch("sh_hasPrimary_e", &sh_hasPrimary_e, "sh_hasPrimary_e[n_showers]/I");
    //fEventTree->Branch("sh_purity", &sh_purity, "sh_purity[n_showers]/D");
    //fEventTree->Branch("sh_completeness", &sh_completeness, "sh_completeness[n_showers]/D");
    fEventTree->Branch("count_primary_e_in_Event", &count_primary_e_in_Event);
    fEventTree->Branch("esh_1_purity", &esh_1_purity);
    fEventTree->Branch("esh_1_completeness", &esh_1_completeness);
    fEventTree->Branch("esh_each_purity", &esh_each_purity, "esh_each_purity[count_primary_e_in_Event]/D");
    fEventTree->Branch("esh_each_completeness", &esh_each_completeness, "esh_each_completeness[count_primary_e_in_Event]/D");
    fEventTree->Branch("esh_purity", &esh_purity);
    fEventTree->Branch("esh_completeness", &esh_completeness);
  }

}

void NuShowerEff::beginRun ( art::Run const &  run )

overrideprivatevirtual

Reimplemented from art::EDAnalyzer.

Definition at line 278 of file NuShowerEff_module.cc.

                                            {
  //cout << "\n===== function: beginRun() =====\n" << endl;
  mf::LogInfo("ShowerEff") << "==== begin run ... ====" << endl;
}

void NuShowerEff::doEfficiencies ( )

private

Definition at line 661 of file NuShowerEff_module.cc.

                                {
  //cout << "\n==== function: doEfficiencies() ====" << endl;

  art::ServiceHandle<art::TFileService const> tfs;

  if(TEfficiency::CheckConsistency(*h_Ev_num,*h_Ev_den)){
    h_Eff_Ev = tfs->make<TEfficiency>(*h_Ev_num,*h_Ev_den);
    TGraphAsymmErrors *grEff_Ev = h_Eff_Ev->CreateGraph();
    grEff_Ev->Write("grEff_Ev");
    h_Eff_Ev->Write("h_Eff_Ev");
  }

  if(TEfficiency::CheckConsistency(*h_Ee_num,*h_Ee_den)){
    h_Eff_Ee = tfs->make<TEfficiency>(*h_Ee_num,*h_Ee_den);
    TGraphAsymmErrors *grEff_Ee = h_Eff_Ee->CreateGraph();
    grEff_Ee->Write("grEff_Ee");
    h_Eff_Ee->Write("h_Eff_Ee");
  }

}

void NuShowerEff::endJob ( )

overrideprivatevirtual

Reimplemented from art::EDAnalyzer.

Definition at line 272 of file NuShowerEff_module.cc.

                        {
  //cout << "\n===== function: endJob() =====\n" << endl;
  doEfficiencies();
}

bool NuShowerEff::insideFV ( double  vertex[4] )

private

Definition at line 683 of file NuShowerEff_module.cc.

                                            {
  //cout << "\n==== function: insideFV() ====" << endl;
  double x = vertex[0];
  double y = vertex[1];
  double z = vertex[2];

  if (x>fFidVolXmin && x<fFidVolXmax&&
      y>fFidVolYmin && y<fFidVolYmax&&
      z>fFidVolZmin && z<fFidVolZmax)
    {return true;}
  else
    {return false;}
}

NuShowerEff& NuShowerEff::operator= ( NuShowerEff const &  )

delete

NuShowerEff& NuShowerEff::operator= ( NuShowerEff &&  )

delete

void NuShowerEff::reset ( )

private

Definition at line 698 of file NuShowerEff_module.cc.

                       {
  //cout << "\n===== function: reset() =====\n" << endl;

  MC_incoming_PDG = -999;
  MC_lepton_PDG =-999;
  MC_isCC =-999;
  MC_channel =-999;
  MC_target =-999;
  MC_Q2 =-999.0;
  MC_W =-999.0;
  MC_lepton_theta = -999.0;
  MC_leptonID = -999;
  MC_LeptonTrack = -999;

  for(int i=0; i<MAX_SHOWERS; i++){
    sh_direction_X[i] = -999.0;
    sh_direction_Y[i] = -999.0;
    sh_direction_Z[i] = -999.0;
    sh_start_X[i] = -999.0;
    sh_start_Y[i] = -999.0;
    sh_start_Z[i] = -999.0;
    sh_length[i] = -999.0;
    sh_hasPrimary_e[i] = -999;
    sh_ehit_Q[i] = 0.0;
    sh_allhit_Q[i] = 0.0;
    sh_purity[i] = 0.0;
    sh_completeness[i] = 0.0;
  }

}

Member Data Documentation

int NuShowerEff::count_primary_e_in_Event

private

Definition at line 147 of file NuShowerEff_module.cc.

double NuShowerEff::esh_1_completeness

private

Definition at line 142 of file NuShowerEff_module.cc.

double NuShowerEff::esh_1_purity

private

Definition at line 141 of file NuShowerEff_module.cc.

double NuShowerEff::esh_completeness

private

Definition at line 146 of file NuShowerEff_module.cc.

double NuShowerEff::esh_each_completeness[MAX_SHOWERS]

private

Definition at line 144 of file NuShowerEff_module.cc.

double NuShowerEff::esh_each_purity[MAX_SHOWERS]

private

Definition at line 143 of file NuShowerEff_module.cc.

double NuShowerEff::esh_purity

private

Definition at line 145 of file NuShowerEff_module.cc.

int NuShowerEff::Event

private

Definition at line 101 of file NuShowerEff_module.cc.

TTree* NuShowerEff::fEventTree

private

Definition at line 150 of file NuShowerEff_module.cc.

float NuShowerEff::fFidVolCutX

private

Definition at line 88 of file NuShowerEff_module.cc.

float NuShowerEff::fFidVolCutY

private

Definition at line 89 of file NuShowerEff_module.cc.

float NuShowerEff::fFidVolCutZ

private

Definition at line 90 of file NuShowerEff_module.cc.

float NuShowerEff::fFidVolXmax

private

Definition at line 94 of file NuShowerEff_module.cc.

float NuShowerEff::fFidVolXmin

private

Definition at line 93 of file NuShowerEff_module.cc.

float NuShowerEff::fFidVolYmax

private

Definition at line 96 of file NuShowerEff_module.cc.

float NuShowerEff::fFidVolYmin

private

Definition at line 95 of file NuShowerEff_module.cc.

float NuShowerEff::fFidVolZmax

private

Definition at line 98 of file NuShowerEff_module.cc.

float NuShowerEff::fFidVolZmin

private

Definition at line 97 of file NuShowerEff_module.cc.

std::string NuShowerEff::fHitModuleLabel

private

Definition at line 77 of file NuShowerEff_module.cc.

bool NuShowerEff::fHitShowerThroughPFParticle

private

Definition at line 85 of file NuShowerEff_module.cc.

int NuShowerEff::fLeptonPDGcode

private

Definition at line 83 of file NuShowerEff_module.cc.

std::string NuShowerEff::fMCTruthModuleLabel

private

Definition at line 76 of file NuShowerEff_module.cc.

double NuShowerEff::fMinCompleteness

private

Definition at line 87 of file NuShowerEff_module.cc.

double NuShowerEff::fMinPurity

private

Definition at line 86 of file NuShowerEff_module.cc.

int NuShowerEff::fNeutrinoPDGcode

private

Definition at line 82 of file NuShowerEff_module.cc.

bool NuShowerEff::fSaveMCTree

private

Definition at line 84 of file NuShowerEff_module.cc.

std::string NuShowerEff::fShowerModuleLabel

private

Definition at line 78 of file NuShowerEff_module.cc.

std::string NuShowerEff::fTruthMatchDataModuleLabel

private

Definition at line 79 of file NuShowerEff_module.cc.

TH1D* NuShowerEff::h_Ee_den

private

Definition at line 155 of file NuShowerEff_module.cc.

TH1D* NuShowerEff::h_Ee_num

private

Definition at line 156 of file NuShowerEff_module.cc.

TEfficiency* NuShowerEff::h_Eff_Ee = 0

private

Definition at line 159 of file NuShowerEff_module.cc.

TEfficiency* NuShowerEff::h_Eff_Ev = 0

private

Definition at line 158 of file NuShowerEff_module.cc.

TH1D* NuShowerEff::h_Ev_den

private

Definition at line 152 of file NuShowerEff_module.cc.

TH1D* NuShowerEff::h_Ev_num

private

Definition at line 153 of file NuShowerEff_module.cc.

int NuShowerEff::MC_channel

private

Definition at line 109 of file NuShowerEff_module.cc.

double NuShowerEff::MC_incoming_E

private

Definition at line 123 of file NuShowerEff_module.cc.

double NuShowerEff::MC_incoming_P[4]

private

Definition at line 114 of file NuShowerEff_module.cc.

int NuShowerEff::MC_incoming_PDG

private

Definition at line 106 of file NuShowerEff_module.cc.

int NuShowerEff::MC_isCC

private

Definition at line 108 of file NuShowerEff_module.cc.

double NuShowerEff::MC_lepton_endMomentum[40]

private

Definition at line 116 of file NuShowerEff_module.cc.

double NuShowerEff::MC_lepton_endXYZT[4]

private

Definition at line 118 of file NuShowerEff_module.cc.

int NuShowerEff::MC_lepton_PDG

private

Definition at line 107 of file NuShowerEff_module.cc.

double NuShowerEff::MC_lepton_startEnergy

private

Definition at line 124 of file NuShowerEff_module.cc.

double NuShowerEff::MC_lepton_startMomentum[4]

private

Definition at line 115 of file NuShowerEff_module.cc.

double NuShowerEff::MC_lepton_startXYZT[4]

private

Definition at line 117 of file NuShowerEff_module.cc.

double NuShowerEff::MC_lepton_theta

private

Definition at line 119 of file NuShowerEff_module.cc.

int NuShowerEff::MC_leptonID

private

Definition at line 120 of file NuShowerEff_module.cc.

int NuShowerEff::MC_LeptonTrack

private

Definition at line 121 of file NuShowerEff_module.cc.

double NuShowerEff::MC_Q2

private

Definition at line 111 of file NuShowerEff_module.cc.

int NuShowerEff::MC_target

private

Definition at line 110 of file NuShowerEff_module.cc.

double NuShowerEff::MC_vertex[4]

private

Definition at line 113 of file NuShowerEff_module.cc.

double NuShowerEff::MC_W

private

Definition at line 112 of file NuShowerEff_module.cc.

int NuShowerEff::n_recoShowers

private

Definition at line 127 of file NuShowerEff_module.cc.

int NuShowerEff::Run

private

Definition at line 102 of file NuShowerEff_module.cc.

double NuShowerEff::sh_allhit_Q[MAX_SHOWERS]

private

Definition at line 138 of file NuShowerEff_module.cc.

double NuShowerEff::sh_completeness[MAX_SHOWERS]

private

Definition at line 140 of file NuShowerEff_module.cc.

double NuShowerEff::sh_direction_X[MAX_SHOWERS]

private

Definition at line 128 of file NuShowerEff_module.cc.

double NuShowerEff::sh_direction_Y[MAX_SHOWERS]

private

Definition at line 129 of file NuShowerEff_module.cc.

double NuShowerEff::sh_direction_Z[MAX_SHOWERS]

private

Definition at line 130 of file NuShowerEff_module.cc.

double NuShowerEff::sh_ehit_Q[MAX_SHOWERS]

private

Definition at line 135 of file NuShowerEff_module.cc.

int NuShowerEff::sh_hasPrimary_e[MAX_SHOWERS]

private

Definition at line 137 of file NuShowerEff_module.cc.

double NuShowerEff::sh_length[MAX_SHOWERS]

private

Definition at line 134 of file NuShowerEff_module.cc.

double NuShowerEff::sh_purity[MAX_SHOWERS]

private

Definition at line 139 of file NuShowerEff_module.cc.

double NuShowerEff::sh_start_X[MAX_SHOWERS]

private

Definition at line 131 of file NuShowerEff_module.cc.

double NuShowerEff::sh_start_Y[MAX_SHOWERS]

private

Definition at line 132 of file NuShowerEff_module.cc.

double NuShowerEff::sh_start_Z[MAX_SHOWERS]

private

Definition at line 133 of file NuShowerEff_module.cc.

int NuShowerEff::sh_TrackId[MAX_SHOWERS]

private

Definition at line 136 of file NuShowerEff_module.cc.

int NuShowerEff::SubRun

private

Definition at line 103 of file NuShowerEff_module.cc.

---

The documentation for this class was generated from the following file:

• larreco/larreco/ShowerFinder/NuShowerEff_module.cc