Inheritance diagram for MichelReco::MichelReco:

Public Member Functions
	MichelReco (fhicl::ParameterSet const &p)

	MichelReco (MichelReco const &)=delete

	MichelReco (MichelReco &&)=delete

MichelReco &	operator= (MichelReco const &)=delete

MichelReco &	operator= (MichelReco &&)=delete

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

void	beginJob () override

void	endJob () override

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

void	reconfigure (fhicl::ParameterSet const &p)

bool	insideFidVol (double pos[3])

bool	hitCloseToTrackEnd (detinfo::DetectorPropertiesData const &detProp, double radius, double end[3], double end2D[2], recob::Hit hit, geo::GeometryCore const &geom)

int	trackMatching (int trackIndex, art::FindManyP< recob::Hit > hitsFromTracks)

bool	isMuonDecaying (simb::MCParticle particle, std::vector< simb::MCParticle > particles)

bool	areHitsMichel (detinfo::DetectorClocksData const &clockData, const std::vector< recob::Hit > &hits)

	MichelReco (fhicl::ParameterSet const &p)

	MichelReco (MichelReco const &)=delete

	MichelReco (MichelReco &&)=delete

MichelReco &	operator= (MichelReco const &)=delete

MichelReco &	operator= (MichelReco &&)=delete

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

void	beginJob () override

void	endJob () override

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

void	reconfigure (fhicl::ParameterSet const &p)

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)

Public Attributes
calo::CalorimetryAlg	fCalorimetryAlg

Private Attributes
size_t	fEvNumber

TCanvas *	c

double	fRadiusThreshold

int	fNumberThreshold

double	fCNNThreshold

int	fNumberCloseHitsStart [3]

int	fNumberCloseHitsEnd [3]

double	fYesSelected

double	fNoSelected

int	fNMichel

int	fBestview

art::InputTag	fTrackModuleLabel

art::InputTag	fNNetModuleLabel

art::InputTag	fParticleModuleLabel

double	fFidVolCut

int	fNMichelHits

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 63 of file MichelEventSelection_module.cc.

Constructor & Destructor Documentation

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

explicit

Definition at line 137 of file MichelEventSelection_module.cc.

   : 
   EDAnalyzer(p),
   fCalorimetryAlg(p.get<fhicl::ParameterSet>("CalorimetryAlg"))
 {
   reconfigure(p);
 }

MichelReco::MichelReco::MichelReco ( MichelReco const & )

delete

MichelReco::MichelReco::MichelReco ( MichelReco && )

delete

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

explicit

MichelReco::MichelReco::MichelReco ( MichelReco const & )

delete

MichelReco::MichelReco::MichelReco ( MichelReco && )

delete

Member Function Documentation

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

overridevirtual

Implements art::EDAnalyzer.

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

overridevirtual

Implements art::EDAnalyzer.

Definition at line 170 of file MichelEventSelection_module.cc.

 {
 
   fEvNumber = evt.id().event();
   std::cout << "Michel event selection is on event " << fEvNumber << std::endl;
 
   // find all decaying muons in truth
   auto const & particles = *evt.getValidHandle<std::vector<simb::MCParticle>>( fParticleModuleLabel );
   for ( auto const & particle : particles ) {
     if ( abs(particle.PdgCode()) == 13 ) {
       double mcEnd[3] = { particle.EndX(), particle.EndY(), particle.EndZ() };
       if ( isMuonDecaying( particle, particles ) && insideFidVol( mcEnd ) ) {
         fNMichel += 1;
       }
     }
   }
   
   // Get the reconstructed objects, tracks and hits
   auto const & trackHandle = evt.getValidHandle<std::vector<recob::Track>>( fTrackModuleLabel );
   auto const & tracks = * trackHandle;
   art::FindManyP<recob::Hit> hitsFromTracks( trackHandle, evt, fTrackModuleLabel );
 
   auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService>()->DataFor(evt);
   auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService>()->DataFor(evt, clockData);
 
   // loop over tracks
   for (auto const & track : tracks) {
     
     // find best true trask ID for this track
     //int bestTrackId = trackMatching( &track - &tracks[0], hitsFromTracks );
     
     // find track start and end
     auto const & trackStart = track.Vertex<TVector3>();
     double start[3] = { trackStart.X(), trackStart.Y(), trackStart.Z() };
     auto const & trackEnd = track.End<TVector3>();
     double end[3] = { trackEnd.X(), trackEnd.Y(), trackEnd.Z() };
 
     // Check if the track starts or ends in the fiducial volume
     bool startInFidVol = insideFidVol( start );
     bool endInFidVol = insideFidVol( end );
     if ( startInFidVol || endInFidVol ) {
 
       // reset hit counts
       for (int plane = 0; plane < 3; plane++) {
         fNumberCloseHitsStart[ plane ] = 0; fNumberCloseHitsEnd[ plane ] = 0;
       }
 
       // get the hit results helper
       anab::MVAReader<recob::Hit,4> hitResults( evt, fNNetModuleLabel );
 
       // loop over hits store the selected hits
       std::vector< recob::Hit > taggedHitsStart;
       std::vector< recob::Hit > taggedHitsEnd;
       for ( size_t h = 0; h < hitResults.size(); ++h ) {
         
         // cnn output vector: [em,trk,none,michel]
         std::array<float,4> cnn_out = hitResults.getOutput( h );
         
         // keep hits with large CNN values
         if ( cnn_out[hitResults.getIndex("michel")] > fCNNThreshold ) {
           
           // need to check at both ends of the track
           if ( startInFidVol ) {
             
             // In each plane check if the michel tagged hits are close to the end of the track
             for ( int plane = 0; plane < 3; plane++ ) {
               
               // project the track onto the plane
               geo::GeometryCore const & geom = *art::ServiceHandle<geo::Geometry>();
               auto const & trackStart2D = pma::GetProjectionToPlane( trackStart, plane, geom.FindTPCAtPosition( start ).TPC, geom.FindCryostatAtPosition( start ) );
               double start2D[2] = {trackStart2D.X(), trackStart2D.Y()};
 
               // get the hit itself
               const recob::Hit & hit = hitResults.item(h);
               if (hit.View() == plane) {
 
                 // check that the hit is close to the track endpoint
                 // and add to tagged hits if so
                 if ( hitCloseToTrackEnd( detProp, fRadiusThreshold, start, start2D, hit, geom ) ) {
                   fNumberCloseHitsStart[ plane ] += 1;
                   if ( plane == 2 ) { taggedHitsStart.push_back( hit ); };
                 }
 
               }
 
             }
 
           }
           if ( endInFidVol ) {
 
             // In each plane check if the michel tagged hits are close to the end of the track
             for ( int plane = 0; plane < 3; plane++ ) {
              
               // project the track onto the plane 
               geo::GeometryCore const & geom = *art::ServiceHandle<geo::Geometry>();
               auto const & trackEnd2D = pma::GetProjectionToPlane( trackEnd, plane, geom.FindTPCAtPosition( end ).TPC, geom.FindCryostatAtPosition( end ) );
               double end2D[2] = {trackEnd2D.X(), trackEnd2D.Y()};
 
               // get the hit itself
               const recob::Hit & hit = hitResults.item(h);
               if (hit.View() == plane) {
 
                 // check that the hit is close to the track endpoint
                 // and add to tagged hits if so
                 if ( hitCloseToTrackEnd( detProp, fRadiusThreshold, end, end2D, hit, geom ) ) {
                   fNumberCloseHitsEnd[ plane ] += 1;
                   if ( plane == 2 ) { taggedHitsEnd.push_back( hit ); }
                 }
 
               }
 
             }
 
           }
 
         }
 
       } // end of loop over hits
 
       // event selection decision: clusters of hits in all planes at one end of the track
       bool startSelected = ( fNumberCloseHitsStart[0] > fNumberThreshold && fNumberCloseHitsStart[1] > fNumberThreshold && fNumberCloseHitsStart[2] > fNumberThreshold );
       bool endSelected = ( fNumberCloseHitsEnd[0] > fNumberThreshold && fNumberCloseHitsEnd[1] > fNumberThreshold && fNumberCloseHitsEnd[2] > fNumberThreshold );
       bool eventSelected = ( startSelected || endSelected );
       if ( startSelected && endSelected ) { eventSelected = false; }
       
       // check if the event was correclty tagged
       //  i.e. do the tagged hits correspond to a michel 
       bool particleIsMichel(false);
       if (eventSelected) {
         if ( startSelected ) { particleIsMichel = areHitsMichel(clockData, taggedHitsStart); }
         if ( endSelected ) { particleIsMichel = areHitsMichel(clockData, taggedHitsEnd); }
       }
       
       // update purity and efficiecny numbers and draw example events
       if ( particleIsMichel && eventSelected ) {
         
         fYesSelected += 1; 
         
         // plot the region around track end to see what the event is
         // if ( startSelected ) {
         //   
         //   for (int plane = 0; plane < 3; plane++) { 
         //   
         //     geo::GeometryCore const & geom = *art::ServiceHandle<geo::Geometry>();
         //     auto const & trackStart2D = pma::GetProjectionToPlane(trackStart, plane, geom.FindTPCAtPosition(start).TPC, geom.FindCryostatAtPosition(start));
         //     double start2D[2] = {trackStart2D.X(), trackStart2D.Y()};
 
         //     TH2D * histCNN = new TH2D("cnn", "Correctly Tagged", 200, start2D[0] - 20.0, start2D[0] +20.0, 200, start2D[1] - 20.0, start2D[1] + 20.0); 
         //     TH2D * histCharge = new TH2D("charge", "Correctly Tagged", 200, start2D[0] - 20.0, start2D[0] + 20.0, 200, start2D[1] - 20.0, start2D[1] + 20.0); 
         //     histCNN -> SetStats(0);
         //     histCharge -> SetStats(0);
 
         //     for ( size_t h = 0; h< hitResults.size(); ++h ) {
         //       const recob::Hit & hit = hitResults.item(h);
         //       auto const & hitLocation = pma::WireDriftToCm(hit.WireID().Wire, hit.PeakTime(), plane, geom.FindTPCAtPosition(start).TPC, geom.FindCryostatAtPosition(start)); 
         //     
         //       if (hit.View() == plane) {
         //         histCNN -> Fill(hitLocation.X(), hitLocation.Y(), hitResults.getOutput(h)[hitResults.getIndex("michel")]);
         //         histCharge -> Fill(hitLocation.X(), hitLocation.Y(), hit.Integral());
         //       }
         //     } 
         //   
         //     histCNN -> Draw("colz");
         //     c->SaveAs(Form("imgs/correctCNN%d_plane%d.png",bestTrackId, plane));
         //     histCharge-> Draw("colz");
         //     c->SaveAs(Form("imgs/correctCharge%d_plane%d.png",bestTrackId, plane));
         // 
         //   }
         // }
         // if ( endSelected ) {
         //   
         //   for (int plane = 0; plane < 3; plane++) { 
         //   
         //     geo::GeometryCore const & geom = *art::ServiceHandle<geo::Geometry>();
         //     auto const & trackEnd2D = pma::GetProjectionToPlane(trackEnd, plane, geom.FindTPCAtPosition(end).TPC, geom.FindCryostatAtPosition(end));
         //     double end2D[2] = {trackEnd2D.X(), trackEnd2D.Y()};
 
         //     TH2D * histCNN = new TH2D("cnn", "Correctly Tagged", 200, end2D[0] - 20.0, end2D[0] +20.0, 200, end2D[1] - 20.0, end2D[1] + 20.0); 
         //     TH2D * histCharge = new TH2D("charge", "Correctly Tagged", 200, end2D[0] - 20.0, end2D[0] + 20.0, 200, end2D[1] - 20.0, end2D[1] + 20.0); 
         //     histCNN -> SetStats(0);
         //     histCharge -> SetStats(0);
 
         //     for ( size_t h = 0; h< hitResults.size(); ++h ) {
         //       const recob::Hit & hit = hitResults.item(h);
         //       auto const & hitLocation = pma::WireDriftToCm(hit.WireID().Wire, hit.PeakTime(), plane, geom.FindTPCAtPosition(end).TPC, geom.FindCryostatAtPosition(end)); 
         //     
         //       if (hit.View() == plane) {
         //         histCNN -> Fill(hitLocation.X(), hitLocation.Y(), hitResults.getOutput(h)[hitResults.getIndex("michel")]);
         //         histCharge -> Fill(hitLocation.X(), hitLocation.Y(), hit.Integral());
         //       }
         //     } 
         //   
         //     histCNN -> Draw("colz");
         //     c->SaveAs(Form("imgs/correctCNN%d_plane%d.png",bestTrackId, plane));
         //     histCharge-> Draw("colz");
         //     c->SaveAs(Form("imgs/correctCharge%d_plane%d.png",bestTrackId, plane));
         // 
         //   }
 
         // }
 
       }
       if ( !particleIsMichel && eventSelected ) {
         
         fNoSelected += 1; 
         
         // plot the region around track end to see what the event is
         // if ( startSelected ) {
         //   
         //   for (int plane = 0; plane < 3; plane++) { 
         //   
         //     geo::GeometryCore const & geom = *art::ServiceHandle<geo::Geometry>();
         //     auto const & trackStart2D = pma::GetProjectionToPlane(trackStart, plane, geom.FindTPCAtPosition(start).TPC, geom.FindCryostatAtPosition(start));
         //     double start2D[2] = {trackStart2D.X(), trackStart2D.Y()};
 
         //     TH2D * histCNN = new TH2D("cnn", "Incorrectly Tagged", 200, start2D[0] - 20.0, start2D[0] +20.0, 200, start2D[1] - 20.0, start2D[1] + 20.0); 
         //     TH2D * histCharge = new TH2D("charge", "Incorrectly Tagged", 200, start2D[0] - 20.0, start2D[0] + 20.0, 200, start2D[1] - 20.0, start2D[1] + 20.0); 
         //     histCNN -> SetStats(0);
         //     histCharge -> SetStats(0);
 
         //     for ( size_t h = 0; h< hitResults.size(); ++h ) {
         //       const recob::Hit & hit = hitResults.item(h);
         //       auto const & hitLocation = pma::WireDriftToCm(hit.WireID().Wire, hit.PeakTime(), plane, geom.FindTPCAtPosition(start).TPC, geom.FindCryostatAtPosition(start)); 
         //     
         //       if (hit.View() == plane) {
         //         histCNN -> Fill(hitLocation.X(), hitLocation.Y(), hitResults.getOutput(h)[hitResults.getIndex("michel")]);
         //         histCharge -> Fill(hitLocation.X(), hitLocation.Y(), hit.Integral());
         //       }
         //     } 
         //   
         //     histCNN -> Draw("colz");
         //     c->SaveAs(Form("imgs/incorrectCNN%d_plane%d.png",bestTrackId, plane));
         //     histCharge-> Draw("colz");
         //     c->SaveAs(Form("imgs/incorrectCharge%d_plane%d.png",bestTrackId, plane));
         // 
         //   }
         // }
         // if ( endSelected ) {
         //   
         //   for (int plane = 0; plane < 3; plane++) { 
         //   
         //     geo::GeometryCore const & geom = *art::ServiceHandle<geo::Geometry>();
         //     auto const & trackEnd2D = pma::GetProjectionToPlane(trackEnd, plane, geom.FindTPCAtPosition(end).TPC, geom.FindCryostatAtPosition(end));
         //     double end2D[2] = {trackEnd2D.X(), trackEnd2D.Y()};
 
         //     TH2D * histCNN = new TH2D("cnn", "Incorrectly Tagged", 200, end2D[0] - 20.0, end2D[0] +20.0, 200, end2D[1] - 20.0, end2D[1] + 20.0); 
         //     TH2D * histCharge = new TH2D("charge", "Incorrectly Tagged", 200, end2D[0] - 20.0, end2D[0] + 20.0, 200, end2D[1] - 20.0, end2D[1] + 20.0); 
         //     histCNN -> SetStats(0);
         //     histCharge -> SetStats(0);
 
         //     for ( size_t h = 0; h< hitResults.size(); ++h ) {
         //       const recob::Hit & hit = hitResults.item(h);
         //       auto const & hitLocation = pma::WireDriftToCm(hit.WireID().Wire, hit.PeakTime(), plane, geom.FindTPCAtPosition(end).TPC, geom.FindCryostatAtPosition(end)); 
         //     
         //       if (hit.View() == plane) {
         //         histCNN -> Fill(hitLocation.X(), hitLocation.Y(), hitResults.getOutput(h)[hitResults.getIndex("michel")]);
         //         histCharge -> Fill(hitLocation.X(), hitLocation.Y(), hit.Integral());
         //       }
         //     } 
         //   
         //     histCNN -> Draw("colz");
         //     c->SaveAs(Form("imgs/incorrectCNN%d_plane%d.png",bestTrackId, plane));
         //     histCharge-> Draw("colz");
         //     c->SaveAs(Form("imgs/incorrectCharge%d_plane%d.png",bestTrackId, plane));
         // 
         //   }
         // 
         // }
 
       }
 
     }
 
   } // end of loop over tracks
 
 }

bool MichelReco::MichelReco::areHitsMichel	(	detinfo::DetectorClocksData const &	clockData,
		const std::vector< recob::Hit > &	hits
	)

Definition at line 576 of file MichelEventSelection_module.cc.

                                                                     {
 
   const simb::MCParticle *  mcParticle = 0;
 
   art::ServiceHandle< cheat::BackTrackerService > bt_serv;
   art::ServiceHandle< cheat::ParticleInventoryService > pi_serv;
   std::unordered_map< int, double > trkIDE;
   for ( auto const & hit : hits ) {
     for ( auto const & ide : bt_serv->HitToTrackIDEs(clockData, hit) ) { trkIDE[ide.trackID] += ide.energy; }
   }
   
   int best_id(0);
   double tot_e(0), max_e(0);
   for ( auto const & contrib : trkIDE ) {
 
     tot_e += contrib.second;
     if ( contrib.second > max_e ) {
       max_e = contrib.second;
       best_id = contrib.first;
     }
   }
 
   if ( (max_e > 0) && (tot_e > 0) ) {
     if ( best_id < 0 ) {
       best_id = -best_id;
     }
     mcParticle = pi_serv->TrackIdToParticle_P( best_id );
   }
  
   if (mcParticle != 0) { 
     return ( abs(mcParticle->PdgCode()) && ( mcParticle->Process() == "Decay" || mcParticle->Process() == "muMinusCaptureAtRest") ) ;
   }
   else {
     std::cout << "No match found" << hits.size() << std::endl;
     return false; 
   }
 
 }

void MichelReco::MichelReco::beginJob ( )

overridevirtual

Reimplemented from art::EDAnalyzer.

void MichelReco::MichelReco::beginJob ( )

overridevirtual

Reimplemented from art::EDAnalyzer.

Definition at line 145 of file MichelEventSelection_module.cc.

 {
   art::ServiceHandle<art::TFileService> tfs;
 }

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

override

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

override

Definition at line 165 of file MichelEventSelection_module.cc.

 {
   art::ServiceHandle<sim::LArG4Parameters> larParameters;
 }

void MichelReco::MichelReco::endJob ( )

overridevirtual

Reimplemented from art::EDAnalyzer.

void MichelReco::MichelReco::endJob ( )

overridevirtual

Reimplemented from art::EDAnalyzer.

Definition at line 150 of file MichelEventSelection_module.cc.

 {
   mf::LogVerbatim("MichelReco") << "MichelReco finished job";
   
   std::ofstream michelRecoOut;
   michelRecoOut.open("test.txt");
   michelRecoOut << "nMichel Purity Efficiency" << std::endl;
 
   double purity =  fYesSelected / ( fYesSelected + fNoSelected );
   double efficiency = fYesSelected/ fNMichel;
 
   michelRecoOut << fNMichel << " " << purity << " " << efficiency << std::endl;
   michelRecoOut.close();
 }

bool MichelReco::MichelReco::hitCloseToTrackEnd	(	detinfo::DetectorPropertiesData const &	detProp,
		double	radius,
		double	end[3],
		double	end2D[2],
		recob::Hit	hit,
		geo::GeometryCore const &	geom
	)

Definition at line 556 of file MichelEventSelection_module.cc.

                                                                                                                                 {
 
   bool close = false;
  
   auto const & hitLocation = pma::WireDriftToCm(detProp,
                                                 hit.WireID().Wire, hit.PeakTime(), hit.View(), geom.FindTPCAtPosition(end).TPC, geom.FindCryostatAtPosition(end));
 
   double deltaXToHit = hitLocation.X() - end2D[0];
   double deltaYToHit = hitLocation.Y() - end2D[1];
 
   double displacementToHit = pow(pow(deltaXToHit,2)+pow(deltaYToHit,2),0.5);
 
   if (displacementToHit < radius) close = true; 
 
   return close;
 
 }

bool MichelReco::MichelReco::insideFidVol ( double pos[3] )

Definition at line 507 of file MichelEventSelection_module.cc.

                                              {
 
   geo::GeometryCore const& geom = *art::ServiceHandle<geo::Geometry>();
   bool inside = false;
 
   geo::TPCID idtpc = geom.FindTPCAtPosition(pos);
 
   if (geom.HasTPC(idtpc)) {
 
     const geo::TPCGeo& tpcgeo = geom.GetElement(idtpc);
     double minx = tpcgeo.MinX(); double maxx = tpcgeo.MaxX();
     double miny = tpcgeo.MinY(); double maxy = tpcgeo.MaxY();
     double minz = tpcgeo.MinZ(); double maxz = tpcgeo.MaxZ();
     
     for (size_t cryo = 0; cryo < geom.Ncryostats(); cryo++) {
       const geo::CryostatGeo& cryostat = geom.Cryostat(cryo);
       for (size_t tpc = 0; tpc < cryostat.NTPC(); tpc++) {
         const geo::TPCGeo& tpcg = cryostat.TPC(tpc);
         if (tpcg.MinX() < minx) minx = tpcg.MinX();
         if (tpcg.MaxX() > maxx) maxx = tpcg.MaxX();
         if (tpcg.MinY() < miny) miny = tpcg.MinY();
         if (tpcg.MaxY() > maxy) maxy = tpcg.MaxY();
         if (tpcg.MinZ() < minz) minz = tpcg.MinZ();
         if (tpcg.MaxZ() > maxz) maxz = tpcg.MaxZ();
       }
     }
 
     double dista = fabs(minx - pos[0]);
     double distb = fabs(pos[0] - maxx); 
     if ((pos[0] > minx) && (pos[0] < maxx) && 
         (dista > fFidVolCut) && (distb > fFidVolCut)) inside = true;
       
     dista = fabs(miny - pos[1]);
     distb = fabs(pos[1]-maxy);
     if (inside && (pos[1] > miny) && (pos[1] < maxy) &&
         (dista > fFidVolCut) && (distb > fFidVolCut)) inside = true;
     else inside = false;
     
     dista = fabs(minz - pos[2]);
     distb = fabs(pos[2] - maxz);
     if (inside && (pos[2] > minz) && (pos[2] < maxz) &&
         (dista > fFidVolCut) && (distb > fFidVolCut)) inside = true;
     else inside = false;
 
   }
   return inside;
 }

bool MichelReco::MichelReco::isMuonDecaying	(	simb::MCParticle	particle,
		std::vector< simb::MCParticle >	particles
	)

Definition at line 484 of file MichelEventSelection_module.cc.

                                                                                                { 
 
   bool hasElectron = false, hasNuMu = false, hasNuE = false;
 
   unsigned int nSec = particle.NumberDaughters();
   for (size_t d = 0; d < nSec; ++d)
   {
     for (auto const & daughter : particles)
     {
       if (daughter.Mother() == particle.TrackId()) 
       {
         int d_pdg = abs(daughter.PdgCode());
         if (d_pdg == 11) hasElectron = true;
         else if (d_pdg == 14) hasNuMu = true;
         else if (d_pdg == 12) hasNuE = true;
       }
     }
   }
   return (hasElectron && hasNuMu && hasNuE);
 
 }

MichelReco& MichelReco::MichelReco::operator= ( MichelReco const & )

delete

MichelReco& MichelReco::MichelReco::operator= ( MichelReco && )

delete

MichelReco& MichelReco::MichelReco::operator= ( MichelReco const & )

delete

MichelReco& MichelReco::MichelReco::operator= ( MichelReco && )

delete

void MichelReco::MichelReco::reconfigure ( fhicl::ParameterSet const & p )

Definition at line 616 of file MichelEventSelection_module.cc.

                                                          {
 
   c = new TCanvas( "canv", "canv" );
 
   fTrackModuleLabel = p.get<std::string>("TrackModuleLabel");
   fNNetModuleLabel = p.get<std::string>("NNetModuleLabel");
   fParticleModuleLabel = p.get<std::string>("ParticleModuleLabel");
   
   fBestview = p.get<int>("Bestview");
 
   fRadiusThreshold = p.get<double>("MichelCloseHitRadius");
   fNumberThreshold = p.get<int>("CloseHitsThreshold"); 
   fCNNThreshold = p.get<double>("CNNThreshSelect");
   fFidVolCut = p.get<double>("FidVolCut"); 
   
   for ( int plane = 0; plane < 3; plane++ ) { 
     fNumberCloseHitsStart[ plane ] = 0; fNumberCloseHitsEnd[ plane ] = 0;
   }
 
   fYesSelected = 0;
   fNoSelected = 0;
   fNMichel = 0;
 
   return;
 }