Creates a simple ROOT tree with tracking and calorimetry information. More...

Inheritance diagram for dune::AnalysisTree:

Public Member Functions
	AnalysisTree (fhicl::ParameterSet const &pset)

virtual	~AnalysisTree ()

void	analyze (const art::Event &evt)
	read access to event More...

void	beginSubRun (const art::SubRun &sr)

void	endSubRun (const art::SubRun &sr)

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	HitsPurity (detinfo::DetectorClocksData const &clockData, std::vector< art::Ptr< recob::Hit > > const &hits, Int_t &trackid, Float_t &purity, double &maxe)

double	length (const recob::Track &track)

double	driftedLength (detinfo::DetectorPropertiesData const &detProp, const simb::MCParticle &part, TLorentzVector &start, TLorentzVector &end, unsigned int &starti, unsigned int &endi)

double	driftedLength (detinfo::DetectorPropertiesData const &detProp, const sim::MCTrack &mctrack, TLorentzVector &tpcstart, TLorentzVector &tpcend, TLorentzVector &tpcmom)

double	length (const simb::MCParticle &part, TLorentzVector &start, TLorentzVector &end, unsigned int &starti, unsigned int &endi)

double	bdist (const TVector3 &pos)

size_t	GetNTrackers () const
	Returns the number of trackers configured. More...

size_t	GetNVertexAlgos () const

size_t	GetNShowerAlgos () const
	Returns the number of shower algorithms configured. More...

std::vector< std::string >	GetShowerAlgos () const
	Returns the name of configured shower algorithms (converted to string) More...

void	CreateData (bool bClearData=false)
	Creates the structure for the tree data; optionally initializes it. More...

void	SetAddresses ()
	Sets the addresses of all the tree branches, creating the missing ones. More...

void	SetTrackerAddresses (size_t iTracker)

void	SetVertexAddresses (size_t iVertexAlg)

void	SetShowerAddresses (size_t iShower)

void	SetPFParticleAddress ()

void	CreateTree (bool bClearData=false)
	Create the output tree and the data structures, if needed. More...

void	DestroyData ()
	Destroy the local buffers (existing branches will point to invalid address!) More...

void	CheckData (std::string caller) const
	Helper function: throws if no data structure is available. More...

void	CheckTree (std::string caller) const
	Helper function: throws if no tree is available. More...

void	FillShower (AnalysisTreeDataStruct::ShowerDataStruct &showerData, size_t iShower, recob::Shower const &showers, const bool fSavePFParticleInfo, const std::map< Short_t, Short_t > &showerIDtoPFParticleIDMap) const
	Stores the information of shower in slot iShower of showerData. More...

void	FillShowers (AnalysisTreeDataStruct::ShowerDataStruct &showerData, std::vector< recob::Shower > const &showers, const bool fSavePFParticleInfo, const std::map< Short_t, Short_t > &showerIDtoPFParticleIDMap) const
	Stores the information of all showers into showerData. More...

Private Attributes
TTree *	fTree

TTree *	fPOT

std::unique_ptr< AnalysisTreeDataStruct >	fData

AnalysisTreeDataStruct::SubRunData_t	SubRunData

std::string	fDigitModuleLabel

std::string	fHitsModuleLabel

std::string	fLArG4ModuleLabel

std::string	fCalDataModuleLabel

std::string	fGenieGenModuleLabel

std::string	fCryGenModuleLabel

std::string	fProtoGenModuleLabel

std::string	fG4ModuleLabel

std::string	fClusterModuleLabel

std::string	fPandoraNuVertexModuleLabel

std::string	fOpFlashModuleLabel

std::string	fExternalCounterModuleLabel

std::string	fMCShowerModuleLabel

std::string	fMCTrackModuleLabel

std::string	fSpacePointSolverModuleLabel

std::string	fCnnModuleLabel

std::vector< std::string >	fTrackModuleLabel

std::string	fPFParticleModuleLabel

std::vector< std::string >	fVertexModuleLabel

std::vector< std::string >	fShowerModuleLabel

std::vector< std::string >	fCalorimetryModuleLabel

std::vector< std::string >	fParticleIDModuleLabel

std::vector< std::string >	fMVAPIDShowerModuleLabel

std::vector< std::string >	fMVAPIDTrackModuleLabel

std::vector< std::string >	fFlashT0FinderLabel

std::vector< std::string >	fMCT0FinderLabel

std::string	fPOTModuleLabel

std::string	fCosmicClusterTaggerAssocLabel

bool	fUseBuffer
	whether to use a permanent buffer (faster, huge memory) More...

bool	fSaveAuxDetInfo
	whether to extract and save auxiliary detector data More...

bool	fSaveCryInfo

bool	fSaveGenieInfo
	whether to extract and save CRY particle data More...

bool	fSaveProtoInfo
	whether to extract and save Genie information More...

bool	fSaveGeantInfo
	whether to extract and save ProtDUNE beam simulation information More...

bool	fSaveMCShowerInfo
	whether to extract and save Geant information More...

bool	fSaveMCTrackInfo
	whether to extract and save MC Shower information More...

bool	fSaveHitInfo
	whether to extract and save MC Track information More...

bool	fSaveRawDigitInfo
	whether to extract and save Hit information More...

bool	fSaveTrackInfo
	whether to extract and save Raw Digit information More...

bool	fSaveVertexInfo
	whether to extract and save Track information More...

bool	fSaveClusterInfo
	whether to extract and save Vertex information More...

bool	fSavePandoraNuVertexInfo
	whether to extract and save Cluster information More...

bool	fSaveFlashInfo
	whether to extract and save nu vertex information from Pandora More...

bool	fSaveExternCounterInfo
	whether to extract and save Flash information More...

bool	fSaveShowerInfo
	whether to extract and save External Counter information More...

bool	fSavePFParticleInfo
	whether to extract and save Shower information More...

bool	fSaveSpacePointSolverInfo
	whether to extract and save PFParticle information More...

bool	fSaveCnnInfo
	whether to extract and save SpacePointSolver information More...

std::vector< std::string >	fCosmicTaggerAssocLabel
	whether to extract and save CNN information More...

std::vector< std::string >	fContainmentTaggerAssocLabel

std::vector< std::string >	fFlashMatchAssocLabel

bool	bIgnoreMissingShowers
	whether to ignore missing shower information More...

bool	isCosmics
	if it contains cosmics More...

bool	fSaveCaloCosmics
	save calorimetry information for cosmics More...

float	fG4minE
	Energy threshold to save g4 particle info. More...

double	ActiveBounds [6]

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

Creates a simple ROOT tree with tracking and calorimetry information.

Configuration parameters

UseBuffers (default: false): if enabled, memory is allocated for tree data for all the run; otherwise, it's allocated on each event, used and freed; use "true" for speed, "false" to save memory
SaveAuxDetInfo (default: false): if enabled, auxiliary detector data will be extracted and included in the tree

Definition at line 1213 of file AnalysisTree_module.cc.

Constructor & Destructor Documentation

AnalysisTree::AnalysisTree ( fhicl::ParameterSet const & pset )

explicit

Definition at line 3326 of file AnalysisTree_module.cc.

                                                             :
   EDAnalyzer(pset),
   fTree(nullptr), fPOT(nullptr),
   fDigitModuleLabel         (pset.get< std::string >("DigitModuleLabel")        ),
   fHitsModuleLabel          (pset.get< std::string >("HitsModuleLabel")         ),
   fLArG4ModuleLabel         (pset.get< std::string >("LArGeantModuleLabel")     ),
   fCalDataModuleLabel       (pset.get< std::string >("CalDataModuleLabel")      ),
   fGenieGenModuleLabel      (pset.get< std::string >("GenieGenModuleLabel")     ),
   fCryGenModuleLabel        (pset.get< std::string >("CryGenModuleLabel")       ),
   fProtoGenModuleLabel      (pset.get< std::string >("ProtoGenModuleLabel")     ),
   fG4ModuleLabel            (pset.get< std::string >("G4ModuleLabel")           ),
   fClusterModuleLabel       (pset.get< std::string >("ClusterModuleLabel")     ),
   fPandoraNuVertexModuleLabel (pset.get< std::string >("PandoraNuVertexModuleLabel")     ),
   fOpFlashModuleLabel       (pset.get< std::string >("OpFlashModuleLabel")      ),
   fExternalCounterModuleLabel (pset.get< std::string >("ExternalCounterModuleLabel")      ),
   fMCShowerModuleLabel      (pset.get< std::string >("MCShowerModuleLabel")     ),
   fMCTrackModuleLabel      (pset.get< std::string >("MCTrackModuleLabel")     ),
   fSpacePointSolverModuleLabel (pset.get< std::string >("SpacePointSolverModuleLabel")),
   fCnnModuleLabel           (pset.get< std::string >("CnnModuleLabel")),
   fTrackModuleLabel         (pset.get< std::vector<std::string> >("TrackModuleLabel")),
   fVertexModuleLabel        (pset.get< std::vector<std::string> >("VertexModuleLabel")),
   fShowerModuleLabel        (pset.get< std::vector<std::string> >("ShowerModuleLabel")),
   fCalorimetryModuleLabel   (pset.get< std::vector<std::string> >("CalorimetryModuleLabel")),
   fParticleIDModuleLabel    (pset.get< std::vector<std::string> >("ParticleIDModuleLabel")   ),
   fMVAPIDShowerModuleLabel  (pset.get< std::vector<std::string> >("MVAPIDShowerModuleLabel")   ),
   fMVAPIDTrackModuleLabel   (pset.get< std::vector<std::string> >("MVAPIDTrackModuleLabel")   ),
   fFlashT0FinderLabel       (pset.get< std::vector<std::string> >("FlashT0FinderLabel")   ),
   fMCT0FinderLabel          (pset.get< std::vector<std::string> >("MCT0FinderLabel")   ),
   fPOTModuleLabel           (pset.get< std::string >("POTModuleLabel")),
   fCosmicClusterTaggerAssocLabel (pset.get< std::string >("CosmicClusterTaggerAssocLabel")),
   fUseBuffer                (pset.get< bool >("UseBuffers", false)),
   fSaveAuxDetInfo           (pset.get< bool >("SaveAuxDetInfo", false)),
   fSaveCryInfo              (pset.get< bool >("SaveCryInfo", false)),
   fSaveGenieInfo            (pset.get< bool >("SaveGenieInfo", false)),
   fSaveProtoInfo            (pset.get< bool >("SaveProtoInfo", false)),
   fSaveGeantInfo            (pset.get< bool >("SaveGeantInfo", false)),
   fSaveMCShowerInfo         (pset.get< bool >("SaveMCShowerInfo", false)),
   fSaveMCTrackInfo          (pset.get< bool >("SaveMCTrackInfo", false)),
   fSaveHitInfo              (pset.get< bool >("SaveHitInfo", false)),
   fSaveRawDigitInfo                 (pset.get< bool >("SaveRawDigitInfo", false)),
   fSaveTrackInfo            (pset.get< bool >("SaveTrackInfo", false)),
   fSaveVertexInfo           (pset.get< bool >("SaveVertexInfo", false)),
   fSaveClusterInfo          (pset.get< bool >("SaveClusterInfo", false)),
   fSavePandoraNuVertexInfo        (pset.get< bool >("SavePandoraNuVertexInfo", false)),
   fSaveFlashInfo            (pset.get< bool >("SaveFlashInfo", false)),
   fSaveExternCounterInfo            (pset.get< bool >("SaveExternCounterInfo", false)),
   fSaveShowerInfo            (pset.get< bool >("SaveShowerInfo", false)),
   fSavePFParticleInfo       (pset.get< bool >("SavePFParticleInfo", false)),
   fSaveSpacePointSolverInfo (pset.get< bool >("SaveSpacePointSolverInfo", false)),
   fSaveCnnInfo              (pset.get< bool >("SaveCnnInfo", false)),
   fCosmicTaggerAssocLabel  (pset.get<std::vector< std::string > >("CosmicTaggerAssocLabel") ),
   fContainmentTaggerAssocLabel  (pset.get<std::vector< std::string > >("ContainmentTaggerAssocLabel") ),
   fFlashMatchAssocLabel (pset.get<std::vector< std::string > >("FlashMatchAssocLabel") ),
   bIgnoreMissingShowers     (pset.get< bool >("IgnoreMissingShowers", false)),
   isCosmics(false),
   fSaveCaloCosmics          (pset.get< bool >("SaveCaloCosmics",false)),
   fG4minE                   (pset.get< float>("G4minE",0.01))
 {
 
   if (fSavePFParticleInfo) fPFParticleModuleLabel = pset.get<std::string>("PFParticleModuleLabel");
 
   if (fSaveAuxDetInfo == true) fSaveGeantInfo = true;
   if (fSaveRawDigitInfo == true) fSaveHitInfo = true;
   mf::LogInfo("AnalysisTree") << "Configuration:"
                               << "\n  UseBuffers: " << std::boolalpha << fUseBuffer
     ;
   if (GetNTrackers() > kMaxTrackers) {
     throw art::Exception(art::errors::Configuration)
       << "AnalysisTree currently supports only up to " << kMaxTrackers
       << " tracking algorithms, but " << GetNTrackers() << " are specified."
       << "\nYou can increase kMaxTrackers and recompile.";
   } // if too many trackers
   if (fTrackModuleLabel.size() != fCalorimetryModuleLabel.size()){
     throw art::Exception(art::errors::Configuration)
       << "fTrackModuleLabel.size() = "<<fTrackModuleLabel.size()<<" does not match "
       << "fCalorimetryModuleLabel.size() = "<<fCalorimetryModuleLabel.size();
   }
   if (fTrackModuleLabel.size() != fParticleIDModuleLabel.size()){
     throw art::Exception(art::errors::Configuration)
       << "fTrackModuleLabel.size() = "<<fTrackModuleLabel.size()<<" does not match "
       << "fParticleIDModuleLabel.size() = "<<fParticleIDModuleLabel.size();
   }
   if (fTrackModuleLabel.size() != fFlashT0FinderLabel.size()){
     throw art::Exception(art::errors::Configuration)
       << "fTrackModuleLabel.size() = "<<fTrackModuleLabel.size()<<" does not match "
       << "fFlashT0FinderLabel.size() = "<<fFlashT0FinderLabel.size();
   }
   if (fTrackModuleLabel.size() != fMCT0FinderLabel.size()){
     throw art::Exception(art::errors::Configuration)
       << "fTrackModuleLabel.size() = "<<fTrackModuleLabel.size()<<" does not match "
       << "fMCT0FinderLabel.size() = "<<fMCT0FinderLabel.size();
   }
   if (fTrackModuleLabel.size() != fCosmicTaggerAssocLabel.size()) {
     throw art::Exception(art::errors::Configuration)
       << "fTrackModuleLabel.size() = "<<fTrackModuleLabel.size()<<" does not match "
       << "fCosmicTaggerAssocLabel.size() = "<<fCosmicTaggerAssocLabel.size();
   }
   if (fTrackModuleLabel.size() != fContainmentTaggerAssocLabel.size()) {
     throw art::Exception(art::errors::Configuration)
       << "fTrackModuleLabel.size() = "<<fTrackModuleLabel.size()<<" does not match "
       << "fCosmicTaggerAssocLabel.size() = "<<fContainmentTaggerAssocLabel.size();
   }
   if (fTrackModuleLabel.size() != fFlashMatchAssocLabel.size()) {
     throw art::Exception(art::errors::Configuration)
       << "fTrackModuleLabel.size() = "<<fTrackModuleLabel.size()<<" does not match "
       << "fCosmicTaggerAssocLabel.size() = "<<fFlashMatchAssocLabel.size();
   }
   if (GetNVertexAlgos() > kMaxVertexAlgos) {
     throw art::Exception(art::errors::Configuration)
       << "AnalysisTree currently supports only up to " << kMaxVertexAlgos
       << " tracking algorithms, but " << GetNVertexAlgos() << " are specified."
       << "\nYou can increase kMaxVertexAlgos and recompile.";
   } // if too many trackers
 
   // Build my Cryostat boundaries array...Taken from Tyler Alion in Geometry Core. Should still return the same values for uBoone.
   ActiveBounds[0] = ActiveBounds[2] = ActiveBounds[4] = DBL_MAX;
   ActiveBounds[1] = ActiveBounds[3] = ActiveBounds[5] = -DBL_MAX;
   // assume single cryostats
   auto const* geom = lar::providerFrom<geo::Geometry>();
   for (geo::TPCGeo const& TPC: geom->IterateTPCs()) {
     // get center in world coordinates
     double origin[3] = {0.};
     double center[3] = {0.};
     TPC.LocalToWorld(origin, center);
     double tpcDim[3] = {TPC.HalfWidth(), TPC.HalfHeight(), 0.5*TPC.Length() };
 
     if( center[0] - tpcDim[0] < ActiveBounds[0] ) ActiveBounds[0] = center[0] - tpcDim[0];
     if( center[0] + tpcDim[0] > ActiveBounds[1] ) ActiveBounds[1] = center[0] + tpcDim[0];
     if( center[1] - tpcDim[1] < ActiveBounds[2] ) ActiveBounds[2] = center[1] - tpcDim[1];
     if( center[1] + tpcDim[1] > ActiveBounds[3] ) ActiveBounds[3] = center[1] + tpcDim[1];
     if( center[2] - tpcDim[2] < ActiveBounds[4] ) ActiveBounds[4] = center[2] - tpcDim[2];
     if( center[2] + tpcDim[2] > ActiveBounds[5] ) ActiveBounds[5] = center[2] + tpcDim[2];
   } // for all TPC
   std::cout << "Active Boundaries: "
             << "\n\tx: " << ActiveBounds[0] << " to " << ActiveBounds[1]
             << "\n\ty: " << ActiveBounds[2] << " to " << ActiveBounds[3]
             << "\n\tz: " << ActiveBounds[4] << " to " << ActiveBounds[5]
             << std::endl;
 } // dune::AnalysisTree::AnalysisTree()

AnalysisTree::~AnalysisTree ( )

virtual

Definition at line 3467 of file AnalysisTree_module.cc.

 {
   DestroyData();
 }

Member Function Documentation

void AnalysisTree::analyze ( const art::Event & evt )

read access to event

transfer the run and subrun data to the tree data object

Definition at line 3538 of file AnalysisTree_module.cc.

 {
   std::cout << "Analysing.\n\n";
   //services
   art::ServiceHandle<cheat::BackTrackerService> bt_serv;
   art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
 
   // collect the sizes which might me needed to resize the tree data structure:
   bool isMC = !evt.isRealData();
 
   // * hits
   std::vector<art::Ptr<recob::Hit> > hitlist;
   auto hitListHandle = evt.getHandle< std::vector<recob::Hit> >(fHitsModuleLabel);
   if (hitListHandle)
     art::fill_ptr_vector(hitlist, hitListHandle);
 
   // * clusters
   art::Handle< std::vector<recob::Cluster> > clusterListHandle;
   std::vector<art::Ptr<recob::Cluster> > clusterlist;
   if (fSaveClusterInfo){
     clusterListHandle = evt.getHandle< std::vector<recob::Cluster> >(fClusterModuleLabel);
     if (clusterListHandle)
       art::fill_ptr_vector(clusterlist, clusterListHandle);
   }
 
   // * spacepoints
   art::Handle<std::vector<recob::SpacePoint>> spacepointListHandle;
   art::Handle<std::vector<recob::PointCharge>> pointchargeListHandle;
   if (fSaveSpacePointSolverInfo) {
     std::cout << "Saving SpacepointSolver info.";
     spacepointListHandle = evt.getHandle<std::vector<recob::SpacePoint>>(fSpacePointSolverModuleLabel);
     pointchargeListHandle = evt.getHandle<std::vector<recob::PointCharge>>(fSpacePointSolverModuleLabel);
     if (spacepointListHandle->size() != pointchargeListHandle->size()) {
       throw cet::exception("tutorial::ReadSpacePointAndCnn")
           << "size of point and charge containers must be equal" << std::endl;
     }
   }
 
   // * flashes
   std::vector<art::Ptr<recob::OpFlash> > flashlist;
   auto flashListHandle = evt.getHandle< std::vector<recob::OpFlash> >(fOpFlashModuleLabel);
   if (flashListHandle)
     art::fill_ptr_vector(flashlist, flashListHandle);
 
   // * External Counters
   std::vector<art::Ptr<raw::ExternalTrigger> > countlist;
   auto countListHandle = evt.getHandle< std::vector<raw::ExternalTrigger> >(fExternalCounterModuleLabel);
   if (countListHandle)
     art::fill_ptr_vector(countlist, countListHandle);
 
   // * MC truth information
   art::Handle< std::vector<simb::MCTruth> > mctruthListHandle;
   std::vector<art::Ptr<simb::MCTruth> > mclist;
   if (isMC){
     mctruthListHandle = evt.getHandle< std::vector<simb::MCTruth> >(fGenieGenModuleLabel);
     if (mctruthListHandle)
       art::fill_ptr_vector(mclist, mctruthListHandle);
   }
 
   // *MC truth cosmic generator information
   art::Handle< std::vector<simb::MCTruth> > mctruthcryListHandle;
   std::vector<art::Ptr<simb::MCTruth> > mclistcry;
   if (isMC && fSaveCryInfo){
     mctruthcryListHandle = evt.getHandle< std::vector<simb::MCTruth> >(fCryGenModuleLabel);
     if (mctruthcryListHandle) {
       art::fill_ptr_vector(mclistcry, mctruthcryListHandle);
     }
     else{
       // If we requested this info but it doesn't exist, don't use it.
       fSaveCryInfo = false;
       mf::LogError("AnalysisTree") << "Requested CRY information but none exists, hence not saving CRY information.";
     }
   }
 
   // ProtoDUNE beam generator information
   art::Handle< std::vector<simb::MCTruth> > mctruthprotoListHandle;
   std::vector<art::Ptr<simb::MCTruth> > mclistproto;
   if (isMC && fSaveProtoInfo){
     mctruthprotoListHandle = evt.getHandle< std::vector<simb::MCTruth> >(fProtoGenModuleLabel);
     if (mctruthprotoListHandle) {
       art::fill_ptr_vector(mclistproto,mctruthprotoListHandle);
     }
     else{
       // If we requested this info but it doesn't exist, don't use it.
       fSaveProtoInfo = false;
       mf::LogError("AnalysisTree") << "Requested protoDUNE beam generator information but none exists, hence not saving protoDUNE beam generator information.";
     }
   }
 
   // *MC Shower information
   art::Handle< std::vector<sim::MCShower> > mcshowerh;
   if (isMC)
     mcshowerh = evt.getHandle< std::vector<sim::MCShower> >(fMCShowerModuleLabel);
 
   int nMCShowers = 0;
   if (fSaveMCShowerInfo && mcshowerh.isValid())
     nMCShowers = mcshowerh->size();
 
   // *MC Track information
   art::Handle< std::vector<sim::MCTrack> > mctrackh;
   if (isMC)
     mctrackh = evt.getHandle< std::vector<sim::MCTrack> >(fMCTrackModuleLabel);
 
   int nMCTracks = 0;
   if (fSaveMCTrackInfo && mctrackh.isValid())
     nMCTracks = mctrackh->size();
 
   art::Ptr<simb::MCTruth> mctruthcry;
   int nCryPrimaries = 0;
 
   if (fSaveCryInfo){
     mctruthcry = mclistcry[0];
     nCryPrimaries = mctruthcry->NParticles();
   }
 
   art::Ptr<simb::MCTruth> mctruthproto;
   int nProtoPrimaries = 0;
   if( fSaveProtoInfo){
     mctruthproto = mclistproto[0];
     nProtoPrimaries = mctruthproto->NParticles();
   }
 
   int nGeniePrimaries = 0, nGEANTparticles = 0;
 
   art::Ptr<simb::MCTruth> mctruth;
 
 
   if (isMC) { //is MC
 
     //find origin
     auto const allmclists = evt.getMany<std::vector<simb::MCTruth>>();
     for(size_t mcl = 0; mcl < allmclists.size(); ++mcl){
       art::Handle< std::vector<simb::MCTruth> > mclistHandle = allmclists[mcl];
       for(size_t m = 0; m < mclistHandle->size(); ++m){
         art::Ptr<simb::MCTruth> mct(mclistHandle, m);
         if (mct->Origin() == simb::kCosmicRay) isCosmics = true;
       }
     }
     if (fSaveCaloCosmics) isCosmics = false; //override to save calo info
 
     // GENIE
     if (!mclist.empty()){//at least one mc record
 
       //        double maxenergy = -1;
       //        int imc0 = 0;
       //        for (std::map<art::Ptr<simb::MCTruth>,double>::iterator ii=mctruthemap.begin(); ii!=mctruthemap.end(); ++ii){
       //          if ((ii->second)>maxenergy){
       //            maxenergy = ii->second;
       //            mctruth = ii->first;
       //            imc = imc0;
       //          }
       //          imc0++;
       //        }
 
       //imc = 0; //set imc to 0 to solve a confusion for BNB+cosmic files where there are two MCTruth
       mctruth = mclist[0];
 
       if (mctruth->NeutrinoSet()) nGeniePrimaries = mctruth->NParticles();
       //} //end (fSaveGenieInfo)
 
       const sim::ParticleList& plist = pi_serv->ParticleList();
       nGEANTparticles = plist.size();
 
       // to know the number of particles in AV would require
       // looking at all of them; so we waste some memory here
     } // if have MC truth
     MF_LOG_DEBUG("AnalysisTree") << "Expected "
                               << nGEANTparticles << " GEANT particles, "
                               << nGeniePrimaries << " GENIE particles";
   } // if MC
 
   // SpacePointSolver information
   int nSpacePoints = 0;
   if (fSaveSpacePointSolverInfo && spacepointListHandle.isValid())
     nSpacePoints = spacepointListHandle->size();
 
   CreateData(); // tracker data is created with default constructor
   if (fSaveGenieInfo)
     fData->ResizeGenie(nGeniePrimaries);
   if (fSaveCryInfo)
     fData->ResizeCry(nCryPrimaries);
   if (fSaveProtoInfo)
     fData->ResizeProto(nProtoPrimaries);
   if (fSaveGeantInfo)
     fData->ResizeGEANT(nGEANTparticles);
   if (fSaveMCShowerInfo)
     fData->ResizeMCShower(nMCShowers);
   if (fSaveMCTrackInfo)
     fData->ResizeMCTrack(nMCTracks);
   if (fSaveSpacePointSolverInfo)
     fData->ResizeSpacePointSolver(nSpacePoints);
 
   fData->ClearLocalData(); // don't bother clearing tracker data yet
 
   const size_t NTrackers = GetNTrackers(); // number of trackers passed into fTrackModuleLabel
   const size_t NShowerAlgos = GetNShowerAlgos(); // number of shower algorithms into fShowerModuleLabel
   const size_t NHits     = hitlist.size(); // number of hits
   const size_t NVertexAlgos = GetNVertexAlgos(); // number of vertex algos
   const size_t NClusters = clusterlist.size(); //number of clusters
   const size_t NFlashes  = flashlist.size(); // number of flashes
   const size_t NExternCounts = countlist.size(); // number of External Counters
   // make sure there is the data, the tree and everything;
   CreateTree();
 
   /// transfer the run and subrun data to the tree data object
   //  fData->RunData = RunData;
   fData->SubRunData = SubRunData;
 
   fData->isdata = int(!isMC);
 
   // * raw trigger
   std::vector<art::Ptr<raw::Trigger>> triggerlist;
   auto triggerListHandle = evt.getHandle< std::vector<raw::Trigger>>(fDigitModuleLabel);
   if (triggerListHandle)
     art::fill_ptr_vector(triggerlist, triggerListHandle);
 
   if (triggerlist.size()){
     fData->triggernumber = triggerlist[0]->TriggerNumber();
     fData->triggertime   = triggerlist[0]->TriggerTime();
     fData->beamgatetime  = triggerlist[0]->BeamGateTime();
     fData->triggerbits   = triggerlist[0]->TriggerBits();
   }
 
   // * vertices
   std::vector< art::Handle< std::vector<recob::Vertex> > > vertexListHandle(NVertexAlgos);
   std::vector< std::vector<art::Ptr<recob::Vertex> > > vertexlist(NVertexAlgos);
   for (unsigned int it = 0; it < NVertexAlgos; ++it){
     vertexListHandle[it] = evt.getHandle< std::vector<recob::Vertex> >(fVertexModuleLabel[it]);
     if (vertexListHandle[it])
       art::fill_ptr_vector(vertexlist[it], vertexListHandle[it]);
   }
 
   // * PFParticles
   lar_pandora::PFParticleVector pfparticlelist;
   lar_pandora::PFParticlesToClusters pfParticleToClusterMap;
   lar_pandora::LArPandoraHelper::CollectPFParticles(evt, fPFParticleModuleLabel, pfparticlelist, pfParticleToClusterMap);
 
   // * tracks
   std::vector< art::Handle< std::vector<recob::Track> > > trackListHandle(NTrackers);
   std::vector< std::vector<art::Ptr<recob::Track> > > tracklist(NTrackers);
   for (unsigned int it = 0; it < NTrackers; ++it){
     trackListHandle[it] = evt.getHandle< std::vector<recob::Track> >(fTrackModuleLabel[it]);
     if (trackListHandle[it])
       art::fill_ptr_vector(tracklist[it], trackListHandle[it]);
   }
 
   // * showers
   // It seems that sometimes the shower data product does not exist;
   // in that case, we store a nullptr in place of the pointer to the vector;
   // the data structure itself is owned by art::Event and we should not try
   // to manage its memory
   std::vector<std::vector<recob::Shower> const*> showerList;
   std::vector< art::Handle< std::vector<recob::Shower> > > showerListHandle;
   showerList.reserve(fShowerModuleLabel.size());
   if (fSaveShowerInfo) {
     for (art::InputTag ShowerInputTag: fShowerModuleLabel) {
       auto ShowerHandle = evt.getHandle<std::vector<recob::Shower>>(ShowerInputTag);
       if (!ShowerHandle) {
         showerList.push_back(nullptr);
         if (!bIgnoreMissingShowers) {
           throw art::Exception(art::errors::ProductNotFound)
             << "Showers with input tag '" << ShowerInputTag.encode()
             << "' were not found in the event."
             " If you really know what you are doing,"
             " set AnalysisTree's configuration parameter IgnoreMissingShowers"
             " to \"true\"; the lack of any shower set will be tolerated,"
             " and the shower list in the corresponding event will be set to"
             " a list of one shower, with an invalid ID.\n";
         } // if bIgnoreMissingShowers
         else {
           // this message is more alarming than strictly necessary; by design.
           mf::LogError("AnalysisTree")
             << "No showers found for input tag '" << ShowerInputTag.encode()
             << "' ; FILLING WITH FAKE DATA AS FOR USER'S REQUEST";
         }
       }
 
       else showerList.push_back(ShowerHandle.product());
 
       showerListHandle.push_back(ShowerHandle); // either way, put it into the handle list
 
     }
 
   } // for shower input tag
 
 
 
   std::vector<const sim::AuxDetSimChannel*> fAuxDetSimChannels;
   if (fSaveAuxDetInfo){
     evt.getView(fLArG4ModuleLabel, fAuxDetSimChannels);
   }
 
   std::vector<const sim::SimChannel*> fSimChannels;
   if (isMC && fSaveGeantInfo)
     evt.getView(fLArG4ModuleLabel, fSimChannels);
 
   fData->run = evt.run();
   fData->subrun = evt.subRun();
   fData->event = evt.id().event();
 
   art::Timestamp ts = evt.time();
   TTimeStamp tts(ts.timeHigh(), ts.timeLow());
   fData->evttime = tts.AsDouble();
 
   //copied from MergeDataPaddles.cxx
   auto beam = evt.getHandle< raw::BeamInfo >("beamdata");
   if (beam){
     fData->beamtime = (double)beam->get_t_ms();
     fData->beamtime/=1000.; //in second
     std::map<std::string, std::vector<double>> datamap = beam->GetDataMap();
     if (datamap["E:TOR860"].size()){
       fData->potbnb = datamap["E:TOR860"][0];
     }
     if (datamap["E:TORTGT"].size()){
       fData->potnumitgt = datamap["E:TORTGT"][0];
     }
     if (datamap["E:TOR101"].size()){
       fData->potnumi101 = datamap["E:TOR101"][0];
     }
   }
 
 
   //  std::cout<<detProp.NumberTimeSamples()<<" "<<detProp.ReadOutWindowSize()<<std::endl;
   //  std::cout<<geom->DetHalfHeight()*2<<" "<<geom->DetHalfWidth()*2<<" "<<geom->DetLength()<<std::endl;
   //  std::cout<<geom->Nwires(0)<<" "<<geom->Nwires(1)<<" "<<geom->Nwires(2)<<std::endl;
   auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
   auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clockData);
 
   //hit information
   if (fSaveHitInfo){
     fData->no_hits = (int) NHits;
     fData->no_hits_stored = TMath::Min( (int) NHits, (int) kMaxHits);
     if (NHits > kMaxHits) {
       // got this error? consider increasing kMaxHits
       // (or ask for a redesign using vectors)
       mf::LogError("AnalysisTree:limits") << "event has " << NHits
                                           << " hits, only kMaxHits=" << kMaxHits << " stored in tree";
     }
     for (size_t i = 0; i < NHits && i < kMaxHits ; ++i){//loop over hits
       fData->hit_channel[i] = hitlist[i]->Channel();
       fData->hit_tpc[i]   = hitlist[i]->WireID().TPC;
       fData->hit_plane[i]   = hitlist[i]->WireID().Plane;
       fData->hit_wire[i]    = hitlist[i]->WireID().Wire;
       fData->hit_peakT[i]   = hitlist[i]->PeakTime();
       fData->hit_charge[i]  = hitlist[i]->Integral();
       fData->hit_ph[i]  = hitlist[i]->PeakAmplitude();
       fData->hit_startT[i] = hitlist[i]->PeakTimeMinusRMS();
       fData->hit_endT[i] = hitlist[i]->PeakTimePlusRMS();
       fData->hit_rms[i] = hitlist[i]->RMS();
       fData->hit_goodnessOfFit[i] = hitlist[i]->GoodnessOfFit();
       fData->hit_multiplicity[i] = hitlist[i]->Multiplicity();
       //std::vector<double> xyz = bt_serv->HitToXYZ(hitlist[i]);
       //when the size of simIDEs is zero, the above function throws an exception
       //and crashes, so check that the simIDEs have non-zero size before
       //extracting hit true XYZ from simIDEs
       if (isMC){
         std::vector<const sim::IDE*> ides;
         try{
           ides= bt_serv->HitToSimIDEs_Ps(clockData, hitlist[i]);
         }
         catch(...){}
           if (ides.size()>0){
             std::vector<double> xyz = bt_serv->SimIDEsToXYZ(ides);
             fData->hit_trueX[i] = xyz[0];
           }
         }
 
       /*
         for (unsigned int it=0; it<fTrackModuleLabel.size();++it){
         art::FindManyP<recob::Track> fmtk(hitListHandle,evt,fTrackModuleLabel[it]);
         if (fmtk.at(i).size()!=0){
         hit_trkid[it][i] = fmtk.at(i)[0]->ID();
         }
         else
         hit_trkid[it][i] = 0;
         }
       */
 
       if (fSaveRawDigitInfo){
         //Hit to RawDigit information
         art::FindManyP<raw::RawDigit> fmrd(hitListHandle,evt,fHitsModuleLabel);
         if (hitlist[i]->WireID().Plane==2)
           {
             int dataSize = fmrd.at(i)[0]->Samples();
             short ped = fmrd.at(i)[0]->GetPedestal();
 
             std::vector<short> rawadc(dataSize);
             raw::Uncompress(fmrd.at(i)[0]->ADCs(), rawadc, fmrd.at(i)[0]->Compression());
             int t0 = hitlist[i]->PeakTime() - 3*(hitlist[i]->RMS());
             if (t0<0) t0 = 0;
             int t1 = hitlist[i]->PeakTime() + 3*(hitlist[i]->RMS());
             if (t1>=dataSize) t1 = dataSize-1;
             fData->rawD_ph[i] = -1;
             fData->rawD_peakT[i] = -1;
             for (int j = t0; j<=t1; ++j){
               if (rawadc[j]-ped>fData->rawD_ph[i]){
                 fData->rawD_ph[i] = rawadc[j]-ped;
                 fData->rawD_peakT[i] = j;
               }
             }
             fData->rawD_charge[i] = 0;
             fData->rawD_fwhh[i] = 0;
             double mean_t = 0.0;
             double mean_t2 = 0.0;
             for (int j = t0; j<=t1; ++j){
               if (rawadc[j]-ped>=0.5*fData->rawD_ph[i]){
                 ++fData->rawD_fwhh[i];
               }
               if (rawadc[j]-ped>=0.1*fData->rawD_ph[i]){
                 fData->rawD_charge[i] += rawadc[j]-ped;
                 mean_t += (double)j*(rawadc[j]-ped);
                 mean_t2 += (double)j*(double)j*(rawadc[j]-ped);
               }
             }
             mean_t/=fData->rawD_charge[i];
             mean_t2/=fData->rawD_charge[i];
             fData->rawD_rms[i] = sqrt(mean_t2-mean_t*mean_t);
           }
       } // Save RawDigitInfo
 
       if (!evt.isRealData()&&!isCosmics){
         fData -> hit_nelec[i] = 0;
         fData -> hit_energy[i] = 0;
         const sim::SimChannel* chan = 0;
         for(size_t sc = 0; sc < fSimChannels.size(); ++sc){
           if(fSimChannels[sc]->Channel() == hitlist[i]->Channel()) chan = fSimChannels[sc];
         }
         if (chan){
     for(auto const& mapitr : chan->TDCIDEMap()){
             // loop over the vector of IDE objects.
       for(auto const& ide : mapitr.second){
               fData -> hit_nelec[i] += ide.numElectrons;
               fData -> hit_energy[i] += ide.energy;
             }
           }
         }
       }
     } // Loop over hits
 
     hitListHandle = evt.getHandle< std::vector<recob::Hit> >(fHitsModuleLabel);
     if (hitListHandle) {
       //Find tracks associated with hits
       art::FindManyP<recob::Track> fmtk(hitListHandle,evt,fTrackModuleLabel[0]);
       for (size_t i = 0; i < NHits && i < kMaxHits ; ++i){//loop over hits
         if (fmtk.isValid()){
           if (fmtk.at(i).size()!=0){
             fData->hit_trkid[i] = fmtk.at(i)[0]->ID();
             fData->hit_trkKey[i] = fmtk.at(i)[0].key();
 
           }
           else
             fData->hit_trkid[i] = -1;
         }
       }
     }
 
     //In the case of ClusterCrawler or linecluster, use "linecluster or clustercrawler" as HitModuleLabel.
     //using cchit will not make this association. In the case of gaushit, just use gaushit
     //Not initializing clusterID to -1 since some clustering algorithms assign negative IDs!
 
     hitListHandle = evt.getHandle< std::vector<recob::Hit> >(fHitsModuleLabel);
     if (hitListHandle) {
       //Find clusters and spacepoints associated with hits
       art::FindManyP<recob::Cluster> fmcl(hitListHandle,evt,fClusterModuleLabel);
       art::FindManyP<recob::SpacePoint> fmsp(hitListHandle,evt,fSpacePointSolverModuleLabel);
       for (size_t i = 0; i < NHits && i < kMaxHits ; ++i){//loop over hits
         if (fmcl.isValid() && fmcl.at(i).size()!=0){
           fData->hit_clusterid[i] = fmcl.at(i)[0]->ID();
           fData->hit_clusterKey[i] = fmcl.at(i)[0].key();
           // std::cout << "ClusterID " <<
         }
         if(fmsp.isValid() && fmsp.at(i).size()!=0){
           fData->hit_spacepointid[i] = fmsp.at(i)[0]->ID();
           fData->hit_spacepointKey[i] = fmsp.at(i)[0].key();
         }
       }
     }
   }// end (fSaveHitInfo)
 
 
   if(fSavePandoraNuVertexInfo) {
     lar_pandora::PFParticleVector particleVector;
     lar_pandora::LArPandoraHelper::CollectPFParticles(evt, fPandoraNuVertexModuleLabel, particleVector);
     lar_pandora::VertexVector vertexVector;
     lar_pandora::PFParticlesToVertices particlesToVertices;
     lar_pandora::LArPandoraHelper::CollectVertices(evt, fPandoraNuVertexModuleLabel, vertexVector, particlesToVertices);
 
     short nprim = 0;
     for (unsigned int n = 0; n < particleVector.size(); ++n) {
       const art::Ptr<recob::PFParticle> particle = particleVector.at(n);
       if(particle->IsPrimary()) nprim++;
     }
 
     if (nprim > kMaxVertices){
       // got this error? consider increasing kMaxClusters
       // (or ask for a redesign using vectors)
       mf::LogError("AnalysisTree:limits") << "event has " << nprim
                                           << " nu neutrino vertices, only kMaxVertices=" << kMaxVertices << " stored in tree";
     }
 
     fData->nnuvtx = nprim;
 
     short iv = 0;
     for (unsigned int n = 0; n < particleVector.size(); ++n) {
       const art::Ptr<recob::PFParticle> particle = particleVector.at(n);
       if(particle->IsPrimary()) {
         lar_pandora::PFParticlesToVertices::const_iterator vIter = particlesToVertices.find(particle);
         if (particlesToVertices.end() != vIter) {
           const lar_pandora::VertexVector &vertexVector = vIter->second;
           if (!vertexVector.empty()) {
             if (vertexVector.size() == 1) {
               const art::Ptr<recob::Vertex> vertex = *(vertexVector.begin());
               double xyz[3] = {0.0, 0.0, 0.0} ;
               vertex->XYZ(xyz);
               fData->nuvtxx[iv] = xyz[0];
               fData->nuvtxy[iv] = xyz[1];
               fData->nuvtxz[iv] = xyz[2];
               fData->nuvtxpdg[iv] = particle->PdgCode();
               iv++;
             }
           }
         }
       }
     }
   } // save PandoraNuVertexInfo
 
 
   if (fSaveClusterInfo){
     fData->nclusters = (int) NClusters;
     if (NClusters > kMaxClusters){
       // got this error? consider increasing kMaxClusters
       // (or ask for a redesign using vectors)
       mf::LogError("AnalysisTree:limits") << "event has " << NClusters
                                           << " clusters, only kMaxClusters=" << kMaxClusters << " stored in tree";
     }
     for(unsigned int ic=0; ic<NClusters;++ic){//loop over clusters
       art::Ptr<recob::Cluster> clusterholder(clusterListHandle, ic);
       const recob::Cluster& cluster = *clusterholder;
       fData->clusterId[ic] = cluster.ID();
       fData->clusterView[ic] = cluster.View();
       fData->cluster_isValid[ic] = cluster.isValid();
       fData->cluster_StartCharge[ic] = cluster.StartCharge();
       fData->cluster_StartAngle[ic] = cluster.StartAngle();
       fData->cluster_EndCharge[ic] = cluster.EndCharge();
       fData->cluster_EndAngle[ic] = cluster.EndAngle();
       fData->cluster_Integral[ic] = cluster.Integral();
       fData->cluster_IntegralAverage[ic] = cluster.IntegralAverage();
       fData->cluster_SummedADC[ic] = cluster.SummedADC();
       fData->cluster_SummedADCaverage[ic] = cluster.SummedADCaverage();
       fData->cluster_MultipleHitDensity[ic] = cluster.MultipleHitDensity();
       fData->cluster_Width[ic] = cluster.Width();
       fData->cluster_NHits[ic] = cluster.NHits();
       fData->cluster_StartWire[ic] = cluster.StartWire();
       fData->cluster_StartTick[ic] = cluster.StartTick();
       fData->cluster_EndWire[ic] = cluster.EndWire();
       fData->cluster_EndTick[ic] = cluster.EndTick();
 
       //Cosmic Tagger information for cluster
       art::FindManyP<anab::CosmicTag> fmcct(clusterListHandle,evt,fCosmicClusterTaggerAssocLabel);
       if (fmcct.isValid()){
         fData->cluncosmictags_tagger[ic]     = fmcct.at(ic).size();
         if (fmcct.at(ic).size()>0){
           if(fmcct.at(ic).size()>1)
             std::cerr << "\n Warning : more than one cosmic tag per cluster in module! assigning the first tag to the cluster" << fCosmicClusterTaggerAssocLabel;
           fData->clucosmicscore_tagger[ic] = fmcct.at(ic).at(0)->CosmicScore();
           fData->clucosmictype_tagger[ic] = fmcct.at(ic).at(0)->CosmicType();
         }
       }
     }//end loop over clusters
   }//end fSaveClusterInfo
 
   if (fSaveSpacePointSolverInfo){
     fData->nspacepoints = (unsigned int) nSpacePoints;
 
     // Largely copied from Robert Sulej's ReadSpacePointAndCnn_module.cc
     if(fSaveCnnInfo) {
       auto cluResults = anab::MVAReader<recob::Cluster, MVA_LENGTH>::create(evt, fCnnModuleLabel);
       if(cluResults) {
         size_t emLikeIdx = cluResults->getIndex("em"); // at which index EM-like is stored in CNN output vector
 
         const art::FindManyP<recob::Hit> hitsFromClusters(cluResults->dataHandle(), evt, cluResults->dataTag());
         const art::FindManyP<recob::SpacePoint> spFromHits(hitListHandle, evt, fSpacePointSolverModuleLabel);
 
         std::vector<size_t> sizeScore(nSpacePoints, 0); // keep track of the max size of a cluster containing hit associated to spacepoint
 
         for(size_t c = 0; c < cluResults->size(); ++c) {
           const std::vector< art::Ptr<recob::Hit> > & hits = hitsFromClusters.at(c);
           std::array<float, MVA_LENGTH> cnn_out = cluResults->getOutput(c);
 
           for (auto& hptr : hits) {
             const std::vector< art::Ptr<recob::SpacePoint> > & sp = spFromHits.at(hptr.key());
             for(const auto & spptr : sp) { // Should always be just one associated spacepoint
               if(hits.size() > sizeScore[spptr.key()]) {
                 sizeScore[spptr.key()] = hits.size();
                 fData->SpacePointEmScore[spptr.key()] = cnn_out[emLikeIdx];
               }
             } // Loop over associated spacepoints
           } // Loop over hits
         } // Loop over cluResults
       } // If cluResults
     }
 
     for (unsigned int is = 0; is < (unsigned int)nSpacePoints;
          ++is) {  // loop over spacepoints
       fData->SpacePointX[is] = (*spacepointListHandle)[is].XYZ()[0];
       fData->SpacePointY[is] = (*spacepointListHandle)[is].XYZ()[1];
       fData->SpacePointZ[is] = (*spacepointListHandle)[is].XYZ()[2];
 
       fData->SpacePointQ[is] = (*pointchargeListHandle)[is].charge();
 
       fData->SpacePointErrX[is] = (*spacepointListHandle)[is].ErrXYZ()[0];
       fData->SpacePointErrY[is] = (*spacepointListHandle)[is].ErrXYZ()[1];
       fData->SpacePointErrZ[is] = (*spacepointListHandle)[is].ErrXYZ()[2];
 
       fData->SpacePointID[is] = (*spacepointListHandle)[is].ID();
 
       fData->SpacePointID[is] = (*spacepointListHandle)[is].Chisq();
     }//end loop over spacepoints
   }//end fSpacePointSolverInfo
 
   if (fSaveFlashInfo){
     fData->no_flashes = (int) NFlashes;
     if (NFlashes > kMaxFlashes) {
       // got this error? consider increasing kMaxHits
       // (or ask for a redesign using vectors)
       mf::LogError("AnalysisTree:limits") << "event has " << NFlashes
                                           << " flashes, only kMaxFlashes=" << kMaxFlashes << " stored in tree";
     }
 
     std::sort(flashlist.begin(), flashlist.end(), recob::OpFlashPtrSortByPE);
 
     for (size_t i = 0; i < NFlashes && i < kMaxFlashes ; ++i){//loop over hits
       fData->flash_time[i]       = flashlist[i]->Time();
       fData->flash_pe[i]         = flashlist[i]->TotalPE();
       fData->flash_ycenter[i]    = flashlist[i]->YCenter();
       fData->flash_zcenter[i]    = flashlist[i]->ZCenter();
       fData->flash_ywidth[i]     = flashlist[i]->YWidth();
       fData->flash_zwidth[i]     = flashlist[i]->ZWidth();
       fData->flash_timewidth[i]  = flashlist[i]->TimeWidth();
     }
   }
 
   if (fSaveExternCounterInfo){
     fData->no_ExternCounts = (int) NExternCounts;
     if (NExternCounts > kMaxExternCounts) {
       // got this error? consider increasing kMaxHits
       // (or ask for a redesign using vectors)
       mf::LogError("AnalysisTree:limits") << "event has " << NExternCounts
                                           << " External Counters, only kMaxExternCounts=" << kMaxExternCounts << " stored in tree";
     }
     for (size_t i = 0; i < NExternCounts && i < kMaxExternCounts ; ++i){//loop over hits
       fData->externcounts_time[i] = countlist[i]->GetTrigTime();
       fData->externcounts_id[i]   = countlist[i]->GetTrigID();
     }
   }
 
 
   // Declare object-ID-to-PFParticleID maps so we can assign hasPFParticle and PFParticleID to the tracks, showers, vertices.
   std::map<Short_t, Short_t> trackIDtoPFParticleIDMap, vertexIDtoPFParticleIDMap, showerIDtoPFParticleIDMap;
 
   //Save PFParticle information
   if (fSavePFParticleInfo){
     AnalysisTreeDataStruct::PFParticleDataStruct& PFParticleData = fData->GetPFParticleData();
     size_t NPFParticles = pfparticlelist.size();
 
     PFParticleData.SetMaxPFParticles(std::max(NPFParticles, (size_t) 1));
     PFParticleData.Clear(); // clear all the data
 
     PFParticleData.nPFParticles = (short) NPFParticles;
 
     // now set the tree addresses to the newly allocated memory;
     // this creates the tree branches in case they are not there yet
     SetPFParticleAddress();
 
     if (NPFParticles > PFParticleData.GetMaxPFParticles()) {
       mf::LogError("AnalysisTree:limits") << "event has " << NPFParticles
                    << " PFParticles, only "
                    << PFParticleData.GetMaxPFParticles() << " stored in tree";
     }
 
     lar_pandora::PFParticleVector neutrinoPFParticles;
     lar_pandora::LArPandoraHelper::SelectNeutrinoPFParticles(pfparticlelist, neutrinoPFParticles);
     PFParticleData.pfp_numNeutrinos = neutrinoPFParticles.size();
 
     for (size_t i = 0; i < std::min(neutrinoPFParticles.size(), (size_t)kMaxNPFPNeutrinos); ++i) {
       PFParticleData.pfp_neutrinoIDs[i] = neutrinoPFParticles[i]->Self();
     }
 
     if (neutrinoPFParticles.size() > kMaxNPFPNeutrinos)
       std::cerr << "Warning: there were " << neutrinoPFParticles.size() << " reconstructed PFParticle neutrinos; only the first " << kMaxNPFPNeutrinos << " being stored in tree" << std::endl;
 
     // Get a PFParticle-to-vertex map.
     lar_pandora::VertexVector allPfParticleVertices;
     lar_pandora::PFParticlesToVertices pfParticleToVertexMap;
     lar_pandora::LArPandoraHelper::CollectVertices(evt, fPFParticleModuleLabel, allPfParticleVertices, pfParticleToVertexMap);
 
     // Get a PFParticle-to-track map.
     lar_pandora::TrackVector allPfParticleTracks;
     lar_pandora::PFParticlesToTracks pfParticleToTrackMap;
     lar_pandora::LArPandoraHelper::CollectTracks(evt, fPFParticleModuleLabel, allPfParticleTracks, pfParticleToTrackMap);
 
     // Get a PFParticle-to-shower map.
     lar_pandora::ShowerVector allPfParticleShowers;
     lar_pandora::PFParticlesToShowers pfParticleToShowerMap;
     lar_pandora::LArPandoraHelper::CollectShowers(evt, fPFParticleModuleLabel, allPfParticleShowers, pfParticleToShowerMap);
 
     for (size_t i = 0; i < NPFParticles && i < PFParticleData.GetMaxPFParticles() ; ++i){
       PFParticleData.pfp_selfID[i] = pfparticlelist[i]->Self();
       PFParticleData.pfp_isPrimary[i] = (Short_t)pfparticlelist[i]->IsPrimary();
       PFParticleData.pfp_numDaughters[i] = pfparticlelist[i]->NumDaughters();
       PFParticleData.pfp_parentID[i] = pfparticlelist[i]->Parent();
       PFParticleData.pfp_pdgCode[i] = pfparticlelist[i]->PdgCode();
       PFParticleData.pfp_isNeutrino[i] = lar_pandora::LArPandoraHelper::IsNeutrino(pfparticlelist[i]);
 
       // Set the daughter IDs.
       std::vector<size_t> daughterIDs = pfparticlelist[i]->Daughters();
 
       for (size_t j = 0; j < daughterIDs.size(); ++j)
         PFParticleData.pfp_daughterIDs[i][j] = daughterIDs[j];
 
       // Set the vertex ID.
       auto vertexMapIter = pfParticleToVertexMap.find(pfparticlelist[i]);
       if (vertexMapIter != pfParticleToVertexMap.end()) {
           lar_pandora::VertexVector pfParticleVertices = vertexMapIter->second;
 
           if (pfParticleVertices.size() > 1)
             std::cerr << "Warning: there was more than one vertex found for PFParticle with ID " << pfparticlelist[i]->Self() << ", storing only one" << std::endl;
 
           if (pfParticleVertices.size() > 0) {
             PFParticleData.pfp_vertexID[i] = pfParticleVertices.at(0)->ID();
             vertexIDtoPFParticleIDMap.insert(std::make_pair(pfParticleVertices.at(0)->ID(), pfparticlelist[i]->Self()));
           }
       }
       else
         std::cerr << "Warning: there was no vertex found for PFParticle with ID " << pfparticlelist[i]->Self() << std::endl;
 
       if (lar_pandora::LArPandoraHelper::IsTrack(pfparticlelist[i])){
         PFParticleData.pfp_isTrack[i] = 1;
 
         // Set the track ID.
         auto trackMapIter = pfParticleToTrackMap.find(pfparticlelist[i]);
         if (trackMapIter != pfParticleToTrackMap.end()) {
             lar_pandora::TrackVector pfParticleTracks = trackMapIter->second;
 
             if (pfParticleTracks.size() > 1)
               std::cerr << "Warning: there was more than one track found for PFParticle with ID " << pfparticlelist[i]->Self() << std::endl;
 
             if (pfParticleTracks.size() > 0) {
               PFParticleData.pfp_trackID[i] = pfParticleTracks.at(0)->ID();
               trackIDtoPFParticleIDMap.insert(std::make_pair(pfParticleTracks.at(0)->ID(), pfparticlelist[i]->Self()));
             }
         }
         else
           std::cerr << "Warning: there was no track found for track-like PFParticle with ID " << pfparticlelist[i]->Self() << std::endl;
       }
       else
         PFParticleData.pfp_isTrack[i] = 0;
 
       if (lar_pandora::LArPandoraHelper::IsShower(pfparticlelist[i])) {
         PFParticleData.pfp_isShower[i] = 1;
         // Set the shower ID.
         auto showerMapIter = pfParticleToShowerMap.find(pfparticlelist[i]);
         if (showerMapIter != pfParticleToShowerMap.end()) {
           lar_pandora::ShowerVector pfParticleShowers = showerMapIter->second;
 
           if (pfParticleShowers.size() > 1)
             std::cerr << "Warning: there was more than one shower found for PFParticle with ID " << pfparticlelist[i]->Self() << std::endl;
 
           if (pfParticleShowers.size() > 0) {
             PFParticleData.pfp_showerID[i] = pfParticleShowers.at(0)->ID();
             showerIDtoPFParticleIDMap.insert(std::make_pair(pfParticleShowers.at(0)->ID(), pfparticlelist[i]->Self()));
           }
         }
         else
           std::cerr << "Warning: there was no shower found for shower-like PFParticle with ID " << pfparticlelist[i]->Self() << std::endl;
       }
       else
         PFParticleData.pfp_isShower[i] = 0;
 
       // Set the cluster IDs.
       auto clusterMapIter = pfParticleToClusterMap.find(pfparticlelist[i]);
       if (clusterMapIter != pfParticleToClusterMap.end()) {
           lar_pandora::ClusterVector pfParticleClusters = clusterMapIter->second;
           PFParticleData.pfp_numClusters[i] = pfParticleClusters.size();
 
           for (size_t j = 0; j < pfParticleClusters.size(); ++j)
             PFParticleData.pfp_clusterIDs[i][j] = pfParticleClusters[j]->ID();
       }
       //else
       //  std::cerr << "Warning: there were no clusters found for PFParticle with ID " << pfparticlelist[i]->Self() << std::endl;
     }
   } // if fSavePFParticleInfo
 
   if (fSaveShowerInfo){
 
     // fill data from all the shower algorithms
     for (size_t iShowerAlgo = 0; iShowerAlgo < NShowerAlgos; ++iShowerAlgo) {
       AnalysisTreeDataStruct::ShowerDataStruct& ShowerData
         = fData->GetShowerData(iShowerAlgo);
       std::vector<recob::Shower> const* pShowers = showerList[iShowerAlgo];
       art::Handle< std::vector<recob::Shower> > showerHandle = showerListHandle[iShowerAlgo];
 
       if (pShowers){
         FillShowers(ShowerData, *pShowers, fSavePFParticleInfo, showerIDtoPFParticleIDMap);
 
         if(fMVAPIDShowerModuleLabel[iShowerAlgo].size()){
           art::FindOneP<anab::MVAPIDResult> fmvapid(showerHandle, evt, fMVAPIDShowerModuleLabel[iShowerAlgo]);
           if(fmvapid.isValid()){
             for(unsigned int iShower=0;iShower<showerHandle->size();++iShower){
               const art::Ptr<anab::MVAPIDResult> pid = fmvapid.at(iShower);
               ShowerData.shwr_pidmvamu[iShower] = pid->mvaOutput.at("muon");
               ShowerData.shwr_pidmvae[iShower] = pid->mvaOutput.at("electron");
               ShowerData.shwr_pidmvapich[iShower] = pid->mvaOutput.at("pich");
               ShowerData.shwr_pidmvaphoton[iShower] = pid->mvaOutput.at("photon");
               ShowerData.shwr_pidmvapr[iShower] = pid->mvaOutput.at("proton");
             }
           } // fmvapid.isValid()
         }
       }
       else ShowerData.MarkMissing(fTree); // tree should reflect lack of data
     } // for iShowerAlgo
 
   } // if fSaveShowerInfo
 
   //track information for multiple trackers
   if (fSaveTrackInfo) {
     for (unsigned int iTracker=0; iTracker < NTrackers; ++iTracker){
       AnalysisTreeDataStruct::TrackDataStruct& TrackerData = fData->GetTrackerData(iTracker);
 
       size_t NTracks = tracklist[iTracker].size();
       // allocate enough space for this number of tracks (but at least for one of them!)
       TrackerData.SetMaxTracks(std::max(NTracks, (size_t) 1));
       TrackerData.Clear(); // clear all the data
 
       TrackerData.ntracks = (int) NTracks;
 
       // now set the tree addresses to the newly allocated memory;
       // this creates the tree branches in case they are not there yet
       SetTrackerAddresses(iTracker);
       if (NTracks > TrackerData.GetMaxTracks()) {
         // got this error? it might be a bug,
         // since we are supposed to have allocated enough space to fit all tracks
         mf::LogError("AnalysisTree:limits") << "event has " << NTracks
                                             << " " << fTrackModuleLabel[iTracker] << " tracks, only "
                                             << TrackerData.GetMaxTracks() << " stored in tree";
       }
 
       //call the track momentum algorithm that gives you momentum based on track range
       // - Should the minimal track length be 50 cm?  The default of 100 has been used.
       trkf::TrackMomentumCalculator trkm{/*100.*/};
 
       for(size_t iTrk=0; iTrk < NTracks; ++iTrk){//loop over tracks
 
         //save t0 from reconstructed flash track matching for every track
         art::FindManyP<anab::T0> fmt0(trackListHandle[iTracker],evt,fFlashT0FinderLabel[iTracker]);
         if (fmt0.isValid()){
           if(fmt0.at(iTrk).size()>0){
             if(fmt0.at(iTrk).size()>1)
               std::cerr << "\n Warning : more than one cosmic tag per track in module! assigning the first tag to the track" << fFlashT0FinderLabel[iTracker];
             TrackerData.trkflashT0[iTrk] = fmt0.at(iTrk).at(0)->Time();
           }
         }
 
         //save t0 from reconstructed flash track matching for every track
         art::FindManyP<anab::T0> fmmct0(trackListHandle[iTracker],evt,fMCT0FinderLabel[iTracker]);
         if (fmmct0.isValid()){
           if(fmmct0.at(iTrk).size()>0){
             if(fmmct0.at(iTrk).size()>1)
               std::cerr << "\n Warning : more than one cosmic tag per track in module! assigning the first tag to the cluster" << fMCT0FinderLabel[iTracker];
             TrackerData.trktrueT0[iTrk] = fmmct0.at(iTrk).at(0)->Time();
           }
         }
 
         //Cosmic Tagger information
         art::FindManyP<anab::CosmicTag> fmct(trackListHandle[iTracker],evt,fCosmicTaggerAssocLabel[iTracker]);
         if (fmct.isValid()){
           TrackerData.trkncosmictags_tagger[iTrk]     = fmct.at(iTrk).size();
           if (fmct.at(iTrk).size()>0){
             if(fmct.at(iTrk).size()>1)
               std::cerr << "\n Warning : more than one cosmic tag per track in module! assigning the first tag to the track" << fCosmicTaggerAssocLabel[iTracker];
             TrackerData.trkcosmicscore_tagger[iTrk] = fmct.at(iTrk).at(0)->CosmicScore();
             TrackerData.trkcosmictype_tagger[iTrk] = fmct.at(iTrk).at(0)->CosmicType();
           }
         }
 
         //Containment Tagger information
         art::FindManyP<anab::CosmicTag> fmcnt(trackListHandle[iTracker],evt,fContainmentTaggerAssocLabel[iTracker]);
         if (fmcnt.isValid()){
           TrackerData.trkncosmictags_containmenttagger[iTrk]     = fmcnt.at(iTrk).size();
           if (fmcnt.at(iTrk).size()>0){
             if(fmcnt.at(iTrk).size()>1)
               std::cerr << "\n Warning : more than one containment tag per track in module! assigning the first tag to the track" << fContainmentTaggerAssocLabel[iTracker];
             TrackerData.trkcosmicscore_containmenttagger[iTrk] = fmcnt.at(iTrk).at(0)->CosmicScore();
             TrackerData.trkcosmictype_containmenttagger[iTrk] = fmcnt.at(iTrk).at(0)->CosmicType();
           }
         }
 
         //Flash match compatibility information
         //Unlike CosmicTagger, Flash match doesn't assign a cosmic tag for every track. For those tracks, AnalysisTree initializes them with -9999 or -99999
         art::FindManyP<anab::CosmicTag> fmbfm(trackListHandle[iTracker],evt,fFlashMatchAssocLabel[iTracker]);
         if (fmbfm.isValid()){
           TrackerData.trkncosmictags_flashmatch[iTrk] = fmbfm.at(iTrk).size();
           if (fmbfm.at(iTrk).size()>0){
             if(fmbfm.at(iTrk).size()>1)
               std::cerr << "\n Warning : more than one cosmic tag per track in module! assigning the first tag to the track" << fFlashMatchAssocLabel[iTracker];
             TrackerData.trkcosmicscore_flashmatch[iTrk] = fmbfm.at(iTrk).at(0)->CosmicScore();
             TrackerData.trkcosmictype_flashmatch[iTrk] = fmbfm.at(iTrk).at(0)->CosmicType();
             //std::cout<<"\n"<<evt.event()<<"\t"<<iTrk<<"\t"<<fmbfm.at(iTrk).at(0)->CosmicScore()<<"\t"<<fmbfm.at(iTrk).at(0)->CosmicType();
           }
         }
 
         art::Ptr<recob::Track> ptrack(trackListHandle[iTracker], iTrk);
         const recob::Track& track = *ptrack;
 
         TVector3 pos, dir_start, dir_end, end;
 
         double tlen = 0., mom = 0.;
         int TrackID = -1;
 
         int ntraj = track.NumberTrajectoryPoints();
         if (ntraj > 0) {
           pos       = track.Vertex<TVector3>();
           dir_start = track.VertexDirection<TVector3>();
           dir_end   = track.EndDirection<TVector3>();
           end       = track.End<TVector3>();
 
           tlen        = track.Length();
           if(track.NumberTrajectoryPoints() > 0)
             mom = track.VertexMomentum();
           // fill non-bezier-track reco branches
           TrackID = track.ID();
 
           double theta_xz = std::atan2(dir_start.X(), dir_start.Z());
           double theta_yz = std::atan2(dir_start.Y(), dir_start.Z());
           double dpos = bdist(pos);  // FIXME - Passing an uncorrected position....
           double dend = bdist(end);  // FIXME - Passing an uncorrected position....
 
           TrackerData.trkId[iTrk]                 = TrackID;
           TrackerData.trkstartx[iTrk]             = pos.X();
           TrackerData.trkstarty[iTrk]             = pos.Y();
           TrackerData.trkstartz[iTrk]             = pos.Z();
           TrackerData.trkstartd[iTrk]             = dpos;
           TrackerData.trkendx[iTrk]               = end.X();
           TrackerData.trkendy[iTrk]               = end.Y();
           TrackerData.trkendz[iTrk]               = end.Z();
           TrackerData.trkendd[iTrk]               = dend;
           TrackerData.trktheta[iTrk]              = dir_start.Theta();
           TrackerData.trkphi[iTrk]                = dir_start.Phi();
           TrackerData.trkstartdcosx[iTrk]         = dir_start.X();
           TrackerData.trkstartdcosy[iTrk]         = dir_start.Y();
           TrackerData.trkstartdcosz[iTrk]         = dir_start.Z();
           TrackerData.trkenddcosx[iTrk]           = dir_end.X();
           TrackerData.trkenddcosy[iTrk]           = dir_end.Y();
           TrackerData.trkenddcosz[iTrk]           = dir_end.Z();
           TrackerData.trkthetaxz[iTrk]            = theta_xz;
           TrackerData.trkthetayz[iTrk]            = theta_yz;
           TrackerData.trkmom[iTrk]                = mom;
           TrackerData.trklen[iTrk]                = tlen;
           TrackerData.trkmomrange[iTrk]           = trkm.GetTrackMomentum(tlen,13);
           //TrackerData.trkmommschi2[iTrk]        = trkm.GetMomentumMultiScatterChi2(ptrack);
           //TrackerData.trkmommsllhd[iTrk]        = trkm.GetMomentumMultiScatterLLHD(ptrack);
 
           if (fSavePFParticleInfo) {
             auto mapIter = trackIDtoPFParticleIDMap.find(TrackID);
             if (mapIter != trackIDtoPFParticleIDMap.end()) {
                 // This track has a corresponding PFParticle.
                 TrackerData.trkhasPFParticle[iTrk] = 1;
                 TrackerData.trkPFParticleID[iTrk] = mapIter->second;
             }
             else
                 TrackerData.trkhasPFParticle[iTrk] = 0;
           }
 
         } // if we have trajectory
 
         // find vertices associated with this track
         /*
           art::FindMany<recob::Vertex> fmvtx(trackListHandle[iTracker], evt, fVertexModuleLabel[iTracker]);
           if(fmvtx.isValid()) {
           std::vector<const recob::Vertex*> verts = fmvtx.at(iTrk);
           // should have two at most
           for(size_t ivx = 0; ivx < verts.size(); ++ivx) {
           verts[ivx]->XYZ(xyz);
           // find the vertex in TrackerData to get the index
           short theVtx = -1;
           for(short jvx = 0; jvx < TrackerData.nvtx; ++jvx) {
           if(TrackerData.vtx[jvx][2] == xyz[2]) {
           theVtx = jvx;
           break;
           }
           } // jvx
           // decide if it should be assigned to the track Start or End.
           // A simple dz test should suffice
           if(fabs(xyz[2] - TrackerData.trkstartz[iTrk]) <
           fabs(xyz[2] - TrackerData.trkendz[iTrk])) {
           TrackerData.trksvtxid[iTrk] = theVtx;
           } else {
           TrackerData.trkevtxid[iTrk] = theVtx;
           }
           } // vertices
           } // fmvtx.isValid()
         */
 
 
         /* //commented out because now have several Vertices
            Float_t minsdist = 10000;
            Float_t minedist = 10000;
            for (int ivx = 0; ivx < NVertices && ivx < kMaxVertices; ++ivx){
            Float_t sdist = sqrt(pow(TrackerData.trkstartx[iTrk]-VertexData.vtxx[ivx],2)+
            pow(TrackerData.trkstarty[iTrk]-VertexData.vtxy[ivx],2)+
            pow(TrackerData.trkstartz[iTrk]-VertexData.vtxz[ivx],2));
            Float_t edist = sqrt(pow(TrackerData.trkendx[iTrk]-VertexData.vtxx[ivx],2)+
            pow(TrackerData.trkendy[iTrk]-VertexData.vtxy[ivx],2)+
            pow(TrackerData.trkendz[iTrk]-VertexData.vtxz[ivx],2));
            if (sdist<minsdist){
            minsdist = sdist;
            if (minsdist<10) TrackerData.trksvtxid[iTrk] = ivx;
            }
            if (edist<minedist){
            minedist = edist;
            if (minedist<10) TrackerData.trkevtxid[iTrk] = ivx;
            }
            }*/
 
         // find particle ID info
         /*
         // Note from Jake Calcutt: There has been a breaking change in the definition
         // of the anab::ParticleID class. If you are using this class and want this info
         // in this tree, these must be updated
         art::FindMany<anab::ParticleID> fmpid(trackListHandle[iTracker], evt, fParticleIDModuleLabel[iTracker]);
         if(fmpid.isValid()) {
           std::vector<const anab::ParticleID*> pids = fmpid.at(iTrk);
           //if(pids.size() > 1) {
           //mf::LogError("AnalysisTree:limits")
           //<< "the " << fTrackModuleLabel[iTracker] << " track #" << iTrk
           //<< " has " << pids.size()
           //<< " set of ParticleID variables. Only one stored in the tree";
           //}
           for (size_t ipid = 0; ipid < pids.size(); ++ipid){
             if (!pids[ipid]->PlaneID().isValid) continue;
             int planenum = pids[ipid]->PlaneID().Plane;
             if (planenum<0||planenum>2) continue;
             TrackerData.trkpidpdg[iTrk][planenum] = pids[ipid]->Pdg();
             TrackerData.trkpidchi[iTrk][planenum] = pids[ipid]->MinChi2();
             TrackerData.trkpidchipr[iTrk][planenum] = pids[ipid]->Chi2Proton();
             TrackerData.trkpidchika[iTrk][planenum] = pids[ipid]->Chi2Kaon();
             TrackerData.trkpidchipi[iTrk][planenum] = pids[ipid]->Chi2Pion();
             TrackerData.trkpidchimu[iTrk][planenum] = pids[ipid]->Chi2Muon();
             TrackerData.trkpidpida[iTrk][planenum] = pids[ipid]->PIDA();
           }
         } // fmpid.isValid()
         */
         if(fMVAPIDTrackModuleLabel[iTracker].size()){
           art::FindOneP<anab::MVAPIDResult> fmvapid(trackListHandle[iTracker], evt, fMVAPIDTrackModuleLabel[iTracker]);
           if(fmvapid.isValid()) {
             const art::Ptr<anab::MVAPIDResult> pid = fmvapid.at(iTrk);
             TrackerData.trkpidmvamu[iTrk] = pid->mvaOutput.at("muon");
             TrackerData.trkpidmvae[iTrk] = pid->mvaOutput.at("electron");
             TrackerData.trkpidmvapich[iTrk] = pid->mvaOutput.at("pich");
             TrackerData.trkpidmvaphoton[iTrk] = pid->mvaOutput.at("photon");
             TrackerData.trkpidmvapr[iTrk] = pid->mvaOutput.at("proton");
           } // fmvapid.isValid()
         }
         art::FindMany<anab::Calorimetry> fmcal(trackListHandle[iTracker], evt, fCalorimetryModuleLabel[iTracker]);
         if (fmcal.isValid()){
           std::vector<const anab::Calorimetry*> calos = fmcal.at(iTrk);
           if (calos.size() > TrackerData.GetMaxPlanesPerTrack(iTrk)) {
             // if you get this message, there is probably a bug somewhere since
             // the calorimetry planes should be 3.
             mf::LogError("AnalysisTree:limits")
               << "the " << fTrackModuleLabel[iTracker] << " track #" << iTrk
               << " has " << calos.size() << " planes for calorimetry , only "
               << TrackerData.GetMaxPlanesPerTrack(iTrk) << " stored in tree";
           }
           for (size_t ical = 0; ical<calos.size(); ++ical){
             if (!calos[ical]) continue;
             if (!calos[ical]->PlaneID().isValid) continue;
             int planenum = calos[ical]->PlaneID().Plane;
             if (planenum<0||planenum>2) continue;
             TrackerData.trkke[iTrk][planenum]    = calos[ical]->KineticEnergy();
             TrackerData.trkrange[iTrk][planenum] = calos[ical]->Range();
             //For now make the second argument as 13 for muons.
             TrackerData.trkpitchc[iTrk][planenum]= calos[ical] -> TrkPitchC();
             const size_t NHits = calos[ical] -> dEdx().size();
             TrackerData.ntrkhits[iTrk][planenum] = (int) NHits;
             if (NHits > TrackerData.GetMaxHitsPerTrack(iTrk, planenum)) {
               // if you get this error, you'll have to increase kMaxTrackHits
               mf::LogError("AnalysisTree:limits")
                 << "the " << fTrackModuleLabel[iTracker] << " track #" << iTrk
                 << " has " << NHits << " hits on calorimetry plane #" << planenum
                 <<", only "
                 << TrackerData.GetMaxHitsPerTrack(iTrk, planenum) << " stored in tree";
             }
             if (!isCosmics){
               for(size_t iTrkHit = 0; iTrkHit < NHits && iTrkHit < TrackerData.GetMaxHitsPerTrack(iTrk, planenum); ++iTrkHit) {
                 TrackerData.trkdedx[iTrk][planenum][iTrkHit]  = (calos[ical] -> dEdx())[iTrkHit];
                 TrackerData.trkdqdx[iTrk][planenum][iTrkHit]  = (calos[ical] -> dQdx())[iTrkHit];
                 TrackerData.trkresrg[iTrk][planenum][iTrkHit] = (calos[ical] -> ResidualRange())[iTrkHit];
                 TrackerData.trktpc[iTrk][planenum][iTrkHit]   = (calos[ical] -> PlaneID()).TPC;
                 const auto& TrkPos = (calos[ical] -> XYZ())[iTrkHit];
                 auto& TrkXYZ = TrackerData.trkxyz[iTrk][planenum][iTrkHit];
                 TrkXYZ[0] = TrkPos.X();
                 TrkXYZ[1] = TrkPos.Y();
                 TrkXYZ[2] = TrkPos.Z();
               } // for track hits
             }
           } // for calorimetry info
           if(TrackerData.ntrkhits[iTrk][0] > TrackerData.ntrkhits[iTrk][1] && TrackerData.ntrkhits[iTrk][0] > TrackerData.ntrkhits[iTrk][2]) TrackerData.trkpidbestplane[iTrk] = 0;
           else if(TrackerData.ntrkhits[iTrk][1] > TrackerData.ntrkhits[iTrk][0] && TrackerData.ntrkhits[iTrk][1] > TrackerData.ntrkhits[iTrk][2]) TrackerData.trkpidbestplane[iTrk] = 1;
           else if(TrackerData.ntrkhits[iTrk][2] > TrackerData.ntrkhits[iTrk][0] && TrackerData.ntrkhits[iTrk][2] > TrackerData.ntrkhits[iTrk][1]) TrackerData.trkpidbestplane[iTrk] = 2;
           else if(TrackerData.ntrkhits[iTrk][2] == TrackerData.ntrkhits[iTrk][0] && TrackerData.ntrkhits[iTrk][2] > TrackerData.ntrkhits[iTrk][1]) TrackerData.trkpidbestplane[iTrk] = 2;
           else if(TrackerData.ntrkhits[iTrk][2] == TrackerData.ntrkhits[iTrk][1] && TrackerData.ntrkhits[iTrk][2] > TrackerData.ntrkhits[iTrk][0]) TrackerData.trkpidbestplane[iTrk] = 2;
           else if(TrackerData.ntrkhits[iTrk][1] == TrackerData.ntrkhits[iTrk][0] && TrackerData.ntrkhits[iTrk][1] > TrackerData.ntrkhits[iTrk][2]) TrackerData.trkpidbestplane[iTrk] = 0;
           else if(TrackerData.ntrkhits[iTrk][1] == TrackerData.ntrkhits[iTrk][0] && TrackerData.ntrkhits[iTrk][1] == TrackerData.ntrkhits[iTrk][2]) TrackerData.trkpidbestplane[iTrk] = 2;
 
           // FIXME - Do i want to add someway to work out the best TPC???....
         } // if has calorimetry info
 
         //track truth information
         if (isMC){
           //get the hits on each plane
           art::FindManyP<recob::Hit>      fmht(trackListHandle[iTracker], evt, fTrackModuleLabel[iTracker]);
           std::vector< art::Ptr<recob::Hit> > allHits = fmht.at(iTrk);
           std::vector< art::Ptr<recob::Hit> > hits[kNplanes];
 
           for(size_t ah = 0; ah < allHits.size(); ++ah){
             if (/* allHits[ah]->WireID().Plane >= 0 && */ // always true
                 allHits[ah]->WireID().Plane <  3){
               hits[allHits[ah]->WireID().Plane].push_back(allHits[ah]);
             }
           }
           for (size_t ipl = 0; ipl < 3; ++ipl){
             double maxe = 0;
             HitsPurity(clockData, hits[ipl],TrackerData.trkidtruth[iTrk][ipl],TrackerData.trkpurtruth[iTrk][ipl],maxe);
             //std::cout<<"\n"<<iTracker<<"\t"<<iTrk<<"\t"<<ipl<<"\t"<<trkidtruth[iTracker][iTrk][ipl]<<"\t"<<trkpurtruth[iTracker][iTrk][ipl]<<"\t"<<maxe;
             if (TrackerData.trkidtruth[iTrk][ipl]>0){
               const art::Ptr<simb::MCTruth> mc = pi_serv->TrackIdToMCTruth_P(TrackerData.trkidtruth[iTrk][ipl]);
               TrackerData.trkorigin[iTrk][ipl] = mc->Origin();
               const simb::MCParticle *particle = pi_serv->TrackIdToParticle_P(TrackerData.trkidtruth[iTrk][ipl]);
               double tote = 0;
               const std::vector<const sim::IDE*> vide=bt_serv->TrackIdToSimIDEs_Ps(TrackerData.trkidtruth[iTrk][ipl]);
               for (auto ide: vide) {
                 tote += ide->energy;
               }
               TrackerData.trkpdgtruth[iTrk][ipl] = particle->PdgCode();
               TrackerData.trkefftruth[iTrk][ipl] = maxe/(tote/kNplanes); //tote include both induction and collection energies
               //std::cout<<"\n"<<trkpdgtruth[iTracker][iTrk][ipl]<<"\t"<<trkefftruth[iTracker][iTrk][ipl];
             }
           }
 
           double maxe = 0;
           HitsPurity(clockData, allHits,TrackerData.trkg4id[iTrk],TrackerData.trkpurity[iTrk],maxe);
           if (TrackerData.trkg4id[iTrk]>0){
             const art::Ptr<simb::MCTruth> mc = pi_serv->TrackIdToMCTruth_P(TrackerData.trkg4id[iTrk]);
             TrackerData.trkorig[iTrk] = mc->Origin();
           }
           if (allHits.size()){
             std::vector<art::Ptr<recob::Hit> > all_hits;
             art::Handle<std::vector<recob::Hit> > hithandle;
             float totenergy = 0;
             if (evt.get(allHits[0].id(), hithandle)){
               art::fill_ptr_vector(all_hits, hithandle);
               for(size_t h = 0; h < all_hits.size(); ++h){
 
                 art::Ptr<recob::Hit> hit = all_hits[h];
                 std::vector<sim::IDE*> ides;
                 //bt_serv->HitToSimIDEs(hit,ides);
                 std::vector<sim::TrackIDE> eveIDs = bt_serv->HitToEveTrackIDEs(clockData, hit);
 
                 for(size_t e = 0; e < eveIDs.size(); ++e){
                   //std::cout<<h<<" "<<e<<" "<<eveIDs[e].trackID<<" "<<eveIDs[e].energy<<" "<<eveIDs[e].energyFrac<<std::endl;
                   if (eveIDs[e].trackID==TrackerData.trkg4id[iTrk]) totenergy += eveIDs[e].energy;
                 }
               }
             }
             if (totenergy) TrackerData.trkcompleteness[iTrk] = maxe/totenergy;
           }
         }//end if (isMC)
       }//end loop over track
     }//end loop over track module labels
   }// end (fSaveTrackInfo)
 
   /*trkf::TrackMomentumCalculator trkm;
     std::cout<<"\t"<<trkm.GetTrackMomentum(200,2212)<<"\t"<<trkm.GetTrackMomentum(-10, 13)<<"\t"<<trkm.GetTrackMomentum(300,-19)<<"\n";
   */
 
   //Save Vertex information for multiple algorithms
   if (fSaveVertexInfo){
     for (unsigned int iVertexAlg=0; iVertexAlg < NVertexAlgos; ++iVertexAlg){
       AnalysisTreeDataStruct::VertexDataStruct& VertexData = fData->GetVertexData(iVertexAlg);
 
       size_t NVertices = vertexlist[iVertexAlg].size();
 
       VertexData.SetMaxVertices(std::max(NVertices, (size_t) 1));
       VertexData.Clear(); // clear all the data
 
       VertexData.nvtx = (short) NVertices;
 
       // now set the tree addresses to the newly allocated memory;
       // this creates the tree branches in case they are not there yet
       SetVertexAddresses(iVertexAlg);
       if (NVertices > VertexData.GetMaxVertices()) {
         // got this error? it might be a bug,
         // since we are supposed to have allocated enough space to fit all tracks
         mf::LogError("AnalysisTree:limits") << "event has " << NVertices
                                             << " " << fVertexModuleLabel[iVertexAlg] << " tracks, only "
                                             << VertexData.GetMaxVertices() << " stored in tree";
       }
 
       for (size_t i = 0; i < NVertices && i < kMaxVertices ; ++i){//loop over hits
             VertexData.vtxId[i] = vertexlist[iVertexAlg][i]->ID();
             Double_t xyz[3] = {};
             vertexlist[iVertexAlg][i] -> XYZ(xyz);
         VertexData.vtxx[i] = xyz[0];
         VertexData.vtxy[i] = xyz[1];
         VertexData.vtxz[i] = xyz[2];
 
         if (fSavePFParticleInfo) {
           auto mapIter = vertexIDtoPFParticleIDMap.find(vertexlist[iVertexAlg][i]->ID());
           if (mapIter != vertexIDtoPFParticleIDMap.end()) {
             // This vertex has a corresponding PFParticle.
             VertexData.vtxhasPFParticle[i] = 1;
             VertexData.vtxPFParticleID[i] = mapIter->second;
           }
           else
             VertexData.vtxhasPFParticle[i] = 0;
         }
 
         // find PFParticle ID info
         art::FindMany<recob::PFParticle> fmPFParticle(vertexListHandle[iVertexAlg], evt, fPFParticleModuleLabel);
         if(fmPFParticle.isValid()) {
           std::vector<const recob::PFParticle*> pfparticles = fmPFParticle.at(i);
           if(pfparticles.size() > 1)
           std::cerr << "Warning: more than one associated PFParticle found for a vertex. Only one stored in tree." << std::endl;
           if (pfparticles.size() == 0)
           VertexData.vtxhasPFParticle[i] = 0;
           else {
             VertexData.vtxhasPFParticle[i] = 1;
             VertexData.vtxPFParticleID[i] = pfparticles.at(0)->Self();
           }
         } // fmPFParticle.isValid()
       }
     }
   }
 
   //mc truth information
   if (isMC){
 
     // Find the simb::MCFlux objects corresponding to
     // each simb::MCTruth object made by the generator with
     // the label fGenieGenModuleLabel
     art::FindOne<simb::MCFlux> find_mcflux(mctruthListHandle,
                                            evt, fGenieGenModuleLabel);
 
     if (fSaveCryInfo){
       //store cry (cosmic generator information)
       fData->mcevts_truthcry = mclistcry.size();
       fData->cry_no_primaries = nCryPrimaries;
       //fData->cry_no_primaries;
       for(Int_t iPartc = 0; iPartc < mctruthcry->NParticles(); ++iPartc){
         const simb::MCParticle& partc(mctruthcry->GetParticle(iPartc));
         fData->cry_primaries_pdg[iPartc]=partc.PdgCode();
         fData->cry_Eng[iPartc]=partc.E();
         fData->cry_Px[iPartc]=partc.Px();
         fData->cry_Py[iPartc]=partc.Py();
         fData->cry_Pz[iPartc]=partc.Pz();
         fData->cry_P[iPartc]=partc.P();
         fData->cry_StartPointx[iPartc] = partc.Vx();
         fData->cry_StartPointy[iPartc] = partc.Vy();
         fData->cry_StartPointz[iPartc] = partc.Vz();
         fData->cry_StartPointt[iPartc] = partc.T();
         fData->cry_status_code[iPartc]=partc.StatusCode();
         fData->cry_mass[iPartc]=partc.Mass();
         fData->cry_trackID[iPartc]=partc.TrackId();
         fData->cry_ND[iPartc]=partc.NumberDaughters();
         fData->cry_mother[iPartc]=partc.Mother();
       } // for cry particles
     }// end fSaveCryInfo
 
     // Save the protoDUNE beam generator information
     if(fSaveProtoInfo){
       fData->proto_no_primaries = nProtoPrimaries;
       for(Int_t iPartp = 0; iPartp < nProtoPrimaries; ++iPartp){
         const simb::MCParticle& partp(mctruthproto->GetParticle(iPartp));
 
         fData->proto_isGoodParticle[iPartp] = (partp.Process() == "primary");
         fData->proto_vx[iPartp] = partp.Vx();
         fData->proto_vy[iPartp] = partp.Vy();
         fData->proto_vz[iPartp] = partp.Vz();
         fData->proto_t[iPartp] = partp.T();
         fData->proto_px[iPartp] = partp.Px();
         fData->proto_py[iPartp] = partp.Py();
         fData->proto_pz[iPartp] = partp.Pz();
         fData->proto_momentum[iPartp] = partp.P();
         fData->proto_energy[iPartp] = partp.E();
         fData->proto_pdg[iPartp] = partp.PdgCode();
         // We will deal with the matching to GEANT later
       }
     }
 
     //save neutrino interaction information
     fData->mcevts_truth = mclist.size();
     if (fData->mcevts_truth > 0){//at least one mc record
       if (fSaveGenieInfo){
         int neutrino_i = 0;
         for(unsigned int iList = 0; (iList < mclist.size()) && (neutrino_i < kMaxTruth) ; ++iList){
           if (mclist[iList]->NeutrinoSet()){
             fData->nuPDG_truth[neutrino_i]  = mclist[iList]->GetNeutrino().Nu().PdgCode();
             fData->ccnc_truth[neutrino_i]   = mclist[iList]->GetNeutrino().CCNC();
             fData->mode_truth[neutrino_i]   = mclist[iList]->GetNeutrino().Mode();
             fData->Q2_truth[neutrino_i]     = mclist[iList]->GetNeutrino().QSqr();
             fData->W_truth[neutrino_i]      = mclist[iList]->GetNeutrino().W();
             fData->X_truth[neutrino_i]      = mclist[iList]->GetNeutrino().X();
             fData->Y_truth[neutrino_i]      = mclist[iList]->GetNeutrino().Y();
             fData->hitnuc_truth[neutrino_i] = mclist[iList]->GetNeutrino().HitNuc();
             fData->enu_truth[neutrino_i]    = mclist[iList]->GetNeutrino().Nu().E();
             fData->nuvtxx_truth[neutrino_i] = mclist[iList]->GetNeutrino().Nu().Vx();
             fData->nuvtxy_truth[neutrino_i] = mclist[iList]->GetNeutrino().Nu().Vy();
             fData->nuvtxz_truth[neutrino_i] = mclist[iList]->GetNeutrino().Nu().Vz();
             if (mclist[iList]->GetNeutrino().Nu().P()){
               fData->nu_dcosx_truth[neutrino_i] = mclist[iList]->GetNeutrino().Nu().Px()/mclist[iList]->GetNeutrino().Nu().P();
               fData->nu_dcosy_truth[neutrino_i] = mclist[iList]->GetNeutrino().Nu().Py()/mclist[iList]->GetNeutrino().Nu().P();
               fData->nu_dcosz_truth[neutrino_i] = mclist[iList]->GetNeutrino().Nu().Pz()/mclist[iList]->GetNeutrino().Nu().P();
             }
             fData->lep_mom_truth[neutrino_i] = mclist[iList]->GetNeutrino().Lepton().P();
             if (mclist[iList]->GetNeutrino().Lepton().P()){
               fData->lep_dcosx_truth[neutrino_i] = mclist[iList]->GetNeutrino().Lepton().Px()/mclist[iList]->GetNeutrino().Lepton().P();
               fData->lep_dcosy_truth[neutrino_i] = mclist[iList]->GetNeutrino().Lepton().Py()/mclist[iList]->GetNeutrino().Lepton().P();
               fData->lep_dcosz_truth[neutrino_i] = mclist[iList]->GetNeutrino().Lepton().Pz()/mclist[iList]->GetNeutrino().Lepton().P();
             }
 
             //flux information
             //
             // Double-check that a simb::MCFlux object is associated with the
             // current simb::MCTruth object. For GENIE events, these should
             // always accompany each other. Other generators (e.g., MARLEY) may
             // create simb::MCTruth objects without corresponding simb::MCFlux
             // objects. -- S. Gardiner
             if (find_mcflux.isValid()) {
               auto flux_maybe_ref = find_mcflux.at(iList);
               if (flux_maybe_ref.isValid()) {
                 auto flux_ref = flux_maybe_ref.ref();
                 fData->vx_flux[neutrino_i]        = flux_ref.fvx;
                 fData->vy_flux[neutrino_i]        = flux_ref.fvy;
                 fData->vz_flux[neutrino_i]        = flux_ref.fvz;
                 fData->pdpx_flux[neutrino_i]      = flux_ref.fpdpx;
                 fData->pdpy_flux[neutrino_i]      = flux_ref.fpdpy;
                 fData->pdpz_flux[neutrino_i]      = flux_ref.fpdpz;
                 fData->ppdxdz_flux[neutrino_i]    = flux_ref.fppdxdz;
                 fData->ppdydz_flux[neutrino_i]    = flux_ref.fppdydz;
                 fData->pppz_flux[neutrino_i]      = flux_ref.fpppz;
 
                 fData->ptype_flux[neutrino_i]      = flux_ref.fptype;
                 fData->ppvx_flux[neutrino_i]       = flux_ref.fppvx;
                 fData->ppvy_flux[neutrino_i]       = flux_ref.fppvy;
                 fData->ppvz_flux[neutrino_i]       = flux_ref.fppvz;
                 fData->muparpx_flux[neutrino_i]    = flux_ref.fmuparpx;
                 fData->muparpy_flux[neutrino_i]    = flux_ref.fmuparpy;
                 fData->muparpz_flux[neutrino_i]    = flux_ref.fmuparpz;
                 fData->mupare_flux[neutrino_i]     = flux_ref.fmupare;
 
                 fData->tgen_flux[neutrino_i]     = flux_ref.ftgen;
                 fData->tgptype_flux[neutrino_i]  = flux_ref.ftgptype;
                 fData->tgppx_flux[neutrino_i]    = flux_ref.ftgppx;
                 fData->tgppy_flux[neutrino_i]    = flux_ref.ftgppy;
                 fData->tgppz_flux[neutrino_i]    = flux_ref.ftgppz;
                 fData->tprivx_flux[neutrino_i]   = flux_ref.ftprivx;
                 fData->tprivy_flux[neutrino_i]   = flux_ref.ftprivy;
                 fData->tprivz_flux[neutrino_i]   = flux_ref.ftprivz;
 
                 fData->dk2gen_flux[neutrino_i]   = flux_ref.fdk2gen;
                 fData->gen2vtx_flux[neutrino_i]   = flux_ref.fgen2vtx;
 
                 fData->tpx_flux[neutrino_i]    = flux_ref.ftpx;
                 fData->tpy_flux[neutrino_i]    = flux_ref.ftpy;
                 fData->tpz_flux[neutrino_i]    = flux_ref.ftpz;
                 fData->tptype_flux[neutrino_i] = flux_ref.ftptype;
               } // flux_maybe_ref.isValid()
             } // find_mcflux.isValid()
             neutrino_i++;
           }//mclist is NeutrinoSet()
         }//loop over mclist
 
         if (mctruth->NeutrinoSet()){
           //genie particles information
           fData->genie_no_primaries = mctruth->NParticles();
 
           size_t StoreParticles = std::min((size_t) fData->genie_no_primaries, fData->GetMaxGeniePrimaries());
           if (fData->genie_no_primaries > (int) StoreParticles) {
             // got this error? it might be a bug,
             // since the structure should have enough room for everything
             mf::LogError("AnalysisTree:limits") << "event has "
                                                 << fData->genie_no_primaries << " MC particles, only "
                                                 << StoreParticles << " stored in tree";
           }
           for(size_t iPart = 0; iPart < StoreParticles; ++iPart){
             const simb::MCParticle& part(mctruth->GetParticle(iPart));
             fData->genie_primaries_pdg[iPart]=part.PdgCode();
             fData->genie_Eng[iPart]=part.E();
             fData->genie_Px[iPart]=part.Px();
             fData->genie_Py[iPart]=part.Py();
             fData->genie_Pz[iPart]=part.Pz();
             fData->genie_P[iPart]=part.P();
             fData->genie_status_code[iPart]=part.StatusCode();
             fData->genie_mass[iPart]=part.Mass();
             fData->genie_trackID[iPart]=part.TrackId();
             fData->genie_ND[iPart]=part.NumberDaughters();
             fData->genie_mother[iPart]=part.Mother();
           } // for particle
           //const simb::MCNeutrino& nu(mctruth->GetNeutrino());
         } //if neutrino set
       }// end (fSaveGenieInfo)
 
       //Extract MC Shower information and fill the Shower branches
       if (fSaveMCShowerInfo){
         fData->no_mcshowers = nMCShowers;
         size_t shwr = 0;
         for(std::vector<sim::MCShower>::const_iterator imcshwr = mcshowerh->begin();
             imcshwr != mcshowerh->end(); ++imcshwr) {
           const sim::MCShower& mcshwr = *imcshwr;
           fData->mcshwr_origin[shwr]          = mcshwr.Origin();
           fData->mcshwr_pdg[shwr]             = mcshwr.PdgCode();
           fData->mcshwr_TrackId[shwr]         = mcshwr.TrackID();
           fData->mcshwr_Process[shwr]         = mcshwr.Process();
           fData->mcshwr_startX[shwr]          = mcshwr.Start().X();
           fData->mcshwr_startY[shwr]          = mcshwr.Start().Y();
           fData->mcshwr_startZ[shwr]          = mcshwr.Start().Z();
           fData->mcshwr_endX[shwr]            = mcshwr.End().X();
           fData->mcshwr_endY[shwr]            = mcshwr.End().Y();
           fData->mcshwr_endZ[shwr]            = mcshwr.End().Z();
           if (mcshwr.DetProfile().E()!= 0){
             fData->mcshwr_isEngDeposited[shwr] = 1;
             fData->mcshwr_CombEngX[shwr]        = mcshwr.DetProfile().X();
             fData->mcshwr_CombEngY[shwr]        = mcshwr.DetProfile().Y();
             fData->mcshwr_CombEngZ[shwr]        = mcshwr.DetProfile().Z();
             fData->mcshwr_CombEngPx[shwr]       = mcshwr.DetProfile().Px();
             fData->mcshwr_CombEngPy[shwr]       = mcshwr.DetProfile().Py();
             fData->mcshwr_CombEngPz[shwr]       = mcshwr.DetProfile().Pz();
             fData->mcshwr_CombEngE[shwr]        = mcshwr.DetProfile().E();
             fData->mcshwr_dEdx[shwr]            = mcshwr.dEdx();
             fData->mcshwr_StartDirX[shwr]       = mcshwr.StartDir().X();
             fData->mcshwr_StartDirY[shwr]       = mcshwr.StartDir().Y();
             fData->mcshwr_StartDirZ[shwr]       = mcshwr.StartDir().Z();
           }
           else
             fData->mcshwr_isEngDeposited[shwr] = 0;
           fData->mcshwr_Motherpdg[shwr]       = mcshwr.MotherPdgCode();
           fData->mcshwr_MotherTrkId[shwr]     = mcshwr.MotherTrackID();
           fData->mcshwr_MotherProcess[shwr]   = mcshwr.MotherProcess();
           fData->mcshwr_MotherstartX[shwr]    = mcshwr.MotherStart().X();
           fData->mcshwr_MotherstartY[shwr]    = mcshwr.MotherStart().Y();
           fData->mcshwr_MotherstartZ[shwr]    = mcshwr.MotherStart().Z();
           fData->mcshwr_MotherendX[shwr]      = mcshwr.MotherEnd().X();
           fData->mcshwr_MotherendY[shwr]      = mcshwr.MotherEnd().Y();
           fData->mcshwr_MotherendZ[shwr]      = mcshwr.MotherEnd().Z();
           fData->mcshwr_Ancestorpdg[shwr]     = mcshwr.AncestorPdgCode();
           fData->mcshwr_AncestorTrkId[shwr]   = mcshwr.AncestorTrackID();
           fData->mcshwr_AncestorProcess[shwr] = mcshwr.AncestorProcess();
           fData->mcshwr_AncestorstartX[shwr]  = mcshwr.AncestorStart().X();
           fData->mcshwr_AncestorstartY[shwr]  = mcshwr.AncestorStart().Y();
           fData->mcshwr_AncestorstartZ[shwr]  = mcshwr.AncestorStart().Z();
           fData->mcshwr_AncestorendX[shwr]    = mcshwr.AncestorEnd().X();
           fData->mcshwr_AncestorendY[shwr]    = mcshwr.AncestorEnd().Y();
           fData->mcshwr_AncestorendZ[shwr]    = mcshwr.AncestorEnd().Z();
           ++shwr;
         }
         fData->mcshwr_Process.resize(shwr);
         fData->mcshwr_MotherProcess.resize(shwr);
         fData->mcshwr_AncestorProcess.resize(shwr);
       }//End if (fSaveMCShowerInfo){
 
       //Extract MC Track information and fill the Shower branches
       if (fSaveMCTrackInfo){
         fData->no_mctracks = nMCTracks;
         size_t trk = 0;
         for(std::vector<sim::MCTrack>::const_iterator imctrk = mctrackh->begin();imctrk != mctrackh->end(); ++imctrk) {
           const sim::MCTrack& mctrk = *imctrk;
           TLorentzVector tpcstart, tpcend, tpcmom;
           double plen = driftedLength(detProp, mctrk, tpcstart, tpcend, tpcmom);
           fData->mctrk_origin[trk]          = mctrk.Origin();
           fData->mctrk_pdg[trk]             = mctrk.PdgCode();
           fData->mctrk_TrackId[trk]         = mctrk.TrackID();
           fData->mctrk_Process[trk]         = mctrk.Process();
           fData->mctrk_startX[trk]          = mctrk.Start().X();
           fData->mctrk_startY[trk]          = mctrk.Start().Y();
           fData->mctrk_startZ[trk]          = mctrk.Start().Z();
           fData->mctrk_endX[trk]            = mctrk.End().X();
           fData->mctrk_endY[trk]            = mctrk.End().Y();
           fData->mctrk_endZ[trk]            = mctrk.End().Z();
           fData->mctrk_Motherpdg[trk]       = mctrk.MotherPdgCode();
           fData->mctrk_MotherTrkId[trk]     = mctrk.MotherTrackID();
           fData->mctrk_MotherProcess[trk]   = mctrk.MotherProcess();
           fData->mctrk_MotherstartX[trk]    = mctrk.MotherStart().X();
           fData->mctrk_MotherstartY[trk]    = mctrk.MotherStart().Y();
           fData->mctrk_MotherstartZ[trk]    = mctrk.MotherStart().Z();
           fData->mctrk_MotherendX[trk]      = mctrk.MotherEnd().X();
           fData->mctrk_MotherendY[trk]      = mctrk.MotherEnd().Y();
           fData->mctrk_MotherendZ[trk]      = mctrk.MotherEnd().Z();
           fData->mctrk_Ancestorpdg[trk]     = mctrk.AncestorPdgCode();
           fData->mctrk_AncestorTrkId[trk]   = mctrk.AncestorTrackID();
           fData->mctrk_AncestorProcess[trk] = mctrk.AncestorProcess();
           fData->mctrk_AncestorstartX[trk]  = mctrk.AncestorStart().X();
           fData->mctrk_AncestorstartY[trk]  = mctrk.AncestorStart().Y();
           fData->mctrk_AncestorstartZ[trk]  = mctrk.AncestorStart().Z();
           fData->mctrk_AncestorendX[trk]    = mctrk.AncestorEnd().X();
           fData->mctrk_AncestorendY[trk]    = mctrk.AncestorEnd().Y();
           fData->mctrk_AncestorendZ[trk]    = mctrk.AncestorEnd().Z();
 
           fData->mctrk_len_drifted[trk]       = plen;
 
           if (plen != 0){
             fData->mctrk_startX_drifted[trk] = tpcstart.X();
             fData->mctrk_startY_drifted[trk] = tpcstart.Y();
             fData->mctrk_startZ_drifted[trk] = tpcstart.Z();
             fData->mctrk_endX_drifted[trk]   = tpcend.X();
             fData->mctrk_endY_drifted[trk]   = tpcend.Y();
             fData->mctrk_endZ_drifted[trk]   = tpcend.Z();
             fData->mctrk_p_drifted[trk]      = tpcmom.Vect().Mag();
             fData->mctrk_px_drifted[trk]     = tpcmom.X();
             fData->mctrk_py_drifted[trk]     = tpcmom.Y();
             fData->mctrk_pz_drifted[trk]     = tpcmom.Z();
           }
           ++trk;
         }
 
         fData->mctrk_Process.resize(trk);
         fData->mctrk_MotherProcess.resize(trk);
         fData->mctrk_AncestorProcess.resize(trk);
       }//End if (fSaveMCTrackInfo){
 
 
       //GEANT particles information
       if (fSaveGeantInfo){
 
         const sim::ParticleList& plist = pi_serv->ParticleList();
 
         std::string pri("primary");
         int primary=0;
         int active = 0;
         size_t geant_particle=0;
         sim::ParticleList::const_iterator itPart = plist.begin(),
           pend = plist.end(); // iterator to pairs (track id, particle)
 
         // helper map track ID => index
         std::map<int, size_t> TrackIDtoIndex;
         std::vector<int> gpdg;
         std::vector<int> gmother;
         for(size_t iPart = 0; (iPart < plist.size()) && (itPart != pend); ++iPart){
           const simb::MCParticle* pPart = (itPart++)->second;
           if (!pPart) {
             throw art::Exception(art::errors::LogicError)
               << "GEANT particle #" << iPart << " returned a null pointer";
           }
 
           //++geant_particle;
           bool isPrimary = pPart->Process() == pri;
           int TrackID = pPart->TrackId();
           TrackIDtoIndex.emplace(TrackID, iPart);
           gpdg.push_back(pPart->PdgCode());
           gmother.push_back(pPart->Mother());
           if (iPart < fData->GetMaxGEANTparticles()) {
             if (pPart->E()<fG4minE&&(!isPrimary)) continue;
             if (isPrimary) ++primary;
 
             TLorentzVector mcstart, mcend, mcstartdrifted, mcenddrifted;
             unsigned int pstarti, pendi, pstartdriftedi, penddriftedi; //mcparticle indices for starts and ends in tpc or drifted volumes
             double plen = length(*pPart, mcstart, mcend, pstarti, pendi);
             double plendrifted = driftedLength(detProp, *pPart, mcstartdrifted, mcenddrifted, pstartdriftedi, penddriftedi);
 
             bool isActive = plen != 0;
             bool isDrifted = plendrifted!= 0;
             if (plen) ++active;
 
             fData->process_primary[geant_particle] = int(isPrimary);
             fData->processname[geant_particle]= pPart->Process();
             fData->Mother[geant_particle]=pPart->Mother();
             fData->TrackId[geant_particle]=TrackID;
             fData->pdg[geant_particle]=pPart->PdgCode();
             fData->status[geant_particle] = pPart->StatusCode();
             fData->Eng[geant_particle]=pPart->E();
             fData->EndE[geant_particle]=pPart->EndE();
             fData->Mass[geant_particle]=pPart->Mass();
             fData->Px[geant_particle]=pPart->Px();
             fData->Py[geant_particle]=pPart->Py();
             fData->Pz[geant_particle]=pPart->Pz();
             fData->P[geant_particle]=pPart->Momentum().Vect().Mag();
             fData->StartPointx[geant_particle]=pPart->Vx();
             fData->StartPointy[geant_particle]=pPart->Vy();
             fData->StartPointz[geant_particle]=pPart->Vz();
             fData->StartT[geant_particle] = pPart->T();
             fData->EndPointx[geant_particle]=pPart->EndPosition()[0];
             fData->EndPointy[geant_particle]=pPart->EndPosition()[1];
             fData->EndPointz[geant_particle]=pPart->EndPosition()[2];
             fData->EndT[geant_particle] = pPart->EndT();
             fData->theta[geant_particle] = pPart->Momentum().Theta();
             fData->phi[geant_particle] = pPart->Momentum().Phi();
             fData->theta_xz[geant_particle] = std::atan2(pPart->Px(), pPart->Pz());
             fData->theta_yz[geant_particle] = std::atan2(pPart->Py(), pPart->Pz());
             fData->pathlen[geant_particle]  = plen;
             fData->pathlen_drifted[geant_particle]  = plendrifted;
             fData->NumberDaughters[geant_particle]=pPart->NumberDaughters();
             fData->inTPCActive[geant_particle] = int(isActive);
             fData->inTPCDrifted[geant_particle] = int(isDrifted);
             art::Ptr<simb::MCTruth> const& mc_truth = pi_serv->ParticleToMCTruth_P(pPart);
             if (mc_truth){
               fData->origin[geant_particle] = mc_truth->Origin();
               fData->MCTruthIndex[geant_particle] = mc_truth.key();
             }
             if (isActive){
               fData->StartPointx_tpcAV[geant_particle] = mcstart.X();
               fData->StartPointy_tpcAV[geant_particle] = mcstart.Y();
               fData->StartPointz_tpcAV[geant_particle] = mcstart.Z();
               fData->StartT_tpcAV[geant_particle] = mcstart.T();
               fData->StartE_tpcAV[geant_particle] = pPart->E(pstarti);
               fData->StartP_tpcAV[geant_particle] = pPart->P(pstarti);
               fData->StartPx_tpcAV[geant_particle] = pPart->Px(pstarti);
               fData->StartPy_tpcAV[geant_particle] = pPart->Py(pstarti);
               fData->StartPz_tpcAV[geant_particle] = pPart->Pz(pstarti);
               fData->EndPointx_tpcAV[geant_particle] = mcend.X();
               fData->EndPointy_tpcAV[geant_particle] = mcend.Y();
               fData->EndPointz_tpcAV[geant_particle] = mcend.Z();
               fData->EndT_tpcAV[geant_particle] = mcend.T();
               fData->EndE_tpcAV[geant_particle] = pPart->E(pendi);
               fData->EndP_tpcAV[geant_particle] = pPart->P(pendi);
               fData->EndPx_tpcAV[geant_particle] = pPart->Px(pendi);
               fData->EndPy_tpcAV[geant_particle] = pPart->Py(pendi);
               fData->EndPz_tpcAV[geant_particle] = pPart->Pz(pendi);
             }
             if (isDrifted){
               fData->StartPointx_drifted[geant_particle] = mcstartdrifted.X();
               fData->StartPointy_drifted[geant_particle] = mcstartdrifted.Y();
               fData->StartPointz_drifted[geant_particle] = mcstartdrifted.Z();
               fData->StartT_drifted[geant_particle] = mcstartdrifted.T();
               fData->StartE_drifted[geant_particle] = pPart->E(pstartdriftedi);
               fData->StartP_drifted[geant_particle] = pPart->P(pstartdriftedi);
               fData->StartPx_drifted[geant_particle] = pPart->Px(pstartdriftedi);
               fData->StartPy_drifted[geant_particle] = pPart->Py(pstartdriftedi);
               fData->StartPz_drifted[geant_particle] = pPart->Pz(pstartdriftedi);
               fData->EndPointx_drifted[geant_particle] = mcenddrifted.X();
               fData->EndPointy_drifted[geant_particle] = mcenddrifted.Y();
               fData->EndPointz_drifted[geant_particle] = mcenddrifted.Z();
               fData->EndT_drifted[geant_particle] = mcenddrifted.T();
               fData->EndE_drifted[geant_particle] = pPart->E(penddriftedi);
               fData->EndP_drifted[geant_particle] = pPart->P(penddriftedi);
               fData->EndPx_drifted[geant_particle] = pPart->Px(penddriftedi);
               fData->EndPy_drifted[geant_particle] = pPart->Py(penddriftedi);
               fData->EndPz_drifted[geant_particle] = pPart->Pz(penddriftedi);
             }
 
             //access auxiliary detector parameters
             if (fSaveAuxDetInfo) {
               unsigned short nAD = 0; // number of cells that particle hit
 
               // find deposit of this particle in each of the detector cells
               for (const sim::AuxDetSimChannel* c: fAuxDetSimChannels) {
 
                 // find if this cell has a contribution (IDE) from this particle,
                 // and which one
                 const std::vector<sim::AuxDetIDE>& setOfIDEs = c->AuxDetIDEs();
                 // using a C++ "lambda" function here; this one:
                 // - sees only TrackID from the current scope
                 // - takes one parameter: the AuxDetIDE to be tested
                 // - returns if that IDE belongs to the track we are looking for
                 std::vector<sim::AuxDetIDE>::const_iterator iIDE
                   = std::find_if(
                                  setOfIDEs.begin(), setOfIDEs.end(),
                                  [TrackID](const sim::AuxDetIDE& IDE){ return IDE.trackID == TrackID; }
                                  );
                 if (iIDE == setOfIDEs.end()) continue;
 
                 // now iIDE points to the energy released by the track #i (TrackID)
 
                 // look for IDE with matching trackID
                 // find trackIDs stored in setOfIDEs with the same trackID, but negative,
                 // this is an untracked particle who's energy should be added as deposited by this original trackID
                 float totalE = 0.; // total energy deposited around by the GEANT particle in this cell
                 for(const auto& adtracks: setOfIDEs) {
                   if( fabs(adtracks.trackID) == TrackID )
                     totalE += adtracks.energyDeposited;
                 } // for
 
                 // fill the structure
                 if (nAD < kMaxAuxDets) {
                   fData->AuxDetID[geant_particle][nAD] = c->AuxDetID();
                   fData->entryX[geant_particle][nAD]   = iIDE->entryX;
                   fData->entryY[geant_particle][nAD]   = iIDE->entryY;
                   fData->entryZ[geant_particle][nAD]   = iIDE->entryZ;
                   fData->entryT[geant_particle][nAD]   = iIDE->entryT;
                   fData->exitX[geant_particle][nAD]    = iIDE->exitX;
                   fData->exitY[geant_particle][nAD]    = iIDE->exitY;
                   fData->exitZ[geant_particle][nAD]    = iIDE->exitZ;
                   fData->exitT[geant_particle][nAD]    = iIDE->exitT;
                   fData->exitPx[geant_particle][nAD]   = iIDE->exitMomentumX;
                   fData->exitPy[geant_particle][nAD]   = iIDE->exitMomentumY;
                   fData->exitPz[geant_particle][nAD]   = iIDE->exitMomentumZ;
                   fData->CombinedEnergyDep[geant_particle][nAD] = totalE;
                 }
                 ++nAD;
               } // for aux det sim channels
               fData->NAuxDets[geant_particle] = nAD;
 
               if (nAD > kMaxAuxDets) {
                 // got this error? consider increasing kMaxAuxDets
                 mf::LogError("AnalysisTree:limits")
                   << "particle #" << iPart
                   << " touches " << nAD << " auxiliary detector cells, only "
                   << kMaxAuxDets << " of them are saved in the tree";
               } // if too many detector cells
             } // if (fSaveAuxDetInfo)
 
             ++geant_particle;
           }
           else if (iPart == fData->GetMaxGEANTparticles()) {
             // got this error? it might be a bug,
             // since the structure should have enough room for everything
             mf::LogError("AnalysisTree:limits") << "event has "
                                                 << plist.size() << " MC particles, only "
                                                 << fData->GetMaxGEANTparticles() << " will be stored in tree";
           }
         } // for particles
 
         fData->geant_list_size_in_tpcAV = active;
         fData->no_primaries = primary;
         fData->geant_list_size = geant_particle;
         fData->processname.resize(geant_particle);
         MF_LOG_DEBUG("AnalysisTree")
           << "Counted "
           << fData->geant_list_size << " GEANT particles ("
           << fData->geant_list_size_in_tpcAV << " in AV), "
           << fData->no_primaries << " primaries, "
           << fData->genie_no_primaries << " GENIE particles";
 
         FillWith(fData->MergedId, 0);
 
         // for each particle, consider all the direct ancestors with the same
         // PDG ID, and mark them as belonging to the same "group"
         // (having the same MergedId)
         /* turn off for now
            int currentMergedId = 1;
            for(size_t iPart = 0; iPart < geant_particle; ++iPart){
            // if the particle already belongs to a group, don't bother
            if (fData->MergedId[iPart]) continue;
            // the particle starts its own group
            fData->MergedId[iPart] = currentMergedId;
            int currentMotherTrackId = fData->Mother[iPart];
            while (currentMotherTrackId > 0) {
            if (TrackIDtoIndex.find(currentMotherTrackId)==TrackIDtoIndex.end()) break;
            size_t gindex = TrackIDtoIndex[currentMotherTrackId];
            if (gindex<0||gindex>=plist.size()) break;
            // if the mother particle is of a different type,
            // don't bother with iPart ancestry any further
            if (gpdg[gindex]!=fData->pdg[iPart]) break;
            if (TrackIDtoIndex.find(currentMotherTrackId)!=TrackIDtoIndex.end()){
            size_t igeantMother = TrackIDtoIndex[currentMotherTrackId];
            if (igeantMother>=0&&igeantMother<geant_particle){
            fData->MergedId[igeantMother] = currentMergedId;
            }
            }
            currentMotherTrackId = gmother[gindex];
            }
            ++currentMergedId;
            }// for merging check
         */
       } // if (fSaveGeantInfo)
 
       // Now we have the GEANT info, see if we can match the protoDUNE generator particles
       if(fSaveProtoInfo){
         for(Int_t prt = 0; prt < nProtoPrimaries; ++prt){
           for(Int_t gnt = 0; gnt < fData->geant_list_size; ++gnt){
 //            if(fData->proto_pdg[prt] == fData->pdg[gnt] && fData->proto_px[prt] == fData->Px[gnt]){
              if(fData->proto_pdg[prt] == fData->pdg[gnt] && std::fabs(fData->proto_px[prt] - fData->Px[gnt]) < 0.0001){
               fData->proto_geantTrackID[prt] = fData->TrackId[gnt];
               fData->proto_geantIndex[prt] = gnt;
               break;
             }
           } // End GEANT loop
         } // End protoDUNE generator loop
       } // End ProtoDUNE generator if statement
 
     }//if (mcevts_truth)
   }//if (isMC){
   fData->taulife = detProp.ElectronLifetime();
   fTree->Fill();
 
   if (mf::isDebugEnabled()) {
     // use mf::LogDebug instead of MF_LOG_DEBUG because we reuse it in many lines;
     // thus, we protect this part of the code with the line above
     mf::LogDebug logStream("AnalysisTreeStructure");
     logStream
       << "Tree data structure contains:"
       << "\n - " << fData->no_hits << " hits (" << fData->GetMaxHits() << ")"
       << "\n - " << fData->genie_no_primaries << " genie primaries (" << fData->GetMaxGeniePrimaries() << ")"
       << "\n - " << fData->geant_list_size << " GEANT particles (" << fData->GetMaxGEANTparticles() << "), "
       << fData->no_primaries << " primaries"
       << "\n - " << fData->geant_list_size_in_tpcAV << " GEANT particles in AV "
       << "\n - " << ((int) fData->kNTracker) << " trackers:"
       ;
 
     size_t iTracker = 0;
     for (auto tracker = fData->TrackData.cbegin();
          tracker != fData->TrackData.cend(); ++tracker, ++iTracker
          ) {
       logStream
         << "\n -> " << tracker->ntracks << " " << fTrackModuleLabel[iTracker]
         << " tracks (" << tracker->GetMaxTracks() << ")"
         ;
       for (int iTrk = 0; iTrk < tracker->ntracks; ++iTrk) {
         logStream << "\n    [" << iTrk << "] "<< tracker->ntrkhits[iTrk][0];
         for (size_t ipl = 1; ipl < tracker->GetMaxPlanesPerTrack(iTrk); ++ipl)
           logStream << " + " << tracker->ntrkhits[iTrk][ipl];
         logStream << " hits (" << tracker->GetMaxHitsPerTrack(iTrk, 0);
         for (size_t ipl = 1; ipl < tracker->GetMaxPlanesPerTrack(iTrk); ++ipl)
           logStream << " + " << tracker->GetMaxHitsPerTrack(iTrk, ipl);
         logStream << ")";
       } // for tracks
     } // for trackers
   } // if logging enabled
 
 
   // if we don't use a permanent buffer (which can be huge),
   // delete the current buffer, and we'll create a new one on the next event
   if (!fUseBuffer) {
     MF_LOG_DEBUG("AnalysisTreeStructure") << "Freeing the tree data structure";
     DestroyData();
   }
 } // dune::AnalysisTree::analyze()

double AnalysisTree::bdist ( const TVector3 & pos )

private

Definition at line 5491 of file AnalysisTree_module.cc.

 {
   double d1 = -ActiveBounds[0] + pos.X();
   double d2 =  ActiveBounds[1] - pos.X();
   double d3 = -ActiveBounds[2] + pos.Y();
   double d4 =  ActiveBounds[3] - pos.Y();
   double d5 = -ActiveBounds[4] + pos.Z();
   double d6 =  ActiveBounds[5] - pos.Z();
 
   double Xmin = std::min(d1, d2);
   double result = 0;
   if (Xmin<0) result = std::min(std::min(std::min( d3, d4), d5), d6);         // - FIXME Passing uncorrected hits means X positions are very wrong ( outside of ActiveVolume )
   else result = std::min(std::min(std::min(std::min(Xmin, d3), d4), d5), d6);
   if (result<-1) result = -1;
   /*
   std::cout << "\n" << std::endl;
   std::cout << "-"<<ActiveBounds[0]<<" + "<<pos.X()<< " = " << d1 << "\n"
             <<      ActiveBounds[1]<<" - "<<pos.X()<< " = " << d2 << "\n"
             << "-"<<ActiveBounds[2]<<" + "<<pos.Y()<< " = " << d3 << "\n"
             <<      ActiveBounds[3]<<" - "<<pos.Y()<< " = " << d4 << "\n"
             << "-"<<ActiveBounds[4]<<" + "<<pos.Z()<< " = " << d5 << "\n"
             <<      ActiveBounds[5]<<" - "<<pos.Z()<< " = " << d6 << "\n"
             << "And the Minimum is " << result << std::endl;
   */
   return result;
 }

void AnalysisTree::beginSubRun ( const art::SubRun & sr )

Definition at line 3490 of file AnalysisTree_module.cc.

 {
 
 //  auto potListHandle = sr.getHandle< sumdata::POTSummary >(fPOTModuleLabel);
 //  if (potListHandle)
 //    SubRunData.pot=potListHandle->totpot;
 //  else
 //    SubRunData.pot=0.;
 
 }

void dune::AnalysisTree::CheckData ( std::string caller ) const

inlineprivate

Helper function: throws if no data structure is available.

Definition at line 1418 of file AnalysisTree_module.cc.

     {
       if (fData) return;
       throw art::Exception(art::errors::LogicError)
         << "AnalysisTree::" << caller << ": no data";
     } // CheckData()

void dune::AnalysisTree::CheckTree ( std::string caller ) const

inlineprivate

Helper function: throws if no tree is available.

Definition at line 1425 of file AnalysisTree_module.cc.

     {
       if (fTree) return;
       throw art::Exception(art::errors::LogicError)
         << "AnalysisTree::" << caller << ": no tree";
     } // CheckData()

void dune::AnalysisTree::CreateData ( bool bClearData = false )

inlineprivate

Creates the structure for the tree data; optionally initializes it.

Definition at line 1321 of file AnalysisTree_module.cc.

     {
       if (!fData) {
         fData.reset
           (new AnalysisTreeDataStruct(GetNTrackers(), GetNVertexAlgos(), GetShowerAlgos()));
         fData->SetBits(AnalysisTreeDataStruct::tdCry,    !fSaveCryInfo);
         fData->SetBits(AnalysisTreeDataStruct::tdGenie,  !fSaveGenieInfo);
         fData->SetBits(AnalysisTreeDataStruct::tdProto,  !fSaveProtoInfo);
         fData->SetBits(AnalysisTreeDataStruct::tdGeant,  !fSaveGeantInfo);
         fData->SetBits(AnalysisTreeDataStruct::tdMCshwr, !fSaveMCShowerInfo);
         fData->SetBits(AnalysisTreeDataStruct::tdMCtrk,  !fSaveMCTrackInfo);
         fData->SetBits(AnalysisTreeDataStruct::tdHit,    !fSaveHitInfo);
         fData->SetBits(AnalysisTreeDataStruct::tdRawDigit,    !fSaveRawDigitInfo);
         fData->SetBits(AnalysisTreeDataStruct::tdFlash,  !fSaveFlashInfo);
         fData->SetBits(AnalysisTreeDataStruct::tdCount,  !fSaveExternCounterInfo);
         fData->SetBits(AnalysisTreeDataStruct::tdShower, !fSaveShowerInfo);
         fData->SetBits(AnalysisTreeDataStruct::tdCluster,!fSaveClusterInfo);
         fData->SetBits(AnalysisTreeDataStruct::tdPandoraNuVertex,!fSavePandoraNuVertexInfo);
         fData->SetBits(AnalysisTreeDataStruct::tdTrack,  !fSaveTrackInfo);
         fData->SetBits(AnalysisTreeDataStruct::tdVertex, !fSaveVertexInfo);
         fData->SetBits(AnalysisTreeDataStruct::tdAuxDet, !fSaveAuxDetInfo);
         fData->SetBits(AnalysisTreeDataStruct::tdPFParticle, !fSavePFParticleInfo);
   fData->SetBits(AnalysisTreeDataStruct::tdSpacePoint, !fSaveSpacePointSolverInfo);
   fData->SetBits(AnalysisTreeDataStruct::tdCnn, !fSaveCnnInfo);
       }
       else {
         fData->SetTrackers(GetNTrackers());
         fData->SetVertexAlgos(GetNVertexAlgos());
         fData->SetShowerAlgos(GetShowerAlgos());
 
         if (bClearData) fData->Clear();
       }
     } // CreateData()

void AnalysisTree::CreateTree ( bool bClearData = false )

private

Create the output tree and the data structures, if needed.

Definition at line 3472 of file AnalysisTree_module.cc.

                                                    {
   if (!fTree) {
     art::ServiceHandle<art::TFileService> tfs;
     fTree = tfs->make<TTree>("anatree","analysis tree");
   }
   if (!fPOT) {
     art::ServiceHandle<art::TFileService> tfs;
     fPOT = tfs->make<TTree>("pottree","pot tree");
     fPOT->Branch("pot",&SubRunData.pot,"pot/D");
     fPOT->Branch("potbnbETOR860",&SubRunData.potbnbETOR860,"potbnbETOR860/D");
     fPOT->Branch("potbnbETOR875",&SubRunData.potbnbETOR875,"potbnbETOR875/D");
     fPOT->Branch("potnumiETORTGT",&SubRunData.potnumiETORTGT,"potnumiETORTGT/D");
   }
   CreateData(bClearData);
   SetAddresses();
 } // dune::AnalysisTree::CreateTree()

void dune::AnalysisTree::DestroyData ( )

inlineprivate

Destroy the local buffers (existing branches will point to invalid address!)

Definition at line 1415 of file AnalysisTree_module.cc.

1415 { fData.reset(); }

dune::AnalysisTree::fData

std::unique_ptr< AnalysisTreeDataStruct > fData

Definition: AnalysisTree_module.cc:1243

double AnalysisTree::driftedLength	(	detinfo::DetectorPropertiesData const &	detProp,
		const simb::MCParticle &	part,
		TLorentzVector &	start,
		TLorentzVector &	end,
		unsigned int &	starti,
		unsigned int &	endi
	)

private

Definition at line 5576 of file AnalysisTree_module.cc.

 {
   auto const* geom = lar::providerFrom<geo::Geometry>();
 
   //compute the drift x range
   double vDrift = detProp.DriftVelocity()*1e-3; //cm/ns
   double xrange[2] = {DBL_MAX, -DBL_MAX };
   for (unsigned int c=0; c<geom->Ncryostats(); ++c) {
     for (unsigned int t=0; t<geom->NTPC(c); ++t) {
       double Xat0 = detProp.ConvertTicksToX(0,0,t,c);
       double XatT = detProp.ConvertTicksToX(detProp.NumberTimeSamples(),0,t,c);
       xrange[0] = std::min(std::min(Xat0, XatT), xrange[0]);
       xrange[1] = std::max(std::max(Xat0, XatT), xrange[1]);
     }
   }
 
   double result = 0.;
   TVector3 disp;
   bool first = true;
 
   for(unsigned int i = 0; i < p.NumberTrajectoryPoints(); ++i) {
     // check if the particle is inside a TPC
     if (p.Vx(i) >= ActiveBounds[0] && p.Vx(i) <= ActiveBounds[1] &&
         p.Vy(i) >= ActiveBounds[2] && p.Vy(i) <= ActiveBounds[3] &&
         p.Vz(i) >= ActiveBounds[4] && p.Vz(i) <= ActiveBounds[5]){
       // Doing some manual shifting to account for
       // an interaction not occuring with the beam dump
       // we will reconstruct an x distance different from
       // where the particle actually passed to to the time
       // being different from in-spill interactions
       double newX = p.Vx(i)+(p.T(i)*vDrift);
       if (newX < xrange[0] || newX > xrange[1]) continue;
       TLorentzVector pos(newX,p.Vy(i),p.Vz(i),p.T());
       if(first){
         start = pos;
         starti=i;
         first = false;
       }
       else {
         disp -= pos.Vect();
         result += disp.Mag();
       }
       disp = pos.Vect();
       end = pos;
       endi = i;
     }
   }
   return result;
 }

double AnalysisTree::driftedLength	(	detinfo::DetectorPropertiesData const &	detProp,
		const sim::MCTrack &	mctrack,
		TLorentzVector &	tpcstart,
		TLorentzVector &	tpcend,
		TLorentzVector &	tpcmom
	)

private

Definition at line 5525 of file AnalysisTree_module.cc.

                                                                                                                                              {
   auto const* geom = lar::providerFrom<geo::Geometry>();
 
   //compute the drift x range
   double vDrift = detProp.DriftVelocity()*1e-3; //cm/ns
   double xrange[2] = {DBL_MAX, -DBL_MAX };
   for (unsigned int c=0; c<geom->Ncryostats(); ++c) {
     for (unsigned int t=0; t<geom->NTPC(c); ++t) {
       double Xat0 = detProp.ConvertTicksToX(0,0,t,c);
       double XatT = detProp.ConvertTicksToX(detProp.NumberTimeSamples(),0,t,c);
       xrange[0] = std::min(std::min(Xat0, XatT), xrange[0]);
       xrange[1] = std::max(std::max(Xat0, XatT), xrange[1]);
     }
   }
 
   double result = 0.;
   TVector3 disp;
   bool first = true;
 
   for(auto step: mctrack) {
     // check if the particle is inside a TPC
     if (step.X() >= ActiveBounds[0] && step.X() <= ActiveBounds[1] &&
         step.Y() >= ActiveBounds[2] && step.Y() <= ActiveBounds[3] &&
         step.Z() >= ActiveBounds[4] && step.Z() <= ActiveBounds[5] ){
       // Doing some manual shifting to account for
       // an interaction not occuring with the beam dump
       // we will reconstruct an x distance different from
       // where the particle actually passed to to the time
       // being different from in-spill interactions
       double newX = step.X()+(step.T()*vDrift);
       if (newX < xrange[0] || newX > xrange[1]) continue;
 
       TLorentzVector pos(newX,step.Y(),step.Z(),step.T());
       if(first){
         tpcstart = pos;
         tpcmom = step.Momentum();
         first = false;
       }
       else {
         disp -= pos.Vect();
         result += disp.Mag();
       }
       disp = pos.Vect();
       tpcend = pos;
     }
   }
   return result;
 }

void AnalysisTree::endSubRun ( const art::SubRun & sr )

Definition at line 3501 of file AnalysisTree_module.cc.

 {
 
   auto potListHandle = sr.getHandle< sumdata::POTSummary >(fPOTModuleLabel);
   if (potListHandle)
     SubRunData.pot=potListHandle->totpot;
   else
     SubRunData.pot=0.;
 
   art::InputTag itag1("beamdata","bnbETOR860");
   auto potSummaryHandlebnbETOR860 = sr.getHandle<sumdata::POTSummary>(itag1);
   if (potSummaryHandlebnbETOR860){
     SubRunData.potbnbETOR860 = potSummaryHandlebnbETOR860->totpot;
   }
   else
     SubRunData.potbnbETOR860 = 0;
 
   art::InputTag itag2("beamdata","bnbETOR875");
   auto potSummaryHandlebnbETOR875 = sr.getHandle<sumdata::POTSummary>(itag2);
   if (potSummaryHandlebnbETOR875){
     SubRunData.potbnbETOR875 = potSummaryHandlebnbETOR875->totpot;
   }
   else
     SubRunData.potbnbETOR875 = 0;
 
   art::InputTag itag3("beamdata","numiETORTGT");
   auto potSummaryHandlenumiETORTGT = sr.getHandle<sumdata::POTSummary>(itag3);
   if (potSummaryHandlenumiETORTGT){
     SubRunData.potnumiETORTGT = potSummaryHandlenumiETORTGT->totpot;
   }
   else
     SubRunData.potnumiETORTGT = 0;
 
   if (fPOT) fPOT->Fill();
 
 }

void AnalysisTree::FillShower	(	AnalysisTreeDataStruct::ShowerDataStruct &	showerData,
		size_t	iShower,
		recob::Shower const &	showers,
		const bool	fSavePFParticleInfo,
		const std::map< Short_t, Short_t > &	showerIDtoPFParticleIDMap
	)		const

private

Stores the information of shower in slot iShower of showerData.

Definition at line 5366 of file AnalysisTree_module.cc.

                                              {
 
   showerData.showerID[iShower]        = shower.ID();
   showerData.shwr_bestplane[iShower]  = shower.best_plane();
   showerData.shwr_length[iShower]     = shower.Length();
 
   TVector3 const& dir_start = shower.Direction();
   showerData.shwr_startdcosx[iShower] = dir_start.X();
   showerData.shwr_startdcosy[iShower] = dir_start.Y();
   showerData.shwr_startdcosz[iShower] = dir_start.Z();
 
   TVector3 const& pos_start = shower.ShowerStart();
   showerData.shwr_startx[iShower]     = pos_start.X();
   showerData.shwr_starty[iShower]     = pos_start.Y();
   showerData.shwr_startz[iShower]     = pos_start.Z();
 
   if (fSavePFParticleInfo) {
     auto mapIter = showerIDtoPFParticleIDMap.find(shower.ID());
     if (mapIter != showerIDtoPFParticleIDMap.end()) {
       // This vertex has a corresponding PFParticle.
       showerData.shwr_hasPFParticle[iShower] = 1;
       showerData.shwr_PFParticleID[iShower] = mapIter->second;
     }
     else
       showerData.shwr_hasPFParticle[iShower] = 0;
   }
 
   if (shower.Energy().size() == kNplanes)
     std::copy_n
       (shower.Energy().begin(),    kNplanes, &showerData.shwr_totEng[iShower][0]);
   if (shower.dEdx().size() == kNplanes)
     std::copy_n
       (shower.dEdx().begin(),      kNplanes, &showerData.shwr_dedx[iShower][0]);
   if (shower.MIPEnergy().size() == kNplanes)
     std::copy_n
       (shower.MIPEnergy().begin(), kNplanes, &showerData.shwr_mipEng[iShower][0]);
 
 } // dune::AnalysisTree::FillShower()

void AnalysisTree::FillShowers	(	AnalysisTreeDataStruct::ShowerDataStruct &	showerData,
		std::vector< recob::Shower > const &	showers,
		const bool	fSavePFParticleInfo,
		const std::map< Short_t, Short_t > &	showerIDtoPFParticleIDMap
	)		const

private

Stores the information of all showers into showerData.

Definition at line 5409 of file AnalysisTree_module.cc.

                                               {
 
   const size_t NShowers = showers.size();
 
   //
   // prepare the data structures, the tree and the connection between them
   //
 
   // allocate enough space for this number of showers
   // (but at least for one of them!)
   showerData.SetMaxShowers(std::max(NShowers, (size_t) 1));
   showerData.Clear(); // clear all the data
 
   // now set the tree addresses to the newly allocated memory;
   // this creates the tree branches in case they are not there yet
   showerData.SetAddresses(fTree);
   if (NShowers > showerData.GetMaxShowers()) {
     // got this error? it might be a bug,
     // since we are supposed to have allocated enough space to fit all showers
     mf::LogError("AnalysisTree:limits") << "event has " << NShowers
                                         << " " << showerData.Name() << " showers, only "
                                         << showerData.GetMaxShowers() << " stored in tree";
   }
 
   //
   // now set the data
   //
   // set the record of the number of showers
   // (data structures are already properly resized)
   showerData.nshowers = (Short_t) NShowers;
 
   // set all the showers one by one
   for (size_t i = 0; i < NShowers; ++i) FillShower(showerData, i, showers[i], fSavePFParticleInfo, showerIDtoPFParticleIDMap);
 
 } // dune::AnalysisTree::FillShowers()

size_t dune::AnalysisTree::GetNShowerAlgos ( ) const

inlineprivate

Returns the number of shower algorithms configured.

Definition at line 1314 of file AnalysisTree_module.cc.

1314 { return fShowerModuleLabel.size(); }

dune::AnalysisTree::fShowerModuleLabel

std::vector< std::string > fShowerModuleLabel

Definition: AnalysisTree_module.cc:1266

size_t dune::AnalysisTree::GetNTrackers ( ) const

inlineprivate

Returns the number of trackers configured.

Definition at line 1309 of file AnalysisTree_module.cc.

1309 { return fTrackModuleLabel.size(); }

dune::AnalysisTree::fTrackModuleLabel

std::vector< std::string > fTrackModuleLabel

Definition: AnalysisTree_module.cc:1263

size_t dune::AnalysisTree::GetNVertexAlgos ( ) const

inlineprivate

Definition at line 1311 of file AnalysisTree_module.cc.

1311 { return fVertexModuleLabel.size(); }

dune::AnalysisTree::fVertexModuleLabel

std::vector< std::string > fVertexModuleLabel

Definition: AnalysisTree_module.cc:1265

std::vector<std::string> dune::AnalysisTree::GetShowerAlgos ( ) const

inlineprivate

Returns the name of configured shower algorithms (converted to string)

Definition at line 1317 of file AnalysisTree_module.cc.

1318 { return { fShowerModuleLabel.begin(), fShowerModuleLabel.end() }; }

dune::AnalysisTree::fShowerModuleLabel

std::vector< std::string > fShowerModuleLabel

Definition: AnalysisTree_module.cc:1266

void AnalysisTree::HitsPurity	(	detinfo::DetectorClocksData const &	clockData,
		std::vector< art::Ptr< recob::Hit > > const &	hits,
		Int_t &	trackid,
		Float_t &	purity,
		double &	maxe
	)

private

Definition at line 5450 of file AnalysisTree_module.cc.

                                                                                                                              {
 
   trackid = -1;
   purity = -1;
 
   art::ServiceHandle<cheat::BackTrackerService> bt_serv;
 
   std::map<int,double> trkide;
 
   for(size_t h = 0; h < hits.size(); ++h){
 
     art::Ptr<recob::Hit> hit = hits[h];
     std::vector<sim::IDE> ides;
     //bt_serv->HitToSimIDEs(hit,ides);
     std::vector<sim::TrackIDE> eveIDs = bt_serv->HitToEveTrackIDEs(clockData, hit);
 
     for(size_t e = 0; e < eveIDs.size(); ++e){
       //std::cout<<h<<" "<<e<<" "<<eveIDs[e].trackID<<" "<<eveIDs[e].energy<<" "<<eveIDs[e].energyFrac<<std::endl;
       trkide[eveIDs[e].trackID] += eveIDs[e].energy;
     }
   }
 
   maxe = -1;
   double tote = 0;
   for (std::map<int,double>::iterator ii = trkide.begin(); ii!=trkide.end(); ++ii){
     tote += ii->second;
     if ((ii->second)>maxe){
       maxe = ii->second;
       trackid = ii->first;
     }
   }
 
   //std::cout << "the total energy of this reco track is: " << tote << std::endl;
 
   if (tote>0){
     purity = maxe/tote;
   }
 }

double AnalysisTree::length ( const recob::Track & track )

private

Definition at line 5519 of file AnalysisTree_module.cc.

 {
   return track.Length();
 }

double AnalysisTree::length	(	const simb::MCParticle &	part,
		TLorentzVector &	start,
		TLorentzVector &	end,
		unsigned int &	starti,
		unsigned int &	endi
	)

private

Definition at line 5628 of file AnalysisTree_module.cc.

 {
   double result = 0.;
   TVector3 disp;
   bool first = true;
 
   for(unsigned int i = 0; i < p.NumberTrajectoryPoints(); ++i) {
     // check if the particle is inside a TPC
     if (p.Vx(i) >= ActiveBounds[0] && p.Vx(i) <= ActiveBounds[1] && p.Vy(i) >= ActiveBounds[2] && p.Vy(i) <= ActiveBounds[3] && p.Vz(i) >= ActiveBounds[4] && p.Vz(i) <= ActiveBounds[5]){
       if(first){
         start = p.Position(i);
         first = false;
         starti = i;
       }else{
         disp -= p.Position(i).Vect();
         result += disp.Mag();
       }
       disp = p.Position(i).Vect();
       end = p.Position(i);
       endi = i;
     }
   }
   return result;
 }

void dune::AnalysisTree::SetAddresses ( )

inlineprivate

Sets the addresses of all the tree branches, creating the missing ones.

Definition at line 1356 of file AnalysisTree_module.cc.

     {
       CheckData(__func__); CheckTree(__func__);
       fData->SetAddresses
         (fTree, fTrackModuleLabel, fVertexModuleLabel, fShowerModuleLabel, isCosmics);
     } // SetAddresses()

void dune::AnalysisTree::SetPFParticleAddress ( )

inlineprivate

Sets the addresses of the tree branch of the PFParticle, creating it if missing

Definition at line 1405 of file AnalysisTree_module.cc.

     {
       CheckData(__func__); CheckTree(__func__);
       fData->GetPFParticleData().SetAddresses(fTree);
     } // SetPFParticleAddress()

void dune::AnalysisTree::SetShowerAddresses ( size_t iShower )

inlineprivate

Sets the addresses of all the tree branches of the specified shower algo, creating the missing ones

Definition at line 1392 of file AnalysisTree_module.cc.

     {
       CheckData(__func__); CheckTree(__func__);
       if (iShower >= fData->GetNShowerAlgos()) {
         throw art::Exception(art::errors::LogicError)
           << "AnalysisTree::SetShowerAddresses(): no shower algo #" << iShower
           << " (" << fData->GetNShowerAlgos() << " available)";
       }
       fData->GetShowerData(iShower).SetAddresses(fTree);
     } // SetShowerAddresses()

void dune::AnalysisTree::SetTrackerAddresses ( size_t iTracker )

inlineprivate

Sets the addresses of all the tree branches of the specified tracking algo, creating the missing ones

Definition at line 1365 of file AnalysisTree_module.cc.

     {
       CheckData(__func__); CheckTree(__func__);
       if (iTracker >= fData->GetNTrackers()) {
         throw art::Exception(art::errors::LogicError)
           << "AnalysisTree::SetTrackerAddresses(): no tracker #" << iTracker
           << " (" << fData->GetNTrackers() << " available)";
       }
       fData->GetTrackerData(iTracker)
         .SetAddresses(fTree, fTrackModuleLabel[iTracker], isCosmics);
     } // SetTrackerAddresses()

void dune::AnalysisTree::SetVertexAddresses ( size_t iVertexAlg )

inlineprivate

Definition at line 1378 of file AnalysisTree_module.cc.

     {
       CheckData(__func__); CheckTree(__func__);
       if (iVertexAlg >= fData->GetNVertexAlgos()) {
         throw art::Exception(art::errors::LogicError)
           << "AnalysisTree::SetVertexAddresses(): no vertex alg #" << iVertexAlg
           << " (" << fData->GetNVertexAlgos() << " available)";
       }
       fData->GetVertexData(iVertexAlg)
         .SetAddresses(fTree, fVertexModuleLabel[iVertexAlg], isCosmics);
     } // SetVertexAddresses()

Member Data Documentation

double dune::AnalysisTree::ActiveBounds[6]

private

Definition at line 1306 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::bIgnoreMissingShowers

private

whether to ignore missing shower information

Definition at line 1300 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fCalDataModuleLabel

private

Definition at line 1250 of file AnalysisTree_module.cc.

std::vector<std::string> dune::AnalysisTree::fCalorimetryModuleLabel

private

Definition at line 1267 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fClusterModuleLabel

private

Definition at line 1255 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fCnnModuleLabel

private

Definition at line 1262 of file AnalysisTree_module.cc.

std::vector<std::string> dune::AnalysisTree::fContainmentTaggerAssocLabel

private

Definition at line 1297 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fCosmicClusterTaggerAssocLabel

private

Definition at line 1274 of file AnalysisTree_module.cc.

std::vector<std::string> dune::AnalysisTree::fCosmicTaggerAssocLabel

private

whether to extract and save CNN information

Definition at line 1296 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fCryGenModuleLabel

private

Definition at line 1252 of file AnalysisTree_module.cc.

std::unique_ptr<AnalysisTreeDataStruct> dune::AnalysisTree::fData

private

Definition at line 1243 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fDigitModuleLabel

private

Definition at line 1247 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fExternalCounterModuleLabel

private

Definition at line 1258 of file AnalysisTree_module.cc.

std::vector<std::string> dune::AnalysisTree::fFlashMatchAssocLabel

private

Definition at line 1298 of file AnalysisTree_module.cc.

std::vector<std::string> dune::AnalysisTree::fFlashT0FinderLabel

private

Definition at line 1271 of file AnalysisTree_module.cc.

float dune::AnalysisTree::fG4minE

private

Energy threshold to save g4 particle info.

Definition at line 1304 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fG4ModuleLabel

private

Definition at line 1254 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fGenieGenModuleLabel

private

Definition at line 1251 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fHitsModuleLabel

private

Definition at line 1248 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fLArG4ModuleLabel

private

Definition at line 1249 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fMCShowerModuleLabel

private

Definition at line 1259 of file AnalysisTree_module.cc.

std::vector<std::string> dune::AnalysisTree::fMCT0FinderLabel

private

Definition at line 1272 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fMCTrackModuleLabel

private

Definition at line 1260 of file AnalysisTree_module.cc.

std::vector<std::string> dune::AnalysisTree::fMVAPIDShowerModuleLabel

private

Definition at line 1269 of file AnalysisTree_module.cc.

std::vector<std::string> dune::AnalysisTree::fMVAPIDTrackModuleLabel

private

Definition at line 1270 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fOpFlashModuleLabel

private

Definition at line 1257 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fPandoraNuVertexModuleLabel

private

Definition at line 1256 of file AnalysisTree_module.cc.

std::vector<std::string> dune::AnalysisTree::fParticleIDModuleLabel

private

Definition at line 1268 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fPFParticleModuleLabel

private

Definition at line 1264 of file AnalysisTree_module.cc.

TTree* dune::AnalysisTree::fPOT

private

Definition at line 1239 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fPOTModuleLabel

private

Definition at line 1273 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fProtoGenModuleLabel

private

Definition at line 1253 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveAuxDetInfo

private

whether to extract and save auxiliary detector data

Definition at line 1276 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveCaloCosmics

private

save calorimetry information for cosmics

Definition at line 1303 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveClusterInfo

private

whether to extract and save Vertex information

Definition at line 1287 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveCnnInfo

private

whether to extract and save SpacePointSolver information

Definition at line 1294 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveCryInfo

private

Definition at line 1277 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveExternCounterInfo

private

whether to extract and save Flash information

Definition at line 1290 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveFlashInfo

private

whether to extract and save nu vertex information from Pandora

Definition at line 1289 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveGeantInfo

private

whether to extract and save ProtDUNE beam simulation information

Definition at line 1280 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveGenieInfo

private

whether to extract and save CRY particle data

Definition at line 1278 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveHitInfo

private

whether to extract and save MC Track information

Definition at line 1283 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveMCShowerInfo

private

whether to extract and save Geant information

Definition at line 1281 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveMCTrackInfo

private

whether to extract and save MC Shower information

Definition at line 1282 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSavePandoraNuVertexInfo

private

whether to extract and save Cluster information

Definition at line 1288 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSavePFParticleInfo

private

whether to extract and save Shower information

Definition at line 1292 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveProtoInfo

private

whether to extract and save Genie information

Definition at line 1279 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveRawDigitInfo

private

whether to extract and save Hit information

Definition at line 1284 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveShowerInfo

private

whether to extract and save External Counter information

Definition at line 1291 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveSpacePointSolverInfo

private

whether to extract and save PFParticle information

Definition at line 1293 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveTrackInfo

private

whether to extract and save Raw Digit information

Definition at line 1285 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveVertexInfo

private

whether to extract and save Track information

Definition at line 1286 of file AnalysisTree_module.cc.

std::vector<std::string> dune::AnalysisTree::fShowerModuleLabel

private

Definition at line 1266 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fSpacePointSolverModuleLabel

private

Definition at line 1261 of file AnalysisTree_module.cc.

std::vector<std::string> dune::AnalysisTree::fTrackModuleLabel

private

Definition at line 1263 of file AnalysisTree_module.cc.

TTree* dune::AnalysisTree::fTree

private

Definition at line 1238 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fUseBuffer

private

whether to use a permanent buffer (faster, huge memory)

Definition at line 1275 of file AnalysisTree_module.cc.

std::vector<std::string> dune::AnalysisTree::fVertexModuleLabel

private

Definition at line 1265 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::isCosmics

private

if it contains cosmics

Definition at line 1302 of file AnalysisTree_module.cc.

AnalysisTreeDataStruct::SubRunData_t dune::AnalysisTree::SubRunData

private

Definition at line 1245 of file AnalysisTree_module.cc.

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

duneana/duneana/AnaTree/AnalysisTree_module.cc

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Configuration parameters

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation