Inheritance diagram for lar::example::AnalysisExample:

Classes
struct	Config

Public Types
using	Parameters = art::EDAnalyzer::Table< Config >

Public Types inherited from art::EDAnalyzer
using	WorkerType = WorkerT< EDAnalyzer >

using	ModuleType = EDAnalyzer

Public Member Functions
	AnalysisExample (Parameters const &config)
	Constructor: configures the module (see the Config structure above) More...

virtual void	beginJob () override

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

virtual void	analyze (const art::Event &event) override

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 Attributes
art::InputTag	fSimulationProducerLabel

art::InputTag	fHitProducerLabel
	The name of the producer that created hits. More...

art::InputTag	fClusterProducerLabel

int	fSelectedPDG
	PDG code of particle we'll focus on. More...

double	fBinSize
	For dE/dx work: the value of dx. More...

TH1D *	fPDGCodeHist
	PDG code of all particles. More...

TH1D *	fMomentumHist
	momentum [GeV] of all selected particles More...

TH1D *	fTrackLengthHist
	true length [cm] of all selected particles More...

TTree *	fSimulationNtuple
	tuple with simulated data More...

TTree *	fReconstructionNtuple
	tuple with reconstructed data More...

geo::GeometryCore const *	fGeometryService
	pointer to Geometry provider More...

double	fElectronsToGeV
	conversion factor More...

int	fTriggerOffset
	(units of ticks) time of expected neutrino event More...

The variables that will go into both n-tuples.
int	fEvent
	number of the event being processed More...

int	fRun
	number of the run being processed More...

int	fSubRun

The variables that will go into the simulation n-tuple.
int	fSimPDG
	PDG ID of the particle being processed. More...

int	fSimTrackID
	GEANT ID of the particle being processed. More...

double	fStartXYZT [4]
	(x,y,z,t) of the true start of the particle More...

double	fEndXYZT [4]
	(x,y,z,t) of the true end of the particle More...

double	fStartPE [4]
	(Px,Py,Pz,E) at the true start of the particle More...

double	fEndPE [4]
	(Px,Py,Pz,E) at the true end of the particle More...

int	fSimNdEdxBins
	Number of dE/dx bins in a given track. More...

std::vector< double >	fSimdEdxBins

Variables used in the reconstruction n-tuple
int	fRecoPDG
	PDG ID of the particle being processed. More...

int	fRecoTrackID
	GEANT ID of the particle being processed. More...

int	fRecoNdEdxBins
	Number of dE/dx bins in a given track. More...

std::vector< double >	fRecodEdxBins

Additional Inherited Members
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

Example analyzer module.

See also: analysis module example overview

This class extracts information from the generated and reconstructed particles.

It produces histograms for the simulated particles in the input file:

PDG ID (flavor) of all particles
momentum of the primary particles selected to have a specific PDG ID
length of the selected particle trajectory

It also produces two ROOT trees.

The first ROOT tree contains information on the simulated particles, including "dEdx", a binned histogram of collected charge as function of track range.

The second ROOT tree contains information on the reconstructed hits.

Configuration parameters

SimulationLabel (string, default: "largeant"): tag of the input data product with the detector simulation information (typically an instance of the LArG4 module)
HitLabel (string, mandatory): tag of the input data product with reconstructed hits
ClusterLabel (string, mandatory): tag of the input data product with reconstructed clusters
PDGcode (integer, mandatory): particle type (PDG ID) to be selected; only primary particles of this type will be considered
BinSize (real, mandatory): [cm] used for the calculation

Definition at line 177 of file AnalysisExample_module.cc.

Member Typedef Documentation

using lar::example::AnalysisExample::Parameters = art::EDAnalyzer::Table<Config>

Definition at line 235 of file AnalysisExample_module.cc.

Constructor & Destructor Documentation

lar::example::AnalysisExample::AnalysisExample ( Parameters const & config )

explicit

Constructor: configures the module (see the Config structure above)

Definition at line 373 of file AnalysisExample_module.cc.

       : EDAnalyzer(config)
       , fSimulationProducerLabel(config().SimulationLabel())
       , fHitProducerLabel(config().HitLabel())
       , fClusterProducerLabel(config().ClusterLabel())
       , fSelectedPDG(config().PDGcode())
       , fBinSize(config().BinSize())
     {
       // Get a pointer to the geometry service provider.
       fGeometryService = lar::providerFrom<geo::Geometry>();
 
       auto const clock_data =
         art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();
       fTriggerOffset = trigger_offset(clock_data);
 
       // Since art 2.8, you can and should tell beforehand, here in the
       // constructor, all the data the module is going to read ("consumes") or
       // might read
       // ("may_consume"). Diligence here will in the future help the framework
       // execute modules in parallel, making sure the order is correct.
       consumes<std::vector<simb::MCParticle>>(fSimulationProducerLabel);
       consumes<std::vector<sim::SimChannel>>(fSimulationProducerLabel);
       consumes<art::Assns<simb::MCTruth, simb::MCParticle>>(fSimulationProducerLabel);
       consumes<std::vector<recob::Hit>>(fHitProducerLabel);
       consumes<std::vector<recob::Cluster>>(fClusterProducerLabel);
       consumes<art::Assns<recob::Cluster, recob::Hit>>(fHitProducerLabel);
     }

Member Function Documentation

void lar::example::AnalysisExample::analyze ( const art::Event & event )

overridevirtual

Definition at line 483 of file AnalysisExample_module.cc.

     {
       // Start by fetching some basic event information for our n-tuple.
       fEvent = event.id().event();
       fRun = event.run();
       fSubRun = event.subRun();
 
       // This is the standard method of reading multiple objects
       // associated with the same event; see
       // <https://cdcvs.fnal.gov/redmine/projects/larsoft/wiki/Using_art_in_LArSoft>
       // for more information.
       //
       // Define a "handle" to point to a vector of the objects.
       art::Handle<std::vector<simb::MCParticle>> particleHandle;
 
       // Then tell the event to fill the vector with all the objects of
       // that type produced by a particular producer.
       //
       // Note that if I don't test whether this is successful, and there
       // aren't any simb::MCParticle objects, then the first time we
       // access particleHandle, art will display a "ProductNotFound"
       // exception message and, depending on art settings, it may skip
       // all processing for the rest of this event (including any
       // subsequent analysis steps) or stop the execution.
       if (!event.getByLabel(fSimulationProducerLabel, particleHandle)) {
         // If we have no MCParticles at all in an event, then we're in
         // big trouble. Throw an exception to force this module to
         // stop. Try to include enough information for the user to
         // figure out what's going on. Note that when we throw a
         // cet::exception, the run and event number will automatically
         // be displayed.
         //
         // __LINE__ and __FILE__ are values computed by the compiler.
 
         throw cet::exception("AnalysisExample")
           << " No simb::MCParticle objects in this event - "
           << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
       }
 
       // Get all the simulated channels for the event. These channels
       // include the energy deposited for each simulated track.
       //
       // Here we use a different method to access objects:
       // art::ValidHandle. Using this method, if there aren't any
       // objects of the given type (sim::SimChannel in this case) in the
       // input file for this art::Event, art will throw a
       // ProductNotFound exception.
       //
       // The "auto" type means that the C++ compiler will determine the
       // appropriate type for us, based on the return type of
       // art::Event::getValidHandle<T>(). The "auto" keyword is a great
       // timesaver, especially with frameworks like LArSoft which often
       // have complicated data product structures.
 
       auto simChannelHandle =
         event.getValidHandle<std::vector<sim::SimChannel>>(fSimulationProducerLabel);
 
       //
       // Let's compute the variables for the simulation n-tuple first.
       //
 
       // The MCParticle objects are not necessarily in any particular
       // order. Since we may have to search the list of particles, let's
       // put them into a map, a sorted container that will make
       // searching fast and easy. To save both space and time, the map
       // will not contain a copy of the MCParticle, but a pointer to it.
       std::map<int, const simb::MCParticle*> particleMap;
 
       // This is a "range-based for loop" in the 2011 version of C++; do
       // a web search on "c++ range based for loop" for more
       // information. Here's how it breaks down:
 
       // - A range-based for loop operates on a container.
       //   "particleHandle" is not a container; it's a pointer to a
       //   container. If we want C++ to "see" a container, we have to
       //   dereference the pointer, like this: *particleHandle.
 
       // - The loop variable that is set to each object in the container
       //   is named "particle". As for the loop variable's type:
 
       //   - To save a little bit of typing, and to make the code easier
       //     to maintain, we're going to let the C++ compiler deduce the
       //     type of what's in the container (simb::MCParticle objects
       //     in this case), so we use "auto".
 
       //   - We do _not_ want to change the contents of the container,
       //     so we use the "const" keyword to make sure.
 
       //   - We don't need to copy each object from the container into
       //     the variable "particle". It's sufficient to refer to the
       //     object by its address. So we use the reference operator "&"
       //     to tell the compiler to just copy the address, not the
       //     entire object.
 
       // It sounds complicated, but the end result is that we loop over
       // the list of particles in the art::Event in the most efficient
       // way possible.
 
       for (auto const& particle : (*particleHandle)) {
         // For the methods you can call to get particle information,
         // see ${NUSIMDATA_INC}/nusimdata/SimulationBase/MCParticle.h.
         fSimTrackID = particle.TrackId();
 
         // Add the address of the MCParticle to the map, with the
         // track ID as the key.
         particleMap[fSimTrackID] = &particle;
 
         // Histogram the PDG code of every particle in the event.
         fSimPDG = particle.PdgCode();
         fPDGCodeHist->Fill(fSimPDG);
 
         // For this example, we want to fill the n-tuples and histograms
         // only with information from the primary particles in the
         // event, whose PDG codes match a value supplied in the .fcl file.
         if (particle.Process() != "primary" || fSimPDG != fSelectedPDG) continue;
 
         // A particle has a trajectory, consisting of a set of
         // 4-positions and 4-mommenta.
         const size_t numberTrajectoryPoints = particle.NumberTrajectoryPoints();
 
         // For trajectories, as for vectors and arrays, the first
         // point is #0, not #1.
         const int last = numberTrajectoryPoints - 1;
         const TLorentzVector& positionStart = particle.Position(0);
         const TLorentzVector& positionEnd = particle.Position(last);
         const TLorentzVector& momentumStart = particle.Momentum(0);
         const TLorentzVector& momentumEnd = particle.Momentum(last);
 
         // Make a histogram of the starting momentum.
         fMomentumHist->Fill(momentumStart.P());
 
         // Fill arrays with the 4-values. (Don't be fooled by
         // the name of the method; it just puts the numbers from
         // the 4-vector into the array.)
         positionStart.GetXYZT(fStartXYZT);
         positionEnd.GetXYZT(fEndXYZT);
         momentumStart.GetXYZT(fStartPE);
         momentumEnd.GetXYZT(fEndPE);
 
         // Use a polar-coordinate view of the 4-vectors to
         // get the track length.
         const double trackLength = (positionEnd - positionStart).Rho();
 
         // Let's print some debug information in the job output to see
         // that everything is fine. MF_LOG_DEBUG() is a messagefacility
         // macro that prints stuff when the message level is set to
         // standard_debug in the .fcl file.
         MF_LOG_DEBUG("AnalysisExample") << "Track length: " << trackLength << " cm";
 
         // Fill a histogram of the track length.
         fTrackLengthHist->Fill(trackLength);
 
         MF_LOG_DEBUG("AnalysisExample")
           << "track ID=" << fSimTrackID << " (PDG ID: " << fSimPDG << ") " << trackLength
           << " cm long, momentum " << momentumStart.P() << " GeV/c, has " << numberTrajectoryPoints
           << " trajectory points";
 
         // Determine the number of dE/dx bins for the n-tuple.
         fSimNdEdxBins = int(trackLength / fBinSize) + 1;
         // Initialize the vector of dE/dx bins to be empty.
         fSimdEdxBins.clear();
 
         // To look at the energy deposited by this particle's track,
         // we loop over the SimChannel objects in the event.
         for (auto const& channel : (*simChannelHandle)) {
           // Get the numeric ID associated with this channel. (The
           // channel number is a 32-bit unsigned int, which normally
           // requires a special data type. Let's use "auto" so we
           // don't have to remember "raw::ChannelID_t".
           auto const channelNumber = channel.Channel();
 
           // A little care: There is more than one plane that reacts
           // to charge in the TPC. We only want to include the
           // energy from the collection plane. Note:
           // geo::kCollection is defined in
           // ${LARCOREOBJ_INC}/larcoreobj/SimpleTypesAndConstants/geo_types.h
           if (fGeometryService->SignalType(channelNumber) != geo::kCollection) continue;
 
           // Each channel has a map inside it that connects a time
           // slice to energy deposits in the detector. We'll use
           // "auto", but it's worth noting that the full type of
           // this map is
           // std::map<unsigned short, std::vector<sim::IDE>>
           auto const& timeSlices = channel.TDCIDEMap();
 
           // For every time slice in this channel:
           for (auto const& timeSlice : timeSlices) {
             // Each entry in a map is a pair<first,second>. For
             // the timeSlices map, the 'first' is a time slice
             // number. The 'second' is a vector of IDE objects.
             auto const& energyDeposits = timeSlice.second;
 
             // Loop over the energy deposits. An "energy deposit"
             // object is something that knows how much
             // charge/energy was deposited in a small volume, by
             // which particle, and where. The type of
             // 'energyDeposit' will be sim::IDE, which is defined
             // in
             // ${LARDATAOBJ_INC}/lardataobj/Simulation/SimChannel.h.
             for (auto const& energyDeposit : energyDeposits) {
               // Check if the track that deposited the
               // energy matches the track of the particle.
               if (energyDeposit.trackID != fSimTrackID) continue;
               // Get the (x,y,z) of the energy deposit.
               TVector3 location(energyDeposit.x, energyDeposit.y, energyDeposit.z);
 
               // Distance from the start of the track.
               const double distance = (location - positionStart.Vect()).Mag();
 
               // Into which bin of the dE/dx array do we add the energy?
               const unsigned int bin = (unsigned int)(distance / fBinSize);
 
               // Is the dE/dx array big enough to include this bin?
               if (fSimdEdxBins.size() < bin + 1) {
                 //  Increase the array size to accomodate
                 //  the new bin, padding it with zeros.
                 fSimdEdxBins.resize(bin + 1, 0.);
               }
 
               // Add the energy deposit to that bin. (If you look at the
               // definition of sim::IDE, you'll see that there's another
               // way to get the energy. Are the two methods equivalent?
               // Compare the results and see!)
               fSimdEdxBins[bin] += energyDeposit.numElectrons * fElectronsToGeV;
 
             } // For each energy deposit
           }   // For each time slice
         }     // For each SimChannel
 
         // At this point we've filled in the values of all the
         // variables and arrays we want to write to the n-tuple
         // for this particle. The following command actually
         // writes those values.
         fSimulationNtuple->Fill();
 
       } // loop over all particles in the event.
 
       //
       // Reconstruction n-tuple
       //
 
       // All of the above is based on objects entirely produced by the
       // simulation. Let's try to do something based on reconstructed
       // objects. A Hit (see ${LARDATAOBJ_INC}/lardataobj/RecoBase/Hit.h)
       // is a 2D object in a plane.
 
       // This code duplicates much of the code in the previous big
       // simulation loop, and will produce the similar results. (It
       // won't be identical, because of shaping and other factors; not
       // every bit of charge in a channel ends up contributing to a
       // hit.) The point is to show different methods of accessing
       // information, not to produce profound physics results -- that
       // part is up to you!
 
       // For the rest of this method, I'm going to assume you've read
       // the comments in previous section; I won't repeat all the C++
       // coding tricks and whatnot.
 
       // Start by reading the Hits. We don't use art::ValidHandle here
       // because there might be no hits in the input; e.g., we ran the
       // simulation but did not run the reconstruction steps yet. If
       // there are no hits we may as well skip the rest of this module.
 
       art::Handle<std::vector<recob::Hit>> hitHandle;
       if (!event.getByLabel(fHitProducerLabel, hitHandle)) return;
 
       // Our goal is to accumulate the dE/dx of any particles associated
       // with the hits that match our criteria: primary particles with
       // the PDG code from the .fcl file. I don't know how many such
       // particles there will be in a given event. We'll use a map, with
       // track ID as the key, to hold the vectors of dE/dx information.
       std::map<int, std::vector<double>> dEdxMap;
 
       auto const clock_data =
         art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
 
       // For every Hit:
       for (auto const& hit : (*hitHandle)) {
         // The channel associated with this hit.
         auto hitChannelNumber = hit.Channel();
 
         // We have a hit. For this example let's just focus on the
         // hits in the collection plane.
         if (fGeometryService->SignalType(hitChannelNumber) != geo::kCollection) continue;
 
         MF_LOG_DEBUG("AnalysisExample") << "Hit in collection plane" << std::endl;
 
         // In a few lines we're going to look for possible energy
         // deposits that correspond to that hit. Determine a
         // reasonable range of times that might correspond to those
         // energy deposits.
 
         // In reconstruction, the channel waveforms are truncated. So
         // we have to adjust the Hit TDC ticks to match those of the
         // SimChannels, which were created during simulation.
 
         // Save a bit of typing, while still allowing for potential
         // changes in the definitions of types in
         // $LARDATAOBJ_DIR/source/lardataobj/Simulation/SimChannel.h
 
         typedef sim::SimChannel::StoredTDC_t TDC_t;
         TDC_t start_tdc = clock_data.TPCTick2TDC(hit.StartTick());
         TDC_t end_tdc = clock_data.TPCTick2TDC(hit.EndTick());
         TDC_t hitStart_tdc = clock_data.TPCTick2TDC(hit.PeakTime() - 3. * hit.SigmaPeakTime());
         TDC_t hitEnd_tdc = clock_data.TPCTick2TDC(hit.PeakTime() + 3. * hit.SigmaPeakTime());
 
         start_tdc = std::max(start_tdc, hitStart_tdc);
         end_tdc = std::min(end_tdc, hitEnd_tdc);
 
         // In the simulation section, we started with particles to find
         // channels with a matching track ID. Now we search in reverse:
         // search the SimChannels for matching channel number, then look
         // at the tracks inside the channel.
 
         for (auto const& channel : (*simChannelHandle)) {
           auto simChannelNumber = channel.Channel();
           if (simChannelNumber != hitChannelNumber) continue;
 
           MF_LOG_DEBUG("AnalysisExample")
             << "SimChannel number = " << simChannelNumber << std::endl;
 
           // The time slices in this channel.
           auto const& timeSlices = channel.TDCIDEMap();
 
           // We want to look over the range of time slices in this
           // channel that correspond to the range of hit times.
 
           // To do this, we're going to use some fast STL search
           // methods; STL algorithms are usually faster than the
           // ones you might write yourself. The price for this speed
           // is a bit of code complexity: in particular, we need to
           // a custom comparison function, TDCIDETimeCompare, to
           // define a "less-than" function for the searches.
 
           // For an introduction to STL algorithms, see
           // <http://www.learncpp.com/cpp-tutorial/16-4-stl-algorithms-overview/>.
           // For a list of available STL algorithms, see
           // <http://en.cppreference.com/w/cpp/algorithm>
 
           // We have to create "dummy" time slices for the search.
           sim::TDCIDE startTime;
           sim::TDCIDE endTime;
           startTime.first = start_tdc;
           endTime.first = end_tdc;
 
           // Here are the fast searches:
 
           // Find a pointer to the first channel with time >= start_tdc.
           auto const startPointer =
             std::lower_bound(timeSlices.begin(), timeSlices.end(), startTime, TDCIDETimeCompare);
 
           // From that time slice, find the last channel with time < end_tdc.
           auto const endPointer =
             std::upper_bound(startPointer, timeSlices.end(), endTime, TDCIDETimeCompare);
 
           // Did we find anything? If not, skip.
           if (startPointer == timeSlices.end() || startPointer == endPointer) continue;
           MF_LOG_DEBUG("AnalysisExample")
             << "Time slice start = " << (*startPointer).first << std::endl;
 
           // Loop over the channel times we found that match the hit
           // times.
           for (auto slicePointer = startPointer; slicePointer != endPointer; slicePointer++) {
             auto const timeSlice = *slicePointer;
             auto time = timeSlice.first;
 
             // How to debug a problem: Lots of print statements. There are
             // debuggers such as gdb, but they can be tricky to use with
             // shared libraries and don't work if you're using software
             // that was compiled somewhere else (e.g., you're accessing
             // LArSoft libraries via CVMFS).
 
             // The MF_LOG_DEBUG statements below are left over from when I
             // was trying to solve a problem about hit timing. I could
             // have deleted them, but decided to leave them to demonsrate
             // what a typical debugging process looks like.
 
             MF_LOG_DEBUG("AnalysisExample")
               << "Hit index = " << hit.LocalIndex() << " channel number = " << hitChannelNumber
               << " start TDC tick = " << hit.StartTick() << " end TDC tick = " << hit.EndTick()
               << " peak TDC tick = " << hit.PeakTime()
               << " sigma peak time = " << hit.SigmaPeakTime()
               << " adjusted start TDC tick = " << clock_data.TPCTick2TDC(hit.StartTick())
               << " adjusted end TDC tick = " << clock_data.TPCTick2TDC(hit.EndTick())
               << " adjusted peak TDC tick = " << clock_data.TPCTick2TDC(hit.PeakTime())
               << " adjusted start_tdc = " << start_tdc << " adjusted end_tdc = " << end_tdc
               << " time = " << time << std::endl;
 
             // Loop over the energy deposits.
             auto const& energyDeposits = timeSlice.second;
             for (auto const& energyDeposit : energyDeposits) {
               // Remember that map of MCParticles we created
               // near the top of this method? Now we can use
               // it. Search the map for the track ID associated
               // with this energy deposit. Since a map is
               // automatically sorted, we can use a fast binary
               // search method, 'find()'.
 
               // By the way, the type of "search" is an iterator
               // (to be specific, it's an
               // std::map<int,simb::MCParticle*>::const_iterator,
               // which makes you thankful for the "auto"
               // keyword). If you're going to work with C++
               // containers, you'll have to learn about
               // iterators eventually; do a web search on "STL
               // iterator" to get started.
               auto search = particleMap.find(energyDeposit.trackID);
 
               // Did we find this track ID in the particle map?
               // It's possible for the answer to be "no"; some
               // particles are too low in kinetic energy to be
               // written by the simulation (see
               // ${LARSIM_DIR}/job/simulationservices.fcl,
               // parameter ParticleKineticEnergyCut).
               if (search == particleMap.end()) continue;
 
               // "search" points to a pair in the map: <track ID, MCParticle*>
               int trackID = (*search).first;
               const simb::MCParticle& particle = *((*search).second);
 
               // Is this a primary particle, with a PDG code that
               // matches the user input?
               if (particle.Process() != "primary" || particle.PdgCode() != fSelectedPDG) continue;
 
               // Determine the dE/dx of this particle.
               const TLorentzVector& positionStart = particle.Position(0);
               TVector3 location(energyDeposit.x, energyDeposit.y, energyDeposit.z);
               double distance = (location - positionStart.Vect()).Mag();
               unsigned int bin = int(distance / fBinSize);
               double energy = energyDeposit.numElectrons * fElectronsToGeV;
 
               // A feature of maps: if we refer to
               // dEdxMap[trackID], and there's no such entry in
               // the map yet, it will be automatically created
               // with a zero-size vector. Test to see if the
               // vector for this track ID is big enough.
               //
               // dEdxMap is a map, which is a slow container
               // compared to a vector. If we are going to access
               // the same element over and over, it is a good
               // idea to find that element once, and then refer
               // to that item directly. Since we don't care
               // about the specific type of dEdxMap[trackID] (a
               // vector, by the way), we can use "auto" to save
               // some time. This must be a reference, since we
               // want to change the original value in the map,
               // and can't be constant.
               auto& track_dEdX = dEdxMap[trackID];
               if (track_dEdX.size() < bin + 1) {
                 // Increase the vector size, padding it with
                 // zeroes.
                 track_dEdX.resize(bin + 1, 0);
               }
 
               // Add the energy to the dE/dx bins for this track.
               track_dEdX[bin] += energy;
 
             } // loop over energy deposits
           }   // loop over time slices
         }     // for each SimChannel
       }       // for each Hit
 
       // We have a map of dE/dx vectors. Write each one of them to the
       // reconstruction n-tuple.
       for (const auto& dEdxEntry : dEdxMap) {
         // Here, the map entries are <first,second>=<track ID, dE/dx vector>
         fRecoTrackID = dEdxEntry.first;
 
         // This is an example of how we'd pick out the PDG code if
         // there are multiple particle types or tracks in a single
         // event allowed in the n-tuple.
         fRecoPDG = particleMap[fRecoTrackID]->PdgCode();
 
         // Get the number of bins for this track.
         const std::vector<double>& dEdx = dEdxEntry.second;
         fRecoNdEdxBins = dEdx.size();
 
         // Copy this track's dE/dx information.
         fRecodEdxBins = dEdx;
 
         // At this point, we've filled in all the reconstruction
         // n-tuple's variables. Write it.
         fReconstructionNtuple->Fill();
       }
 
       // Think about the two big loops above, One starts from the
       // particles then looks at the channels; the other starts with the
       // hits and backtracks to the particles. What links the
       // information in simb::MCParticle and sim::SimChannel is the
       // track ID number assigned by the LArG4 simulation; what links
       // sim::SimChannel and recob::Hit is the channel ID.
 
       // In general, that is not how objects in the LArSoft
       // reconstruction chain are linked together. Most of them are
       // linked using associations and the art::Assns class:
       // <https://cdcvs.fnal.gov/redmine/projects/larsoft/wiki/Using_art_in_LArSoft#artAssns>
 
       // The web page may a bit difficult to understand (at least, it is
       // for me), so let's try to simplify things:
 
       // - If you're doing an analysis task, you don't have to worry
       //   about creating art::Assns objects.
 
       // - You don't have to read the art:Assns objects on your
       //   own. There are helper classes (FindXXX) which will do that for
       //   you.
 
       // - There's only one helper you need: art::FindManyP. It will do
       //   what you want with a minimum of fuss.
 
       // - Associations are symmetric. If you see an
       //   art::Assns<ClassA,ClassB>, the helper classes will find all of
       //   ClassA associated with ClassB or all of ClassB associated with
       //   ClassA.
 
       // - To know what associations exist, you have to be a 'code
       //   detective'. The important clue is to look for a 'produces' line
       //   in the code that writes objects to an art::Event. For example,
       //   in ${LARSIM_DIR}/source/larsim/LArG4/LArG4_module.cc, you'll see this
       //   line:
 
       //   produces< art::Assns<simb::MCTruth, simb::MCParticle> >();
 
       //   That means a simulated event will have an association between
       //   simb::MCTruth (the primary particle produced by the event
       //   generator) and the simb::MCParticle objects (the secondary
       //   particles produced in the detector simulation).
 
       // Let's try it. The following statement will find the
       // simb::MCTruth objects associated with the simb::MCParticle
       // objects in the event (referenced by particleHandle):
 
       const art::FindManyP<simb::MCTruth> findManyTruth(
         particleHandle, event, fSimulationProducerLabel);
 
       // Note that we still have to supply the module label of the step
       // that created the association. Also note that we did not have to
       // explicitly read in the simb::MCTruth objects from the
       // art::Event object 'event'; FindManyP did that for us.
 
       // Also note that at this point art::FindManyP has already found
       // all the simb::MCTruth associated with each of the particles in
       // particleHandle. This is a slow process, so in general you want
       // to do it only once. If we had a loop over the particles, we
       // would still do this outside that loop.
 
       // Now we can query the 'findManyTruth' object to access the
       // information. First, check that there wasn't some kind of error:
 
       if (!findManyTruth.isValid()) {
         mf::LogError("AnalysisExample")
           << "findManyTruth simb::MCTruth for simb::MCParticle failed!";
       }
 
       // I'm going to be lazy, and just look at the simb::MCTruth object
       // associated with the first simb::MCParticle we read in. (The
       // main reason I'm being lazy is that if I used the
       // single-particle generator in prodsingle.fcl, every particle in
       // the event is going to be associated with just the one primary
       // particle from the event generator.)
 
       size_t particle_index = 0; // look at first particle in
                                  // particleHandle's vector.
 
       // I'm using "auto" to save on typing. The result of
       // FindManyP::at() is a vector of pointers, in this case
       // simb::MCTruth*. In this case it will be a vector with just one
       // entry; I could have used art::FindOneP instead. (This will be a
       // vector of art::Ptr, which is a type of smart pointer; see
       // <https://cdcvs.fnal.gov/redmine/projects/larsoft/wiki/Using_art_in_LArSoft#artPtrltTgt-and-artPtrVectorltTgt>
       // To avoid unnecessary copying, and since art::FindManyP returns
       // a constant reference, use "auto const&".
 
       auto const& truth = findManyTruth.at(particle_index);
 
       // Make sure there's no problem.
       if (truth.empty()) {
         mf::LogError("AnalysisExample")
           << "Particle ID=" << particleHandle->at(particle_index).TrackId() << " has no primary!";
       }
 
       // Use the message facility to write output. I don't have to write
       // the event, run, or subrun numbers; the message facility takes
       // care of that automatically. I'm "going at warp speed" with the
       // vectors, pointers, and methods; see
       // ${NUSIMDATA_INC}/nusimdata/SimulationBase/MCTruth.h and
       // ${NUSIMDATA_INC}/nusimdata/SimulationBase/MCParticle.h for the
       // nested calls I'm using.
       mf::LogInfo("AnalysisExample")
         << "Particle ID=" << particleHandle->at(particle_index).TrackId()
         << " primary PDG code=" << truth[0]->GetParticle(0).PdgCode();
 
       // Let's try a slightly more realistic example. Suppose I want to
       // read in the clusters, and learn what hits are associated with
       // them. Then I could backtrack from those hits to determine the
       // dE/dx of the particles in the clusters. (Don't worry; I won't
       // put you through the n-tuple creation procedure for a third
       // time.)
 
       // First, read in the clusters (if there are any).
       art::Handle<std::vector<recob::Cluster>> clusterHandle;
       if (!event.getByLabel(fClusterProducerLabel, clusterHandle)) return;
 
       // Now use the associations to find which hits are associated with
       // which clusters. Note that this is not as trivial a query as the
       // one with MCTruth, since it's possible for a hit to be assigned
       // to more than one cluster.
 
       // We have to include fClusterProducerLabel, since that's the step
       // that created the art::Assns<recob::Hit,recob::Cluster> object;
       // look at the modules in ${LARRECO_DIR}/source/larreco/ClusterFinder/
       // and search for the 'produces' lines. (I did not know this before I
       // wrote these comments. I had to be a code detective and use UNIX
       // tools like 'grep' and 'find' to locate those routines.)
       const art::FindManyP<recob::Hit> findManyHits(clusterHandle, event, fClusterProducerLabel);
 
       if (!findManyHits.isValid()) {
         mf::LogError("AnalysisExample") << "findManyHits recob::Hit for recob::Cluster failed;"
                                         << " cluster label='" << fClusterProducerLabel << "'";
       }
 
       // Now we'll loop over the clusters to see the hits associated
       // with each one. Note that we're not using a range-based for
       // loop. That's because FindManyP::at() expects a number as an
       // argument, so we might as well go through the cluster objects
       // via numerical index instead.
       for (size_t cluster_index = 0; cluster_index != clusterHandle->size(); cluster_index++) {
         // In this case, FindManyP::at() will return a vector of
         // pointers to recob::Hit that corresponds to the
         // "cluster_index"-th cluster.
         auto const& hits = findManyHits.at(cluster_index);
 
         // We have a vector of pointers to the hits associated
         // with the cluster, but for this example I'm not going to
         // do anything fancy with them. I'll just print out how
         // many there are.
         mf::LogInfo("AnalysisExample") << "Cluster ID=" << clusterHandle->at(cluster_index).ID()
                                        << " has " << hits.size() << " hits";
       }
 
     } // AnalysisExample::analyze()

void lar::example::AnalysisExample::beginJob ( )

overridevirtual

Reimplemented from art::EDAnalyzer.

Definition at line 403 of file AnalysisExample_module.cc.

     {
       // Get the detector length, to determine the maximum bin edge of one
       // of the histograms.
       const double detectorLength = DetectorDiagonal(*fGeometryService);
 
       // Access ART's TFileService, which will handle creating and writing
       // histograms and n-tuples for us.
       art::ServiceHandle<art::TFileService const> tfs;
 
       // For TFileService, the arguments to 'make<whatever>' are the
       // same as those passed to the 'whatever' constructor, provided
       // 'whatever' is a ROOT class that TFileService recognizes.
 
       // Define the histograms. Putting semi-colons around the title
       // causes it to be displayed as the x-axis label if the histogram
       // is drawn (the format is "title;label on abscissae;label on ordinates").
       fPDGCodeHist = tfs->make<TH1D>("pdgcodes", ";PDG Code;", 5000, -2500, 2500);
       fMomentumHist = tfs->make<TH1D>("mom", ";particle Momentum (GeV);", 100, 0., 10.);
       fTrackLengthHist =
         tfs->make<TH1D>("length", ";particle track length (cm);", 200, 0, detectorLength);
 
       // Define our n-tuples, which are limited forms of ROOT
       // TTrees. Start with the TTree itself.
       fSimulationNtuple =
         tfs->make<TTree>("AnalysisExampleSimulation", "AnalysisExampleSimulation");
       fReconstructionNtuple =
         tfs->make<TTree>("AnalysisExampleReconstruction", "AnalysisExampleReconstruction");
 
       // Define the branches (columns) of our simulation n-tuple. To
       // write a variable, we give the address of the variable to
       // TTree::Branch.
       fSimulationNtuple->Branch("Event", &fEvent, "Event/I");
       fSimulationNtuple->Branch("SubRun", &fSubRun, "SubRun/I");
       fSimulationNtuple->Branch("Run", &fRun, "Run/I");
       fSimulationNtuple->Branch("TrackID", &fSimTrackID, "TrackID/I");
       fSimulationNtuple->Branch("PDG", &fSimPDG, "PDG/I");
       // When we write arrays, we give the address of the array to
       // TTree::Branch; in C++ this is simply the array name.
       fSimulationNtuple->Branch("StartXYZT", fStartXYZT, "StartXYZT[4]/D");
       fSimulationNtuple->Branch("EndXYZT", fEndXYZT, "EndXYZT[4]/D");
       fSimulationNtuple->Branch("StartPE", fStartPE, "StartPE[4]/D");
       fSimulationNtuple->Branch("EndPE", fEndPE, "EndPE[4]/D");
       // For a variable-length array: include the number of bins.
       fSimulationNtuple->Branch("NdEdx", &fSimNdEdxBins, "NdEdx/I");
       // ROOT branches can contain std::vector objects.
       fSimulationNtuple->Branch("dEdx", &fSimdEdxBins);
 
       // A similar definition for the reconstruction n-tuple. Note that we
       // use some of the same variables in both n-tuples.
       fReconstructionNtuple->Branch("Event", &fEvent, "Event/I");
       fReconstructionNtuple->Branch("SubRun", &fSubRun, "SubRun/I");
       fReconstructionNtuple->Branch("Run", &fRun, "Run/I");
       fReconstructionNtuple->Branch("TrackID", &fRecoTrackID, "TrackID/I");
       fReconstructionNtuple->Branch("PDG", &fRecoPDG, "PDG/I");
       fReconstructionNtuple->Branch("NdEdx", &fRecoNdEdxBins, "NdEdx/I");
       fReconstructionNtuple->Branch("dEdx", &fRecodEdxBins);
     }

void lar::example::AnalysisExample::beginRun ( const art::Run & run )

overridevirtual

Definition at line 471 of file AnalysisExample_module.cc.

     {
       // How to convert from number of electrons to GeV. The ultimate
       // source of this conversion factor is
       // ${LARCOREOBJ_INC}/larcoreobj/SimpleTypesAndConstants/PhysicalConstants.h.
       // But sim::LArG4Parameters might in principle ask a database for it.
       art::ServiceHandle<sim::LArG4Parameters const> larParameters;
       fElectronsToGeV = 1. / larParameters->GeVToElectrons();
     }

Member Data Documentation

double lar::example::AnalysisExample::fBinSize

private

For dE/dx work: the value of dx.

Definition at line 296 of file AnalysisExample_module.cc.

art::InputTag lar::example::AnalysisExample::fClusterProducerLabel

private

The name of the producer that created clusters

Definition at line 293 of file AnalysisExample_module.cc.

double lar::example::AnalysisExample::fElectronsToGeV

private

conversion factor

Definition at line 352 of file AnalysisExample_module.cc.

double lar::example::AnalysisExample::fEndPE[4]

private

(Px,Py,Pz,E) at the true end of the particle

Definition at line 326 of file AnalysisExample_module.cc.

double lar::example::AnalysisExample::fEndXYZT[4]

private

(x,y,z,t) of the true end of the particle

Definition at line 324 of file AnalysisExample_module.cc.

int lar::example::AnalysisExample::fEvent

private

number of the event being processed

Definition at line 310 of file AnalysisExample_module.cc.

geo::GeometryCore const* lar::example::AnalysisExample::fGeometryService

private

pointer to Geometry provider

Definition at line 351 of file AnalysisExample_module.cc.

art::InputTag lar::example::AnalysisExample::fHitProducerLabel

private

The name of the producer that created hits.

Definition at line 292 of file AnalysisExample_module.cc.

TH1D* lar::example::AnalysisExample::fMomentumHist

private

momentum [GeV] of all selected particles

Definition at line 300 of file AnalysisExample_module.cc.

TH1D* lar::example::AnalysisExample::fPDGCodeHist

private

PDG code of all particles.

Definition at line 299 of file AnalysisExample_module.cc.

std::vector<double> lar::example::AnalysisExample::fRecodEdxBins

private

The vector that will be used to accumulate dE/dx values as a function of range.

Definition at line 346 of file AnalysisExample_module.cc.

int lar::example::AnalysisExample::fRecoNdEdxBins

private

Number of dE/dx bins in a given track.

Definition at line 342 of file AnalysisExample_module.cc.

TTree* lar::example::AnalysisExample::fReconstructionNtuple

private

tuple with reconstructed data

Definition at line 305 of file AnalysisExample_module.cc.

int lar::example::AnalysisExample::fRecoPDG

private

PDG ID of the particle being processed.

Definition at line 338 of file AnalysisExample_module.cc.

int lar::example::AnalysisExample::fRecoTrackID

private

GEANT ID of the particle being processed.

Definition at line 339 of file AnalysisExample_module.cc.

int lar::example::AnalysisExample::fRun

private

number of the run being processed

Definition at line 311 of file AnalysisExample_module.cc.

int lar::example::AnalysisExample::fSelectedPDG

private

PDG code of particle we'll focus on.

Definition at line 295 of file AnalysisExample_module.cc.

std::vector<double> lar::example::AnalysisExample::fSimdEdxBins

private

The vector that will be used to accumulate dE/dx values as a function of range.

Definition at line 333 of file AnalysisExample_module.cc.

int lar::example::AnalysisExample::fSimNdEdxBins

private

Number of dE/dx bins in a given track.

Definition at line 329 of file AnalysisExample_module.cc.

int lar::example::AnalysisExample::fSimPDG

private

PDG ID of the particle being processed.

Definition at line 317 of file AnalysisExample_module.cc.

int lar::example::AnalysisExample::fSimTrackID

private

GEANT ID of the particle being processed.

Definition at line 318 of file AnalysisExample_module.cc.

TTree* lar::example::AnalysisExample::fSimulationNtuple

private

tuple with simulated data

Definition at line 304 of file AnalysisExample_module.cc.

art::InputTag lar::example::AnalysisExample::fSimulationProducerLabel

private

The name of the producer that tracked simulated particles through the detector

Definition at line 290 of file AnalysisExample_module.cc.

double lar::example::AnalysisExample::fStartPE[4]

private

(Px,Py,Pz,E) at the true start of the particle

Definition at line 325 of file AnalysisExample_module.cc.

double lar::example::AnalysisExample::fStartXYZT[4]

private

(x,y,z,t) of the true start of the particle

Definition at line 323 of file AnalysisExample_module.cc.

int lar::example::AnalysisExample::fSubRun

private

number of the sub-run being processed

Definition at line 312 of file AnalysisExample_module.cc.

TH1D* lar::example::AnalysisExample::fTrackLengthHist

private

true length [cm] of all selected particles

Definition at line 301 of file AnalysisExample_module.cc.

int lar::example::AnalysisExample::fTriggerOffset

private

(units of ticks) time of expected neutrino event

Definition at line 353 of file AnalysisExample_module.cc.

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

larexamples/larexamples/AnalysisExample/AnalysisExample_module.cc

Classes

Public Types

Public Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Configuration parameters

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation