dune::AnaRootParser Class Reference

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

Inheritance diagram for dune::AnaRootParser:

Public Member Functions
	AnaRootParser (fhicl::ParameterSet const &pset)
virtual	~AnaRootParser ()
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	HitsBackTrack (detinfo::DetectorClocksData const &clockData, art::Ptr< recob::Hit > const &hit, int &trackid, float &maxe, float &purity)
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)
int	CountHits (const art::Event &evt, const art::InputTag &which, unsigned int cryostat, unsigned int tpc, unsigned int plane)
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 (AnaRootParserDataStruct::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 (AnaRootParserDataStruct::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
std::unique_ptr< AnaRootParserDataStruct >	fData
AnaRootParserDataStruct::SubRunData_t	SubRunData
int	fLogLevel
short	fEventsPerSubrun
std::string	fRawDigitModuleLabel
std::string	fHitsModuleLabel
std::string	fPhotonPropS1ModuleLabel
std::string	fElecDriftModuleLabel
std::string	fCalDataModuleLabel
std::string	fGenieGenModuleLabel
std::string	fCryGenModuleLabel
std::string	fProtoGenModuleLabel
std::string	fG4ModuleLabel
std::string	fSimEnergyDepositTPCActiveLabel
std::string	fClusterModuleLabel
std::string	fPandoraNuVertexModuleLabel
std::string	fOpFlashModuleLabel
std::string	fExternalCounterModuleLabel
std::string	fMCShowerModuleLabel
std::string	fMCTrackModuleLabel
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	fIsMC
	whether to use a permanent buffer (faster, huge memory) More...
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	fSavePhotonInfo
	whether to extract and save ProtDUNE beam simulation information More...
bool	fSaveGeneratorInfo
	whether to extract and save collected photons More...
bool	fSaveGeantInfo
	whether to extract and save Geant information More...
bool	fSaveGeantInAVInfo
	whether to extract and save Geant information More...
bool	fSaveGeantTrajectoryInfo
	whether to extract and save Geant information More...
bool	fSaveSimEnergyDepositTPCActiveInfo
	whether to extract and save Geant information More...
bool	fSaveMCShowerInfo
	whether to extract and save Cluster 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	fSaveRecobWireInfo
	whether to extract and save raw digit information More...
bool	fSaveTrackInfo
	whether to extract and save recob wire 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...
std::vector< std::string >	fCosmicTaggerAssocLabel
	whether to extract and save PFParticle 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...
trkf::TrajectoryMCSFitter	fMCSFitter
double	ActiveBounds [6]
	objet holding the configuration of the uBoone MCS fitting alg More...

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 1374 of file AnaRootParser_module.cc.

Constructor & Destructor Documentation

dune::AnaRootParser::AnaRootParser ( fhicl::ParameterSet const &  pset )

explicit

Definition at line 3998 of file AnaRootParser_module.cc.

                                                              :
  EDAnalyzer(pset),
  fTree(nullptr),
  //  fPOT(nullptr),

  fLogLevel                 (pset.get< short >("LogLevel")        ),
  fEventsPerSubrun          (pset.get< short >("EventsPerSubrun")        ),
  fRawDigitModuleLabel         (pset.get< std::string >("RawDigitModuleLabel")        ),
  fHitsModuleLabel          (pset.get< std::string >("HitsModuleLabel")         ),
  fPhotonPropS1ModuleLabel  (pset.get< std::string >("PhotonPropS1ModuleLabel")     ),
  fElecDriftModuleLabel     (pset.get< std::string >("ElecDriftModuleLabel")     ),
  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")           ),
  fSimEnergyDepositTPCActiveLabel    (pset.get< std::string >("SimEnergyDepositTPCActiveLabel")           ),
  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")     ),
  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")),
  fIsMC                     (pset.get< bool >("IsMC", false)),
  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)),
  fSavePhotonInfo           (pset.get< bool >("SavePhotonInfo", false)),
  fSaveGeneratorInfo        (pset.get< bool >("SaveGeneratorInfo", false)),
  fSaveGeantInfo            (pset.get< bool >("SaveGeantInfo", false)),
  fSaveGeantInAVInfo        (pset.get< bool >("SaveGeantInAVInfo", false)),
  fSaveGeantTrajectoryInfo  (pset.get< bool >("SaveGeantTrajectoryInfo", false)),
  fSaveSimEnergyDepositTPCActiveInfo (pset.get< bool >("SaveSimEnergyDepositTPCActiveInfo", 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)),
  fSaveRecobWireInfo        (pset.get< bool >("SaveRecobWireInfo", 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)),
  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)),
  fMCSFitter( pset.get<fhicl::ParameterSet>("TrajMCSFitter") )

{

  if (fSavePFParticleInfo) fPFParticleModuleLabel = pset.get<std::string>("PFParticleModuleLabel");

  if (fSaveAuxDetInfo == true) fSaveGeantInfo = true;
//  if (fSaveRawDigitInfo == true) fSaveHitInfo = true;
  mf::LogInfo("AnaRootParser") << "Configuration:"
    << "\n  UseBuffers: " << std::boolalpha << fUseBuffer
    ;
  if (GetNTrackers() > kMaxTrackers) {
    throw art::Exception(art::errors::Configuration)
      << "AnaRootParser 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)
      << "AnaRootParser 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::AnaRootParser::AnaRootParser()

dune::AnaRootParser::~AnaRootParser ( )

virtual

Definition at line 4150 of file AnaRootParser_module.cc.

{
  DestroyData();
}

Member Function Documentation

void dune::AnaRootParser::analyze

(

const art::Event &

evt

)

read access to event

transfer the run and subrun data to the tree data object

Definition at line 4221 of file AnaRootParser_module.cc.

{

  //services
//  art::ServiceHandle<cheat::BackTrackerService> bt_serv;
  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;

  // collect the sizes which might be needed to resize the tree data structure:

  // RawDigit
  std::vector<art::Ptr<raw::RawDigit> > rawdigitlist;
  auto rawdigitVecHandle = evt.getHandle< std::vector<raw::RawDigit> >(fRawDigitModuleLabel);
  if (rawdigitVecHandle)
    art::fill_ptr_vector(rawdigitlist, rawdigitVecHandle);

  // RecobWire
  std::vector<art::Ptr<recob::Wire> > recobwirelist;
  auto recobwireVecHandle = evt.getHandle< std::vector<recob::Wire> >(fCalDataModuleLabel);
  if (recobwireVecHandle)
    art::fill_ptr_vector(recobwirelist, recobwireVecHandle);


  // recob::Wire
//  auto wireVecHandle = evt.getHandle< std::vector<recob::Wire> >(fCalDataModuleLabel);

  // * photons
  //std::vector<art::Ptr<sim::SimPhotonsLite> > photonlist;
  art::Handle< std::vector<sim::SimPhotonsLite> > photonHandle;
  if(fSavePhotonInfo) photonHandle = evt.getHandle< std::vector<sim::SimPhotonsLite> >(fPhotonPropS1ModuleLabel);
//    art::fill_ptr_vector(photonlist, photonHandle);

  // * 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);
  }

  // * 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 (fIsMC){
    mctruthListHandle = evt.getHandle< std::vector<simb::MCTruth> >(fGenieGenModuleLabel);
    if (mctruthListHandle)
      art::fill_ptr_vector(mclist, mctruthListHandle);
  }

  // SimEnergyDepositTPCActive
  int nSimEnergyDepositsTPCActive = 0;
  int nParticlesWithSimEnergyDepositsTPCActive = 0;
  int nLowestParticleIDSimEnergyDepositsTPCActive = 0;
  int nHighestParticleIDSimEnergyDepositsTPCActive = 0;

  std::vector<art::Ptr<sim::SimEnergyDeposit> > energyDepositTPCActivelist;
  if (fIsMC && fSaveSimEnergyDepositTPCActiveInfo)
  {
    auto energyDepositTPCActiveHandle = evt.getHandle< std::vector<sim::SimEnergyDeposit> >(fSimEnergyDepositTPCActiveLabel);
    if (energyDepositTPCActiveHandle)
    {
      art::fill_ptr_vector(energyDepositTPCActivelist,energyDepositTPCActiveHandle);
      nSimEnergyDepositsTPCActive = energyDepositTPCActivelist.size();

      std::vector<int> tempparticleID;
      tempparticleID.resize(nSimEnergyDepositsTPCActive);
      FillWith(tempparticleID, 0);

      for(int sedavit = 0; sedavit < nSimEnergyDepositsTPCActive; sedavit++)
      {
        if( std::find(tempparticleID.begin(), tempparticleID.end(), energyDepositTPCActivelist[sedavit]->TrackID()) == tempparticleID.end() ) //check if current particle ID is already in the vector. if true: not in vector
        {
          tempparticleID[nParticlesWithSimEnergyDepositsTPCActive] = energyDepositTPCActivelist[sedavit]->TrackID();
          nParticlesWithSimEnergyDepositsTPCActive++;
        }
        if(energyDepositTPCActivelist[sedavit]->TrackID() < nLowestParticleIDSimEnergyDepositsTPCActive) nLowestParticleIDSimEnergyDepositsTPCActive=energyDepositTPCActivelist[sedavit]->TrackID();
        if(energyDepositTPCActivelist[sedavit]->TrackID() > nHighestParticleIDSimEnergyDepositsTPCActive) nHighestParticleIDSimEnergyDepositsTPCActive=energyDepositTPCActivelist[sedavit]->TrackID();
      }

    }
  }

  // *MC truth cosmic generator information
  std::vector<art::Ptr<simb::MCTruth> > mclistcry;
  if (fIsMC && fSaveCryInfo){
    auto 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("AnaRootParser") << "Requested CRY information but none exists, hence not saving CRY information.";
    }
  }

  // ProtoDUNE beam generator information
  std::vector<art::Ptr<simb::MCTruth> > mclistproto;
  if(fIsMC && fSaveProtoInfo){
    auto 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("AnaRootParser") << "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 (fIsMC)
    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 (fIsMC)
    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 nPhotons=0;
  int nGeniePrimaries = 0, nGeneratorParticles = 0, nGEANTparticles = 0, nGEANTparticlesInAV = 0, nGEANTtrajectorysteps=0;

  art::Ptr<simb::MCTruth> mctruth;


  if (fIsMC) { //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();

      sim::ParticleList::const_iterator itPart = plist.begin(),
      pend = plist.end(); // iterator to pairs (track id, particle)

      std::string pri("primary");
      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";
        }
        nGEANTtrajectorysteps += pPart->NumberTrajectoryPoints();

        if(pPart->Mother() == 0 && pPart->Process() == pri) nGeneratorParticles++;

        TLorentzVector mcstart, mcend;
        unsigned int pstarti, pendi;
        double plen = length(*pPart, mcstart, mcend, pstarti, pendi);
        if( plen != 0) nGEANTparticlesInAV += 1;

      } // for particles

      // counting photons
      if(fSavePhotonInfo)
      {
        for ( auto const& pmt : (*photonHandle) )
        {
          std::map<int, int> PhotonsMap = pmt.DetectedPhotons;

          for(auto itphoton = PhotonsMap.begin(); itphoton!= PhotonsMap.end(); itphoton++)
          {
            for(int i = 0; i < itphoton->second ; i++)
            {
              nPhotons++;
            }
          }
        }
      }
    } // if have MC truth
    MF_LOG_DEBUG("AnaRootParser") << "Expected "
    << nGEANTparticles << " GEANT particles, "
    << nGeniePrimaries << " GENIE particles";
  } // if MC

CreateData(); // tracker data is created with default constructor
if (fSaveGenieInfo){  fData->ResizeGenie(nGeniePrimaries);}
if (fSaveCryInfo){  fData->ResizeCry(nCryPrimaries);}
if (fSaveProtoInfo){  fData->ResizeProto(nProtoPrimaries);}
if (fSaveGeneratorInfo){  fData->ResizeGenerator(nGeneratorParticles);}
if (fSaveGeantInfo){  fData->ResizeGEANT(nGEANTparticles);}
if (fSaveGeantInfo && fSaveGeantInAVInfo){  fData->ResizeGEANTInAV(nGEANTparticlesInAV);}
if (fSaveGeantTrajectoryInfo){  fData->ResizeGEANTTrajectory(nGEANTtrajectorysteps);}
if (fSaveSimEnergyDepositTPCActiveInfo) {fData->ResizeSimEnergyDepositsTPCActive(nSimEnergyDepositsTPCActive, nParticlesWithSimEnergyDepositsTPCActive);}
if (fSaveMCShowerInfo){  fData->ResizeMCShower(nMCShowers);}
if (fSaveMCTrackInfo){  fData->ResizeMCTrack(nMCTracks);}
if (fSavePhotonInfo){  fData->ResizePhotons(nPhotons);}

if(fLogLevel >= 1)
{
  std::cout << "nGeneratorParticles: " << nGeneratorParticles << std::endl;
  std::cout << "nGEANTparticles: " << nGEANTparticles << std::endl;
  std::cout << "nGEANTtrajectorysteps: " << nGEANTtrajectorysteps << std::endl;
  std::cout << "nGEANTparticlesInAV: " << nGEANTparticlesInAV << std::endl;
  std::cout << "nSimEnergyDepositsTPCActive: " << nSimEnergyDepositsTPCActive << std::endl;
  std::cout << "nParticlesWithSimEnergyDepositsTPCActive: " << nParticlesWithSimEnergyDepositsTPCActive << std::endl;
  std::cout << "nLowestParticleIDSimEnergyDepositsTPCActive: " << nLowestParticleIDSimEnergyDepositsTPCActive << std::endl;
  std::cout << "nHighestParticleIDSimEnergyDepositsTPCActive: " << nHighestParticleIDSimEnergyDepositsTPCActive << std::endl;
  std::cout << "nPhotons: " << nPhotons << std::endl;
}

  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 NChannels = rawdigitlist.size(); //number of channels holding raw waveforms
  const size_t NRecoChannels = recobwirelist.size(); //number of channels holding ROI's
  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(!fIsMC);

  //  raw trigger
  std::vector<art::Ptr<raw::Trigger> > triggerlist;
  auto triggerListHandle = evt.getHandle< std::vector<raw::Trigger>>(fRawDigitModuleLabel);
  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 AnaRootParser'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("AnaRootParser")
          << "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


// 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);

std::vector<const sim::AuxDetSimChannel*> fAuxDetSimChannels;
if (fSaveAuxDetInfo){
  evt.getView(fElecDriftModuleLabel, fAuxDetSimChannels);
}

std::vector<const sim::SimChannel*> fSimChannels;
if (fIsMC && fSaveGeantInfo){  evt.getView(fElecDriftModuleLabel, fSimChannels);}

  fData->run = evt.run();
  fData->subrun = evt.subRun();
  fData->event = evt.id().event() + fEventsPerSubrun*evt.subRun();

  art::Timestamp ts = evt.time();
  fData->evttime_seconds = ts.timeHigh();
  fData->evttime_nanoseconds = ts.timeLow();

  //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];
    }
  }

  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clockData);
//  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;


//RawDigit info
if (fSaveRawDigitInfo){
  fData->no_channels = (int) NChannels;
  if (NChannels>0)
  {
    fData->no_ticks = (int) rawdigitlist[0]->Samples();
    fData->no_ticksinallchannels = fData->no_channels*fData->no_ticks;
  }
  else
  {
  fData->no_ticks = 0;
  fData->no_ticksinallchannels = 0;
  }

  for (int i = 0; i < fData->no_channels ; i++) //loop over channels holding raw waveforms
  {
    fData->rawD_Channel[i] = (int) rawdigitlist[i]->Channel();
    int k=0;

    for (int j = fData->no_ticks*i; j < (fData->no_ticks*(i+1)) && j < kMaxReadoutTicksInAllChannels ; j++) //loop over ticks
    {
      fData->rawD_ADC[j] = rawdigitlist[i]->ADC(k);
      k++;
    }
  }
}

//RecobWire info
if (fSaveRecobWireInfo){
  fData->no_recochannels = (int) NRecoChannels;
  int RecoWTick=0;

  for(int i = 0; i < fData->no_recochannels; i++)  //loop over channels holding reco waveforms
  {
    fData->recoW_NTicks[i]=0;
    fData->recoW_Channel[i] = recobwirelist[i]->Channel();
    const recob::Wire::RegionsOfInterest_t& signalROI = recobwirelist[i]->SignalROI();

    for(const auto& range : signalROI.get_ranges())  //loop over ROIs
    {
      const std::vector<float>& signal = range.data();
      int NTicksInThisROI = range.end_index() - range.begin_index();

      for(int j = 0; j < NTicksInThisROI; j++) //loop over ticks
      {
        fData->recoW_Tick[RecoWTick] = j+range.begin_index();
        fData->recoW_ADC[RecoWTick] = signal.at(j);
        RecoWTick++;
      }
      fData->recoW_NTicks[i] += NTicksInThisROI;
    }
  }
  fData->no_recoticksinallchannels = RecoWTick;
}

//hit information
if (fSaveHitInfo){

  int trueID = -9999;
  float maxe = -9999;
  float purity = -9999;

  fData->no_hits = (int) NHits;
  fData->NHitsInAllTracks = (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("AnaRootParser:limits") << "event has " << NHits
      << " hits, only kMaxHits=" << kMaxHits << " stored in tree";
  }

// trying to assign a space point to those recob::Hits t hat were assigned to a track, checking for every hit. Not really working yet.
//    art::FindManyP<recob::SpacePoint> fmsp(hitlist,evt,"pmtrack");

  auto hitResults = anab::FVectorReader<recob::Hit, 4>::create(evt, "dprawhit");
  const auto & fitParams = hitResults->vectors();

  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_view[i]   = hitlist[i]->WireID().Plane;
    fData->hit_wire[i]    = hitlist[i]->WireID().Wire;
    fData->hit_peakT[i]   = hitlist[i]->PeakTime();
    fData->hit_chargesum[i]  = hitlist[i]->SummedADC();
    fData->hit_chargeintegral[i]  = hitlist[i]->Integral();
    fData->hit_ph[i]  = hitlist[i]->PeakAmplitude();
    fData->hit_startT[i] = hitlist[i]->StartTick();
    fData->hit_endT[i] = hitlist[i]->EndTick();
    fData->hit_rms[i] = hitlist[i]->RMS();
    fData->hit_goodnessOfFit[i] = hitlist[i]->GoodnessOfFit();

    fData->hit_fitparamt0[i] = fitParams[i][0];
    fData->hit_fitparamtau1[i] = fitParams[i][1];
    fData->hit_fitparamtau2[i] = fitParams[i][2];
    fData->hit_fitparamampl[i] = fitParams[i][3];

    fData->hit_multiplicity[i] = hitlist[i]->Multiplicity();

    //quantities from backtracker are different from the real value in MC

    if( fIsMC )
        HitsBackTrack(clockData, hitlist[i], trueID, maxe, purity );

    fData->hit_trueID[i] = trueID;
    fData->hit_trueEnergyMax[i] = maxe;
    fData->hit_trueEnergyFraction[i] = purity;

/*
    std::vector< art::Ptr<recob::SpacePoint> > sptv = fmsp.at(i);
        if (fmsp.at(i).size()!=0){
    std::cout << "sptv[i]->XYZ()[0]: " << sptv[i]->XYZ()[0] << std::endl;
    std::cout << "sptv[i]->XYZ()[1]: " << sptv[i]->XYZ()[1] << std::endl;
    std::cout << "sptv[i]->XYZ()[2]: " << sptv[i]->XYZ()[2] << std::endl;
      }
    std::cout << "test3" << std::endl;
//    std::cout << "test" << std::endl;
//    art::FindManyP<recob::SpacePoint> fmsp(hitListHandle,evt,"pmtrack");
//    art::FindManyP<recob::SpacePoint> fmsp(allHits, evt, fSpacePointModuleLabel);
//    art::FindManyP<recob::SpacePoint> fmsp(hitlist,evt,"pmtrack");
    for (size_t i = 0; i < NHits && i < kMaxHits ; ++i){//loop over hits
      if (fmsp.isValid()){
        if (fmsp.at(i).size()!=0){

        std::cout << "Spacepoint in x for hit" << i << ": " << fmsp.at(i)[0]->XYZ()[0] << std::endl;
        //  fData->hit_clusterid[i] = fmcl.at(i)[0]->ID();
          //        fData->hit_clusterKey[i] = fmcl.at(i)[0].key();
        }
        else std::cout << "No spacepoint for this hit" << std::endl;
      }
    }
    //Get space points associated with the hit
    std::vector< art::Ptr<recob::SpacePoint> > sptv = fmsp.at(hits[ipl][i]);
*/



    //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 (fIsMC){
    //        std::vector<sim::IDE> ides;
    //        try{bt_serv->HitToSimIDEs(hitlist[i], ides); }
    //        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 associated with hits
    art::FindManyP<recob::Cluster> fmcl(hitListHandle,evt,fClusterModuleLabel);
    for (size_t i = 0; i < NHits && i < kMaxHits ; ++i){//loop over hits
      if (fmcl.isValid()){
        if (fmcl.at(i).size()!=0){
          fData->hit_clusterid[i] = fmcl.at(i)[0]->ID();
          //        fData->hit_clusterKey[i] = fmcl.at(i)[0].key();
        }
        else
          fData->hit_clusterid[i] = -1;
      }
    }
  }
}// 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("AnaRootParser: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("AnaRootParser: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();
    // View 3 in LArSoft corresponds to View 0 in real life (3m long channels). View 2 in LArSoft corresponds to View 1 in real life (1m long channels).
    if(fData->clusterView[ic] == 3) {fData->clusterView[ic] = 0;}
    if(fData->clusterView[ic] == 2) {fData->clusterView[ic] = 1;}
    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();

    // Transform StartWire and EndWire to channels
    if(fData->clusterView[ic] == 1){fData->cluster_StartWire[ic]+=320; fData->cluster_EndWire[ic] += 320;}

    //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 (fSaveFlashInfo){
  fData->no_flashes = (int) NFlashes;
  if (NFlashes > kMaxFlashes) {
    // got this error? consider increasing kMaxHits
    // (or ask for a redesign using vectors)
    mf::LogError("AnaRootParser:limits") << "event has " << NFlashes
      << " flashes, only kMaxFlashes=" << kMaxFlashes << " stored in tree";
  }
  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("AnaRootParser: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){
  AnaRootParserDataStruct::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("AnaRootParser: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) {
    AnaRootParserDataStruct::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) {

  // Get Geometry
  art::ServiceHandle<geo::Geometry> geomhandle;

  for (unsigned int iTracker=0; iTracker < NTrackers; ++iTracker){
    AnaRootParserDataStruct::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;
    /*
    //First loop over tracks to get array sizes (number of hits in each track)
    for(size_t iTrk=0; iTrk < NTracks; ++iTrk){//loop over tracks
    art::FindMany<anab::Calorimetry> fmcaltemp(trackListHandle[iTracker], evt, fCalorimetryModuleLabel[iTracker]);
    if (fmcaltemp.isValid()){
    std::vector<const anab::Calorimetry*> calostemp = fmcaltemp.at(iTrk);
    for (size_t ical = 0; ical<calostemp.size(); ++ical){
    if (!calostemp[ical]) continue;
    if (!calostemp[ical]->PlaneID().isValid) continue;
    int planenumtemp = calostemp[ical]->PlaneID().Plane;
    if (planenumtemp<0||planenumtemp>2) continue;
    //For now make the second argument as 13 for muons.
    //      const size_t NHitsTemp = calostemp[ical] -> dEdx().size();
    //      TrackerData.NHitsInAllTracks+=(int) NHitsTemp;
    //      fData->NHitsInAllTracks+=(int) NHitsTemp;
    } // for calorimetry info
    } // if has calorimetry info
    }//loop over tracks
    */

    // 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("AnaRootParser: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
    // - change the minimal track length requirement to 50 cm
    trkf::TrackMomentumCalculator trkm{50.};

    recob::MCSFitResult res;

    int HitIterator=0;
    int HitIterator2=0;

    int trueID = -9999;
    float maxe = -9999.;
    float purity = -9999.;

    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, AnaRootParser initializes them with -999 or -9999
      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;
      TVector3 dir_start_flipped, dir_end_flipped;
      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>();

        dir_start_flipped.SetXYZ(dir_start.Z(), dir_start.Y(), dir_start.X());
        dir_end_flipped.SetXYZ(dir_end.Z(), dir_end.Y(), dir_end.X());

        tlen        = length(track);
        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.trkstarttheta[iTrk]         = (180.0/3.14159)*dir_start_flipped.Theta();
        TrackerData.trkstartphi[iTrk]           = (180.0/3.14159)*dir_start_flipped.Phi();
        TrackerData.trkstartdirectionx[iTrk]    = dir_start.X();
        TrackerData.trkstartdirectiony[iTrk]    = dir_start.Y();
        TrackerData.trkstartdirectionz[iTrk]    = dir_start.Z();

        if(fLogLevel >= 2)
        {
          std::cout << std::endl;
          std::cout << "start.X(): " << pos.X() << "\t" << "start.Y(): " << pos.Y() << "\t" << "start.Z(): " << pos.Z() << std::endl;
          std::cout << "end.X(): " << end.X() << "\t" << "end.Y(): " << end.Y() << "\t" << "end.Z(): " << end.Z() << std::endl;
          std::cout << "dir_start.X(): " << dir_start.X() << "\t" << "dir_start.Y(): " << dir_start.Y() << "\t" << "dir_start.Z(): " << dir_start.Z() << std::endl;
          std::cout << "dir_end.X(): " << dir_end.X() << "\t" << "dir_end.Y(): " << dir_end.Y() << "\t" << "dir_end.Z(): " << dir_end.Z() << std::endl;
          std::cout << "dir_start_flipped.Theta(): " << (180.0/3.14159)*dir_start_flipped.Theta() << "\t" << "dir_start_flipped.Phi(): " << (180.0/3.14159)*dir_start_flipped.Phi() << std::endl;
          std::cout << "dir_end_flipped.Theta(): " << (180.0/3.14159)*dir_end_flipped.Theta() << "\t" << "dir_end_flipped.Phi(): " << (180.0/3.14159)*dir_end_flipped.Phi() << std::endl;
          std::cout << std::endl;
        }

        TrackerData.trkendtheta[iTrk]     = (180.0/3.14159)*dir_end_flipped.Theta();
        TrackerData.trkendphi[iTrk]       = (180.0/3.14159)*dir_end_flipped.Phi();
        TrackerData.trkenddirectionx[iTrk]        = dir_end.X();
        TrackerData.trkenddirectiony[iTrk]        = dir_end.Y();
        TrackerData.trkenddirectionz[iTrk]        = dir_end.Z();

        TrackerData.trkthetaxz[iTrk]         = theta_xz;
        TrackerData.trkthetayz[iTrk]         = theta_yz;
        TrackerData.trkmom[iTrk]             = mom;
        TrackerData.trkchi2PerNDF[iTrk]      = track.Chi2PerNdof();
        TrackerData.trkNDF[iTrk]             = track.Ndof();
        TrackerData.trklen[iTrk]             = tlen;
        TrackerData.trklenstraightline[iTrk] = sqrt(pow(pos.X()-end.X(),2) + pow(pos.Y()-end.Y(),2) + pow(pos.Z()-end.Z(),2));
        TrackerData.trkmomrange[iTrk]        = trkm.GetTrackMomentum(tlen,13);
        TrackerData.trkmommschi2[iTrk]       = trkm.GetMomentumMultiScatterChi2(ptrack);
        TrackerData.trkmommsllhd[iTrk]       = trkm.GetMomentumMultiScatterLLHD(ptrack);

        //uBoone MCS
        res = fMCSFitter.fitMcs(*ptrack);
        TrackerData.trkmommscmic[iTrk]    = res.bestMomentum();
        TrackerData.trkmommscfwd[iTrk]    = res.fwdMomentum();
        TrackerData.trkmommscbwd[iTrk]    = res.bwdMomentum();
        TrackerData.trkmommscllfwd[iTrk]  = res.fwdLogLikelihood();
        TrackerData.trkmommscllbwd[iTrk]  = res.bwdLogLikelihood();


        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
      /*
      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("AnaRootParser: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::FindManyP<recob::Hit, recob::TrackHitMeta> fmthm(trackListHandle[iTracker], evt, "pmtrack");

//      if (fmthm.isValid()){
        auto vhit = fmthm.at(iTrk);
        auto vmeta = fmthm.data(iTrk);

        TrackerData.NHitsPerTrack[iTrk] = vhit.size();
        art::FindManyP<recob::SpacePoint> fmspts(vhit, evt, "pmtrack");

        int NHitsView0 = 0;
        int NHitsView1 = 0;

        if(fLogLevel >= 2)
        {
          std::cout << "track.NumberTrajectoryPoints(): " << track.NumberTrajectoryPoints() << std::endl;
          std::cout << "track.NPoints(): " << track.NPoints() << std::endl;
          std::cout << "vhit.size(): " << vhit.size() << std::endl;
          std::cout << "vmeta.size(): " << vmeta.size() << std::endl;
          std::cout << "fmspts.size(): " << fmspts.size() << std::endl;
        }

          for (unsigned int h = 0; h < vhit.size(); h++)
          {
            //corrected pitch
            double angleToVert = geomhandle->WireAngleToVertical(vhit[h]->View(), vhit[h]->WireID().TPC, vhit[h]->WireID().Cryostat) - 0.5*::util::pi<>();
            const TVector3& dir = tracklist[iTracker][iTrk]->DirectionAtPoint<TVector3>(h);
            const TVector3& loc = tracklist[iTracker][iTrk]->LocationAtPoint<TVector3>(h);
            double cosgamma = std::abs(std::sin(angleToVert)*dir.Y() + std::cos(angleToVert)*dir.Z());


            TrackerData.hittrklocaltrackdirectionx[HitIterator2] = dir.X();
            TrackerData.hittrklocaltrackdirectiony[HitIterator2] = dir.Y();
            TrackerData.hittrklocaltrackdirectionz[HitIterator2] = dir.Z();


            //XYZ
            std::vector< art::Ptr<recob::SpacePoint> > sptv = fmspts.at(h);
            TrackerData.hittrkx[HitIterator2] = sptv[0]->XYZ()[0];
            TrackerData.hittrky[HitIterator2] = sptv[0]->XYZ()[1];
            TrackerData.hittrkz[HitIterator2] = sptv[0]->XYZ()[2];

            TVector3 dir_hit_flipped;
            dir_hit_flipped.SetXYZ(dir.Z(), dir.Y(), dir.X());

            TrackerData.hittrklocaltrackdirectiontheta[HitIterator2] = (180.0/3.14159)*dir_hit_flipped.Theta();
            TrackerData.hittrklocaltrackdirectionphi[HitIterator2] = (180.0/3.14159)*dir_hit_flipped.Phi();

            //dx
            if(vhit[h]->WireID().Plane == 0) TrackerData.hittrkpitchC[HitIterator2] = std::abs(geomhandle->WirePitch()/( sin(dir_hit_flipped.Theta())*sin(dir_hit_flipped.Phi()) ));
            if(vhit[h]->WireID().Plane == 1) TrackerData.hittrkpitchC[HitIterator2] = std::abs(geomhandle->WirePitch()/( sin(dir_hit_flipped.Theta())*cos(dir_hit_flipped.Phi()) ));

            TrackerData.hittrkds[HitIterator2] = vmeta[h]->Dx();

            if(fLogLevel >= 2)
            {
              std::cout << "pos.X(): " << sptv[0]->XYZ()[0] << "\t" << "pos.Y(): " << sptv[0]->XYZ()[1] << "\t" << "pos.Z(): " << sptv[0]->XYZ()[2] << std::endl;
              std::cout << "pos2.X(): " << loc.X() << "\t" << "pos2.Y(): " << loc.Y() << "\t" << "pos2.Z(): " << loc.Z() << std::endl;
              std::cout << "dir.X(): " << dir.X() << "\t" << "dir.Y(): " << dir.Y() << "\t" << "dir.Z(): " << dir.Z() << std::endl;
              std::cout << "dir_hit_flipped.Theta(): " << (180.0/3.14159)*dir_hit_flipped.Theta() << "\t" << "dir_hit_flipped.Phi(): " << (180.0/3.14159)*dir_hit_flipped.Phi() << std::endl;
              std::cout << "vmeta[h]->Dx(): " << vmeta[h]->Dx() << std::endl;
              std::cout << "Dx corrected pitch old: " << geomhandle->WirePitch()/cosgamma << std::endl;
              std::cout << "Dx corrected pitch new: " << TrackerData.hittrkpitchC[HitIterator2] << std::endl;
              std::cout << "view: " << vhit[h]->WireID().Plane << std::endl;
            }

            //hit variables
            TrackerData.hittrkchannel[HitIterator2] = vhit[h]->Channel();
            TrackerData.hittrktpc[HitIterator2] = vhit[h]->WireID().TPC;
            TrackerData.hittrkview[HitIterator2] = vhit[h]->WireID().Plane;
            TrackerData.hittrkwire[HitIterator2] = vhit[h]->WireID().Wire;
            TrackerData.hittrkpeakT[HitIterator2] = vhit[h]->PeakTime();
            TrackerData.hittrkchargeintegral[HitIterator2] = vhit[h]->Integral();
            TrackerData.hittrkph[HitIterator2] = vhit[h]->PeakAmplitude();
            TrackerData.hittrkchargesum[HitIterator2] = vhit[h]->SummedADC();
            TrackerData.hittrkstarT[HitIterator2] = vhit[h]->StartTick();
            TrackerData.hittrkendT[HitIterator2] = vhit[h]->EndTick();
            TrackerData.hittrkrms[HitIterator2] = vhit[h]->RMS();
            TrackerData.hittrkgoddnessofFit[HitIterator2] = vhit[h]->GoodnessOfFit();
            TrackerData.hittrkmultiplicity[HitIterator2] = vhit[h]->Multiplicity();

      //quantities from backtracker are different from the real value in MC
      if( fIsMC )
              HitsBackTrack(clockData, vhit[h], trueID, maxe, purity );

      TrackerData.hittrktrueID[HitIterator2] = trueID;
      TrackerData.hittrktrueEnergyMax[HitIterator2] = maxe;
      TrackerData.hittrktrueEnergyFraction[HitIterator2] = purity;

            HitIterator2++;

            if(vhit[h]->WireID().Plane == 0) NHitsView0++;
            if(vhit[h]->WireID().Plane == 1) NHitsView1++;
          }
        TrackerData.ntrkhitsperview[iTrk][0] = NHitsView0;
        TrackerData.ntrkhitsperview[iTrk][1] = NHitsView1;

//      }

/*
        std::cout << "tracklist[iTracker][iTrk]->NumberTrajectoryPoints(): " << tracklist[iTracker][iTrk]->NumberTrajectoryPoints() << std::endl;
      for(size_t itp = 0; itp < tracklist[iTracker][iTrk]->NumberTrajectoryPoints(); ++itp)
      {
        const TVector3& pos = tracklist[iTracker][iTrk]->LocationAtPoint(itp);
      }
*/
//
//      art::FindManyP<recob::Hit> fmht(trackListHandle[iTracker], evt, fTrackModuleLabel);
//      std::vector< art::Ptr<recob::Hit> > allHits = fmht.at(trkIter);
//      art::FindManyP<recob::SpacePoint> fmspts(allHits, evt, fSpacePointModuleLabel);
//
/*
      art::FindManyP<recob::SpacePoint> fmspts(vhit, evt, "pmtrack");
        for (size_t h = 0; h < vhit.size(); ++h)
        {
          std::vector< art::Ptr<recob::SpacePoint> > sptv = fmspts.at(h);

          for (size_t j = 0; j < sptv.size(); ++j)
          {
            std::cout << "sptv[j]->XYZ()[0]: " << sptv[j]->XYZ()[0] << std::endl;
            std::cout << "sptv[j]->XYZ()[1]: " << sptv[j]->XYZ()[1] << std::endl;
            std::cout << "sptv[j]->XYZ()[2]: " << sptv[j]->XYZ()[2] << std::endl;
            std::cout << "sptv[j]->ErrXYZ()[0]: " << sptv[j]->ErrXYZ()[0] << std::endl;
            std::cout << "sptv[j]->ErrXYZ()[1]: " << sptv[j]->ErrXYZ()[1] << std::endl;
            std::cout << "sptv[j]->ErrXYZ()[2]: " << sptv[j]->ErrXYZ()[2] << std::endl;
          }
        }
*/
//


/////////////////////
      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("AnaRootParser: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){
        for (size_t ical = calos.size() - 1; ical <= 1 ; --ical){  //reverse order so that we access info on plane 0 (which is view 0 in real life) first
          if (!calos[ical]) continue;
          if (!calos[ical]->PlaneID().isValid) continue;
          int planenum = calos[ical]->PlaneID().Plane;
          if (planenum<0||planenum>1) continue; //only two planes (0 and 1) in dual phase
          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.ntrkhitsperview[iTrk][planenum] = (int) NHits;
          if (NHits > TrackerData.GetMaxHitsPerTrack(iTrk, planenum)) {
            // if you get this error, you'll have to increase kMaxTrackHits
            mf::LogError("AnaRootParser: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();
              //TrackerData.trkdedx2.push_back((calos[ical] -> dEdx())[iTrkHit]);
              TrackerData.trkdedx2[HitIterator] = (calos[ical] -> dEdx())[iTrkHit];
              TrackerData.trkdqdx2[HitIterator] = (calos[ical] -> dQdx())[iTrkHit];
              TrackerData.trktpc2[HitIterator] = (calos[ical] -> PlaneID()).TPC;
              TrackerData.trkx2[HitIterator] = TrkPos.X();
              TrackerData.trky2[HitIterator] = TrkPos.Y();
              TrackerData.trkz2[HitIterator] = TrkPos.Z();
              HitIterator++;
            } // for track hits
          }
        } // for calorimetry info


        //best plane
        if(TrackerData.ntrkhitsperview[iTrk][0] > TrackerData.ntrkhitsperview[iTrk][1] && TrackerData.ntrkhitsperview[iTrk][0] > TrackerData.ntrkhitsperview[iTrk][2]) TrackerData.trkpidbestplane[iTrk] = 0;
        else if(TrackerData.ntrkhitsperview[iTrk][1] > TrackerData.ntrkhitsperview[iTrk][0] && TrackerData.ntrkhitsperview[iTrk][1] > TrackerData.ntrkhitsperview[iTrk][2]) TrackerData.trkpidbestplane[iTrk] = 1;
        else if(TrackerData.ntrkhitsperview[iTrk][2] > TrackerData.ntrkhitsperview[iTrk][0] && TrackerData.ntrkhitsperview[iTrk][2] > TrackerData.ntrkhitsperview[iTrk][1]) TrackerData.trkpidbestplane[iTrk] = 2;
        else if(TrackerData.ntrkhitsperview[iTrk][2] == TrackerData.ntrkhitsperview[iTrk][0] && TrackerData.ntrkhitsperview[iTrk][2] > TrackerData.ntrkhitsperview[iTrk][1]) TrackerData.trkpidbestplane[iTrk] = 2;
        else if(TrackerData.ntrkhitsperview[iTrk][2] == TrackerData.ntrkhitsperview[iTrk][1] && TrackerData.ntrkhitsperview[iTrk][2] > TrackerData.ntrkhitsperview[iTrk][0]) TrackerData.trkpidbestplane[iTrk] = 2;
        else if(TrackerData.ntrkhitsperview[iTrk][1] == TrackerData.ntrkhitsperview[iTrk][0] && TrackerData.ntrkhitsperview[iTrk][1] > TrackerData.ntrkhitsperview[iTrk][2]) TrackerData.trkpidbestplane[iTrk] = 0;
        else if(TrackerData.ntrkhitsperview[iTrk][1] == TrackerData.ntrkhitsperview[iTrk][0] && TrackerData.ntrkhitsperview[iTrk][1] == TrackerData.ntrkhitsperview[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 (fIsMC){
        //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
             if( 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;
        HitsBackTrack(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;
        HitsBackTrack(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 (fIsMC)
*/
    }//end loop over track
      //std::cout << "HitIterator: " << HitIterator << std::endl;
  }//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){
      AnaRootParserDataStruct::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("AnaRootParser: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 (fIsMC)
  {

    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();
      }
    }

    //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("AnaRootParser: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){


      //Photon particles information
      if (fSavePhotonInfo)
      {
        int photoncounter=0;
        for ( auto const& pmt : (*photonHandle) )
        {
          std::map<int, int> PhotonsMap = pmt.DetectedPhotons;

          for(auto itphoton = PhotonsMap.begin(); itphoton!= PhotonsMap.end(); itphoton++)
          {
            for(int iphotonatthistime = 0; iphotonatthistime < itphoton->second ; iphotonatthistime++)
            {
              fData->photons_time[photoncounter]=itphoton->first;
              fData->photons_channel[photoncounter]=pmt.OpChannel;
              photoncounter++;
            }
          }
        }
        fData->numberofphotons=photoncounter;
      }

          //Generator particles information
      if (fSaveGeneratorInfo)
      {
        const sim::ParticleList& plist = pi_serv->ParticleList();

        size_t generator_particle=0;
        sim::ParticleList::const_iterator itPart = plist.begin(),
        pend = plist.end(); // iterator to pairs (track id, particle)

        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";
          }

          if (iPart < fData->GetMaxGeneratorparticles()) {

            std::string pri("primary");
            if( pPart->Mother() == 0 && pPart->Process() == pri )
            {
              fData->TrackId[generator_particle]=pPart->TrackId();
              fData->pdg[generator_particle]=pPart->PdgCode();
              fData->status[generator_particle] = pPart->StatusCode();
              fData->Eng[generator_particle]=pPart->E();
              fData->Mass[generator_particle]=pPart->Mass();
              fData->Px[generator_particle]=pPart->Px();
              fData->Py[generator_particle]=pPart->Py();
              fData->Pz[generator_particle]=pPart->Pz();
              fData->P[generator_particle]=pPart->Momentum().Vect().Mag();
              fData->StartPointx[generator_particle]=pPart->Vx();
              fData->StartPointy[generator_particle]=pPart->Vy();
              fData->StartPointz[generator_particle]=pPart->Vz();
              fData->StartT[generator_particle] = pPart->T();

              TVector3 momentum_start_flipped;
              momentum_start_flipped.SetXYZ(pPart->Pz(), pPart->Py(), pPart->Px());

              fData->theta[generator_particle] = (180.0/3.14159)*momentum_start_flipped.Theta();
              fData->phi[generator_particle] = (180.0/3.14159)*momentum_start_flipped.Phi();
            }
            ++generator_particle;
          }
          else if (iPart == fData->GetMaxGeneratorparticles()) {
            // got this error? it might be a bug,
            // since the structure should have enough room for everything
            mf::LogError("AnaRootParser:limits") << "event has "
              << plist.size() << " MC particles, only "
              << fData->GetMaxGeneratorparticles() << " will be stored in tree";
          }
        } // for particles
        fData->generator_list_size = generator_particle;
        fData->processname.resize(generator_particle);
        MF_LOG_DEBUG("AnaRootParser")
          << "Counted "
          << fData->generator_list_size << " Generator particles";
      }//fSaveGeneratorInfo




      //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";
          }

          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;
            unsigned int pstarti, pendi;
            double plen = length(*pPart, mcstart, mcend, pstarti, pendi);
            bool isActive = plen != 0;

//            TLorentzVector mcstartdrifted, mcenddrifted;
//            unsigned int pstartdriftedi, penddriftedi;
//            double plendrifted = driftedLength(*pPart, mcstartdrifted, mcenddrifted, pstartdriftedi, penddriftedi);
//            bool isDrifted = plendrifted!= 0;

            if(fLogLevel >= 3)
            {
              std::cout << "pPart->TrackId():" << pPart->TrackId() << std::endl;
              std::cout << "pPart->Mother():" << pPart->Mother() << std::endl;
              std::cout << "pPart->PdgCode():" << pPart->PdgCode() << std::endl;
              std::cout << std::endl;
            }


            fData->process_primary[geant_particle] = int(isPrimary);
            fData->processname[geant_particle]= pPart->Process();
            fData->Mother[geant_particle]=pPart->Mother();
            fData->TrackId[geant_particle]=pPart->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->NTrajectoryPointsPerParticle[geant_particle] = pPart->NumberTrajectoryPoints();

            //fData->theta[geant_particle] = pPart->Momentum().Theta();
            //fData->phi[geant_particle] = pPart->Momentum().Phi();
            //Change definition of theta and phi (swap x and z coordinate since x is "up" in dual phase)
            TVector3 momentum_start_flipped;
            momentum_start_flipped.SetXYZ(pPart->Pz(), pPart->Py(), pPart->Px());

            fData->theta[geant_particle] = (180.0/3.14159)*momentum_start_flipped.Theta();
            fData->phi[geant_particle] = (180.0/3.14159)*momentum_start_flipped.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_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 && fSaveGeantInAVInfo){
              fData->pathlen_tpcAV[active]  = plen;
              fData->TrackId_tpcAV[active] =  pPart->TrackId();
              fData->PDGCode_tpcAV[active] = pPart->PdgCode();

              fData->StartPointx_tpcAV[active] = mcstart.X();
              fData->StartPointy_tpcAV[active] = mcstart.Y();
              fData->StartPointz_tpcAV[active] = mcstart.Z();
              fData->StartT_tpcAV[active] = mcstart.T();
              fData->StartE_tpcAV[active] = pPart->E(pstarti);
              fData->StartP_tpcAV[active] = pPart->P(pstarti);
              fData->StartPx_tpcAV[active] = pPart->Px(pstarti);
              fData->StartPy_tpcAV[active] = pPart->Py(pstarti);
              fData->StartPz_tpcAV[active] = pPart->Pz(pstarti);
              fData->EndPointx_tpcAV[active] = mcend.X();
              fData->EndPointy_tpcAV[active] = mcend.Y();
              fData->EndPointz_tpcAV[active] = mcend.Z();
              fData->EndT_tpcAV[active] = mcend.T();
              fData->EndE_tpcAV[active] = pPart->E(pendi);
              fData->EndP_tpcAV[active] = pPart->P(pendi);
              fData->EndPx_tpcAV[active] = pPart->Px(pendi);
              fData->EndPy_tpcAV[active] = pPart->Py(pendi);
              fData->EndPz_tpcAV[active] = pPart->Pz(pendi);

              //Change definition of theta and phi (swap x and z coordinate since x is "up" in dual phase)
              TVector3 momentum_start_tpcAv_flipped;
              momentum_start_tpcAv_flipped.SetXYZ(pPart->Pz(pstarti), pPart->Py(pstarti), pPart->Px(pstarti));
              TVector3 momentum_end_tpcAv_flipped;
              momentum_end_tpcAv_flipped.SetXYZ(pPart->Pz(pendi), pPart->Py(pendi), pPart->Px(pendi));

              fData->thetastart_tpcAV[active] = (180.0/3.14159)*momentum_start_tpcAv_flipped.Theta();
              fData->phistart_tpcAV[active] = (180.0/3.14159)*momentum_start_tpcAv_flipped.Phi();

              fData->thetaend_tpcAV[active] = (180.0/3.14159)*momentum_end_tpcAv_flipped.Theta();
              fData->phiend_tpcAV[active] = (180.0/3.14159)*momentum_end_tpcAv_flipped.Phi();

              active++;
            }

/*
            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("AnaRootParser: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("AnaRootParser: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("AnaRootParser")
          << "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)


        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;
        unsigned int gindex = TrackIDtoIndex[currentMotherTrackId];
        if (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()){
        unsigned int igeantMother = TrackIDtoIndex[currentMotherTrackId];
        if (igeantMother<geant_particle){
        fData->MergedId[igeantMother] = currentMergedId;
        }
        }
        currentMotherTrackId = gmother[gindex];
        }
        ++currentMergedId;
        }// for merging check

      } // if (fSaveGeantInfo)

      //GEANT trajectory information
      if (fSaveGeantTrajectoryInfo)
      {
        const sim::ParticleList& plist = pi_serv->ParticleList();
        sim::ParticleList::const_iterator itPart = plist.begin(),
          pend = plist.end(); // iterator to pairs (track id, particle)
        int trajpointcounter = 0;

        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";
          }

          unsigned int numberTrajectoryPoints = pPart->NumberTrajectoryPoints();
          for(unsigned int trajpoint=0; trajpoint < numberTrajectoryPoints; ++trajpoint)
          {
            const TLorentzVector& trajpointPosition= pPart->Position(trajpoint);
            //const TLorentzVector& trajpointMomentum= pPart->Momentum(trajpoint);

            fData->TrajTrackId[trajpointcounter] = pPart->TrackId();
            fData->TrajPDGCode[trajpointcounter] = pPart->PdgCode();

            fData->TrajX[trajpointcounter] = trajpointPosition.X();
            fData->TrajY[trajpointcounter] = trajpointPosition.Y();
            fData->TrajZ[trajpointcounter] = trajpointPosition.Z();
            fData->TrajT[trajpointcounter] = pPart->T(trajpoint);
            fData->TrajE[trajpointcounter] = pPart->E(trajpoint);
            fData->TrajP[trajpointcounter] = pPart->P(trajpoint);
            fData->TrajPx[trajpointcounter] = pPart->Px(trajpoint);
            fData->TrajPy[trajpointcounter] = pPart->Py(trajpoint);
            fData->TrajPz[trajpointcounter] = pPart->Pz(trajpoint);

            TVector3 trajpointMomentum_flipped;
            trajpointMomentum_flipped.SetXYZ(pPart->Pz(trajpoint), pPart->Py(trajpoint), pPart->Px(trajpoint));

            fData->TrajTheta[trajpointcounter] = (180.0/3.14159)*trajpointMomentum_flipped.Theta();
            fData->TrajPhi[trajpointcounter] = (180.0/3.14159)*trajpointMomentum_flipped.Phi();

            if(fLogLevel >= 4)
            {
            std::cout << std::endl;
            std::cout << "trajpointcounter: " << trajpointcounter << std::endl;
            std::cout << "fData->TrajTrackId[trajpointcounter]: " << fData->TrajTrackId[trajpointcounter] << std::endl;
            std::cout << "fData->TrajPDGCode[trajpointcounter]: " << fData->TrajPDGCode[trajpointcounter] << std::endl;
            std::cout << "fData->TrajX[trajpointcounter]: " << fData->TrajX[trajpointcounter] << std::endl;
            std::cout << "fData->TrajY[trajpointcounter]: " << fData->TrajY[trajpointcounter] << std::endl;
            std::cout << "fData->TrajZ[trajpointcounter]: " << fData->TrajZ[trajpointcounter] << std::endl;
            std::cout << "fData->TrajT[trajpointcounter]: " << fData->TrajT[trajpointcounter] << std::endl;
            std::cout << "fData->TrajE[trajpointcounter]: " << fData->TrajE[trajpointcounter] << std::endl;
            std::cout << "fData->TrajP[trajpointcounter]: " << fData->TrajP[trajpointcounter] << std::endl;
            std::cout << "fData->TrajPx[trajpointcounter]: " << fData->TrajPx[trajpointcounter] << std::endl;
            std::cout << "fData->TrajPy[trajpointcounter]: " << fData->TrajPy[trajpointcounter] << std::endl;
            std::cout << "fData->TrajPz[trajpointcounter]: " << fData->TrajPz[trajpointcounter] << std::endl;
            std::cout << "fData->TrajTheta[trajpointcounter]: " << fData->TrajTheta[trajpointcounter] << std::endl;
            std::cout << "fData->TrajPhi[trajpointcounter]: " << fData->TrajPhi[trajpointcounter] << std::endl;
            }

            trajpointcounter++;
          }
        } // for particles

        fData->geant_trajectory_size = trajpointcounter;
      }// if (fSaveGeantTrajectoryInfo)


      if (fSaveSimEnergyDepositTPCActiveInfo)
      {
        fData->simenergydeposit_size = nSimEnergyDepositsTPCActive;
        fData->particleswithsimenergydeposit_size = nParticlesWithSimEnergyDepositsTPCActive;

        // Find particle ID's for each particle with sim energy deposit in this event
        int particlewithsedit=0;
        for(int sedavit = 0; sedavit < nSimEnergyDepositsTPCActive; sedavit++)
        {
          if (std::find(fData->ParticleIDSimEnergyDepositsTPCActivePerParticle.begin(), fData->ParticleIDSimEnergyDepositsTPCActivePerParticle.end(), energyDepositTPCActivelist[sedavit]->TrackID()) == fData->ParticleIDSimEnergyDepositsTPCActivePerParticle.end() ) //check if current particle ID is already in the vector. if true: not in vector
          {
             fData->ParticleIDSimEnergyDepositsTPCActivePerParticle[particlewithsedit] = energyDepositTPCActivelist[sedavit]->TrackID();

             particlewithsedit++;
          }
        }

        //Sort particle ID's for each particle with sim energy deposit in this event by particle ID.
        std::sort(fData->ParticleIDSimEnergyDepositsTPCActivePerParticle.begin(), fData->ParticleIDSimEnergyDepositsTPCActivePerParticle.end());


        // Find number of sim energy deposits for each particle in this event and sort all SEDs by Particle ID.
        std::vector<int> it_sortedbyparticleID;

        int totalsed = 0;
        for(int psedavit = 0; psedavit < nParticlesWithSimEnergyDepositsTPCActive; psedavit++)
        {
          int NSEDForThisParticle = 0;
          for(int sedavit = 0; sedavit < nSimEnergyDepositsTPCActive; sedavit++)
          {
            if(energyDepositTPCActivelist[sedavit]->TrackID() == fData->ParticleIDSimEnergyDepositsTPCActivePerParticle[psedavit])
            {
              NSEDForThisParticle++;
              it_sortedbyparticleID.push_back(sedavit);
            }
          }

          fData->NSimEnergyDepositsTPCActivePerParticle[psedavit] = NSEDForThisParticle;
          totalsed += NSEDForThisParticle;
        }

/*
        std::cout << std::endl;
        std::cout << "it_sortedbyparticleID.size(): " << it_sortedbyparticleID.size() << std::endl;

        for(int sedavit=0; sedavit < nSimEnergyDepositsTPCActive; sedavit++)
        {
          std::cout << std::endl;
          std::cout << "energyDepositTPCActivelist[" << sedavit << "]->TrackID(): " << energyDepositTPCActivelist[it_sortedbyparticleID[sedavit]]->TrackID() << std::endl;
          std::cout << "energyDepositTPCActivelist[" << sedavit << "]->Time(): " << energyDepositTPCActivelist[it_sortedbyparticleID[sedavit]]->Time() << std::endl;
        }
*/


//Algorithm to sort SEDs by time
/*
        //Now also sort SEDs by time
        //Sort by time (after sorting by particle ID)
        std::vector<int> it_sortedbyparticleIDandtime;

        int sedavit_currentparticle=0;
        int sedavit_mintime_currentparticle = -1;


        for(int psedavit = 0; psedavit < nParticlesWithSimEnergyDepositsTPCActive; psedavit++)
        {
          for(int NSEDForThisParticle = 0; NSEDForThisParticle < fData->NSimEnergyDepositsTPCActivePerParticle[psedavit]; NSEDForThisParticle++)
          {
            double mintime_currentparticle = 1e9;

            for(int sedavit = sedavit_currentparticle; sedavit < sedavit_currentparticle + fData->NSimEnergyDepositsTPCActivePerParticle[psedavit]; sedavit++)
            {
              if( energyDepositTPCActivelist[it_sortedbyparticleID[sedavit]]->Time() < mintime_currentparticle && ( std::find(it_sortedbyparticleIDandtime.begin(), it_sortedbyparticleIDandtime.end(), sedavit) == it_sortedbyparticleIDandtime.end() ) )
              {
                mintime_currentparticle = energyDepositTPCActivelist[it_sortedbyparticleID[sedavit]]->Time();
                sedavit_mintime_currentparticle = sedavit;
              }
            }
            it_sortedbyparticleIDandtime.push_back(sedavit_mintime_currentparticle);
          }
          sedavit_currentparticle+=fData->NSimEnergyDepositsTPCActivePerParticle[psedavit];
        }

*/


//Alternative algorithm to sort SEDs by time
/*
        //Sort by time (after sorting by particle ID)
        std::vector<int> it_sortedbyparticleIDandtime;
        it_sortedbyparticleIDandtime.resize(nSimEnergyDepositsTPCActive);
        FillWith(it_sortedbyparticleIDandtime, -1);

        int sedavit_sortedbyparticleIDandtime = 0;
        int mintime_sedavit = -1;
        int psedavit=0;

        //for(int sedavit_sortedbyparticleIDandtime = 0; sedavit_sortedbyparticleIDandtime < nSimEnergyDepositsTPCActive; i++)
        //{
        while(sedavit_sortedbyparticleIDandtime < nSimEnergyDepositsTPCActive )
        {
          double mintime = 1e9;
          bool FoundSEDForThisParticleID = false;
          int NSEDForThisParticle = 0;

          for(int sedavit = 0; sedavit < nSimEnergyDepositsTPCActive; sedavit++)
          {
            if( energyDepositTPCActivelist[sedavit]->TrackID() == fData->ParticleIDSimEnergyDepositsTPCActivePerParticle[psedavit] )
            {
              if( energyDepositTPCActivelist[sedavit]->Time() < mintime && ( std::find(it_sortedbyparticleIDandtime.begin(), it_sortedbyparticleIDandtime.end(), sedavit) == it_sortedbyparticleIDandtime.end() ) )
              {
                mintime = energyDepositTPCActivelist[sedavit]->Time();
                mintime_sedavit = sedavit;
                FoundSEDForThisParticleID = true;
              }
            }
          }

          if( FoundSEDForThisParticleID )
          {
            it_sortedbyparticleIDandtime[sedavit_sortedbyparticleIDandtime] = mintime_sedavit;
            sedavit_sortedbyparticleIDandtime++;
          }
          else
          {
            psedavit++;
          }
        }
*/

        // For each energy deposit in this event (sorted by particle ID, not by time)

        for(int sedavit = 0; sedavit < nSimEnergyDepositsTPCActive; sedavit++)
        {
          fData->SEDTPCAVTrackID[sedavit] = energyDepositTPCActivelist[it_sortedbyparticleID[sedavit]]->TrackID();
          fData->SEDTPCAVPDGCode[sedavit] = energyDepositTPCActivelist[it_sortedbyparticleID[sedavit]]->PdgCode();
          fData->SEDTPCAVEnergy[sedavit] = energyDepositTPCActivelist[it_sortedbyparticleID[sedavit]]->Energy();
          fData->SEDTPCAVNumPhotons[sedavit] = energyDepositTPCActivelist[it_sortedbyparticleID[sedavit]]->NumPhotons();
          fData->SEDTPCAVNumElectrons[sedavit] = energyDepositTPCActivelist[it_sortedbyparticleID[sedavit]]->NumElectrons();
          fData->SEDTPCAVLength[sedavit] = energyDepositTPCActivelist[it_sortedbyparticleID[sedavit]]->StepLength();

          fData->SEDTPCAVStartTime[sedavit] = energyDepositTPCActivelist[it_sortedbyparticleID[sedavit]]->StartT();
          fData->SEDTPCAVStartX[sedavit] = energyDepositTPCActivelist[it_sortedbyparticleID[sedavit]]->StartX();
          fData->SEDTPCAVStartY[sedavit] = energyDepositTPCActivelist[it_sortedbyparticleID[sedavit]]->StartY();
          fData->SEDTPCAVStartZ[sedavit] = energyDepositTPCActivelist[it_sortedbyparticleID[sedavit]]->StartZ();

          fData->SEDTPCAVMidTime[sedavit] = energyDepositTPCActivelist[it_sortedbyparticleID[sedavit]]->Time();
          fData->SEDTPCAVMidX[sedavit] = energyDepositTPCActivelist[it_sortedbyparticleID[sedavit]]->MidPointX();
          fData->SEDTPCAVMidY[sedavit] = energyDepositTPCActivelist[it_sortedbyparticleID[sedavit]]->MidPointY();
          fData->SEDTPCAVMidZ[sedavit] = energyDepositTPCActivelist[it_sortedbyparticleID[sedavit]]->MidPointZ();

          fData->SEDTPCAVEndTime[sedavit] = energyDepositTPCActivelist[it_sortedbyparticleID[sedavit]]->EndT();
          fData->SEDTPCAVEndX[sedavit] = energyDepositTPCActivelist[it_sortedbyparticleID[sedavit]]->EndX();
          fData->SEDTPCAVEndY[sedavit] = energyDepositTPCActivelist[it_sortedbyparticleID[sedavit]]->EndY();
          fData->SEDTPCAVEndZ[sedavit] = energyDepositTPCActivelist[it_sortedbyparticleID[sedavit]]->EndZ();


          if(fLogLevel >= 5)
          {
            std::cout << std::endl;
            std::cout << "sedavit: " << sedavit << std::endl;
            std::cout << "fData->SEDTPCAVTrackID[" << sedavit << "]: " << fData->SEDTPCAVTrackID[sedavit] << std::endl;
            std::cout << "fData->SEDTPCAVPDGCode[" << sedavit << "]: " << fData->SEDTPCAVPDGCode[sedavit] << std::endl;
            std::cout << "fData->SEDTPCAVEnergy[" << sedavit << "]: " << fData->SEDTPCAVEnergy[sedavit] << std::endl;
            std::cout << "fData->SEDTPCAVNumPhotons[" << sedavit << "]: " << fData->SEDTPCAVNumPhotons[sedavit] << std::endl;
            std::cout << "fData->SEDTPCAVNumElectrons[" << sedavit << "]: " << fData->SEDTPCAVNumElectrons[sedavit] << std::endl;

            std::cout << "fData->SEDTPCAVLength[" << sedavit << "]: " << fData->SEDTPCAVLength[sedavit] << std::endl;
            std::cout << "fData->SEDTPCAVStartTime[" << sedavit << "]: " << fData->SEDTPCAVStartTime[sedavit] << std::endl;
            std::cout << "fData->SEDTPCAVStartX[" << sedavit << "]: " << fData->SEDTPCAVStartX[sedavit] << std::endl;
            std::cout << "fData->SEDTPCAVStartY[" << sedavit << "]: " << fData->SEDTPCAVStartY[sedavit] << std::endl;
            std::cout << "fData->SEDTPCAVStartZ[" << sedavit << "]: " << fData->SEDTPCAVStartZ[sedavit] << std::endl;
            std::cout << "fData->SEDTPCAVMidTime[" << sedavit << "]: " << fData->SEDTPCAVMidTime[sedavit] << std::endl;
            std::cout << "fData->SEDTPCAVMidX[" << sedavit << "]: " << fData->SEDTPCAVMidX[sedavit] << std::endl;
            std::cout << "fData->SEDTPCAVMidY[" << sedavit << "]: " << fData->SEDTPCAVMidY[sedavit] << std::endl;
            std::cout << "fData->SEDTPCAVMidZ[" << sedavit << "]: " << fData->SEDTPCAVMidZ[sedavit] << std::endl;
            std::cout << "fData->SEDTPCAVEndTime[" << sedavit << "]: " << fData->SEDTPCAVEndTime[sedavit] << std::endl;
            std::cout << "fData->SEDTPCAVEndX[" << sedavit << "]: " << fData->SEDTPCAVEndX[sedavit] << std::endl;
            std::cout << "fData->SEDTPCAVEndY[" << sedavit << "]: " << fData->SEDTPCAVEndY[sedavit] << std::endl;
            std::cout << "fData->SEDTPCAVEndZ[" << sedavit << "]: " << fData->SEDTPCAVEndZ[sedavit] << std::endl;
          }
        }
      }
    }//if (mcevts_truth)
  }//if (fIsMC){

    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("AnaRootParserStructure");
      logStream
        << "Tree data structure contains:"
        << "\n - " << fData->no_hits << " hits (" << fData->GetMaxHits() << ")"
        << "\n - " << fData->NHitsInAllTracks << " 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->ntrkhitsperview[iTrk][0];
          for (size_t ipl = 1; ipl < tracker->GetMaxPlanesPerTrack(iTrk); ++ipl)
            logStream << " + " << tracker->ntrkhitsperview[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("AnaRootParserStructure") << "Freeing the tree data structure";
      DestroyData();
    }
    } // dune::AnaRootParser::analyze()

double dune::AnaRootParser::bdist ( const TVector3 &  pos )

private

Definition at line 6891 of file AnaRootParser_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 dune::AnaRootParser::beginSubRun

(

const art::SubRun &

sr

)

Definition at line 4173 of file AnaRootParser_module.cc.

{

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

}

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

inlineprivate

Helper function: throws if no data structure is available.

Definition at line 1594 of file AnaRootParser_module.cc.

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

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

inlineprivate

Helper function: throws if no tree is available.

Definition at line 1601 of file AnaRootParser_module.cc.

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

int dune::AnaRootParser::CountHits ( const art::Event &  evt, const art::InputTag &  which, unsigned int  cryostat, unsigned int  tpc, unsigned int  plane  )

private

Definition at line 7054 of file AnaRootParser_module.cc.

    {
      std::vector<const recob::Hit*> temp;
      int NumberOfHitsBeforeThisPlane=0;
      evt.getView(which, temp);   //temp.size() = total number of hits for this event (number of all hits in all Cryostats, TPC's, planes and wires)
      for(size_t t = 0; t < temp.size(); ++t){
        if( temp[t]->WireID().Cryostat == cryostat&& temp[t]->WireID().TPC == tpc && temp[t]->WireID().Plane == plane ) break;
        NumberOfHitsBeforeThisPlane++;
      }
      return NumberOfHitsBeforeThisPlane;
    }

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

inlineprivate

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

Definition at line 1494 of file AnaRootParser_module.cc.

      {
        if (!fData) {
          fData.reset
            (new AnaRootParserDataStruct(GetNTrackers(), GetNVertexAlgos(), GetShowerAlgos()));
          fData->SetBits(AnaRootParserDataStruct::tdCry,    !fSaveCryInfo);
          fData->SetBits(AnaRootParserDataStruct::tdGenie,  !fSaveGenieInfo);
          fData->SetBits(AnaRootParserDataStruct::tdProto,  !fSaveProtoInfo);
          fData->SetBits(AnaRootParserDataStruct::tdPhotons,  !fSavePhotonInfo);
          fData->SetBits(AnaRootParserDataStruct::tdGenerator,  !fSaveGeneratorInfo);
          fData->SetBits(AnaRootParserDataStruct::tdGeant,  !fSaveGeantInfo);
          fData->SetBits(AnaRootParserDataStruct::tdGeantInAV,  !fSaveGeantInAVInfo);
          fData->SetBits(AnaRootParserDataStruct::tdGeantTrajectory,  !fSaveGeantTrajectoryInfo);
          fData->SetBits(AnaRootParserDataStruct::tdSimEnergyDepositTPCActive,  !fSaveSimEnergyDepositTPCActiveInfo);
          fData->SetBits(AnaRootParserDataStruct::tdMCshwr, !fSaveMCShowerInfo);
          fData->SetBits(AnaRootParserDataStruct::tdMCtrk,  !fSaveMCTrackInfo);
          fData->SetBits(AnaRootParserDataStruct::tdHit,    !fSaveHitInfo);
          fData->SetBits(AnaRootParserDataStruct::tdRawDigit,    !fSaveRawDigitInfo);
          fData->SetBits(AnaRootParserDataStruct::tdRecobWire,    !fSaveRecobWireInfo);
          fData->SetBits(AnaRootParserDataStruct::tdFlash,  !fSaveFlashInfo);
          fData->SetBits(AnaRootParserDataStruct::tdCount,  !fSaveExternCounterInfo);
          fData->SetBits(AnaRootParserDataStruct::tdShower, !fSaveShowerInfo);
          fData->SetBits(AnaRootParserDataStruct::tdCluster,!fSaveClusterInfo);
          fData->SetBits(AnaRootParserDataStruct::tdPandoraNuVertex,!fSavePandoraNuVertexInfo);
          fData->SetBits(AnaRootParserDataStruct::tdTrack,  !fSaveTrackInfo);
          fData->SetBits(AnaRootParserDataStruct::tdVertex, !fSaveVertexInfo);
          fData->SetBits(AnaRootParserDataStruct::tdAuxDet, !fSaveAuxDetInfo);
          fData->SetBits(AnaRootParserDataStruct::tdPFParticle, !fSavePFParticleInfo);
        }
        else {
          fData->SetTrackers(GetNTrackers());
          fData->SetVertexAlgos(GetNVertexAlgos());
          fData->SetShowerAlgos(GetShowerAlgos());
          if (bClearData) fData->Clear();
        }
      } // CreateData()

void dune::AnaRootParser::CreateTree ( bool  bClearData = false )

private

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

Definition at line 4155 of file AnaRootParser_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::AnaRootParser::CreateTree()

void dune::AnaRootParser::DestroyData ( )

inlineprivate

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

Definition at line 1591 of file AnaRootParser_module.cc.

{ fData.reset(); }

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

private

Definition at line 6976 of file AnaRootParser_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 dune::AnaRootParser::driftedLength ( detinfo::DetectorPropertiesData const &  detProp, const sim::MCTrack &  mctrack, TLorentzVector &  tpcstart, TLorentzVector &  tpcend, TLorentzVector &  tpcmom  )

private

Definition at line 6925 of file AnaRootParser_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 dune::AnaRootParser::endSubRun

(

const art::SubRun &

sr

)

Definition at line 4184 of file AnaRootParser_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 dune::AnaRootParser::FillShower ( AnaRootParserDataStruct::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 6772 of file AnaRootParser_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::AnaRootParser::FillShower()

void dune::AnaRootParser::FillShowers ( AnaRootParserDataStruct::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 6815 of file AnaRootParser_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("AnaRootParser: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::AnaRootParser::FillShowers()

size_t dune::AnaRootParser::GetNShowerAlgos ( ) const

inlineprivate

Returns the number of shower algorithms configured.

Definition at line 1487 of file AnaRootParser_module.cc.

{ return fShowerModuleLabel.size(); }

size_t dune::AnaRootParser::GetNTrackers ( ) const

inlineprivate

Returns the number of trackers configured.

Definition at line 1482 of file AnaRootParser_module.cc.

{ return fTrackModuleLabel.size(); }

size_t dune::AnaRootParser::GetNVertexAlgos ( ) const

inlineprivate

Definition at line 1484 of file AnaRootParser_module.cc.

{ return fVertexModuleLabel.size(); }

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

inlineprivate

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

Definition at line 1490 of file AnaRootParser_module.cc.

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

void dune::AnaRootParser::HitsBackTrack ( detinfo::DetectorClocksData const &  clockData, art::Ptr< recob::Hit > const &  hit, int &  trackid, float &  maxe, float &  purity  )

private

Definition at line 6856 of file AnaRootParser_module.cc.

                                                                                                                                                         {

      trackid = -1;
      purity = -1;

      art::ServiceHandle<cheat::BackTrackerService> bt_serv;

      std::map<int,double> trkide;
      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<<" "<<e<<" "<<eveIDs[e].trackID<<" "<<eveIDs[e].energy<<" "<<eveIDs[e].energyFrac<<std::endl;
          trkide[eveIDs[e].trackID] += eveIDs[e].energy;
      }

      maxe = 0;
      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 dune::AnaRootParser::length ( const recob::Track &  track )

private

Definition at line 6919 of file AnaRootParser_module.cc.

    {
      return track.Length();
    }

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

private

Definition at line 7028 of file AnaRootParser_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::AnaRootParser::SetAddresses ( )

inlineprivate

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

Definition at line 1532 of file AnaRootParser_module.cc.

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

void dune::AnaRootParser::SetPFParticleAddress ( )

inlineprivate

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

Definition at line 1581 of file AnaRootParser_module.cc.

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

void dune::AnaRootParser::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 1568 of file AnaRootParser_module.cc.

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

void dune::AnaRootParser::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 1541 of file AnaRootParser_module.cc.

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

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

inlineprivate

Definition at line 1554 of file AnaRootParser_module.cc.

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

Member Data Documentation

double dune::AnaRootParser::ActiveBounds[6]

private

objet holding the configuration of the uBoone MCS fitting alg

Definition at line 1479 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::bIgnoreMissingShowers

private

whether to ignore missing shower information

Definition at line 1471 of file AnaRootParser_module.cc.

std::string dune::AnaRootParser::fCalDataModuleLabel

private

Definition at line 1415 of file AnaRootParser_module.cc.

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

private

Definition at line 1431 of file AnaRootParser_module.cc.

std::string dune::AnaRootParser::fClusterModuleLabel

private

Definition at line 1421 of file AnaRootParser_module.cc.

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

private

Definition at line 1468 of file AnaRootParser_module.cc.

std::string dune::AnaRootParser::fCosmicClusterTaggerAssocLabel

private

Definition at line 1438 of file AnaRootParser_module.cc.

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

private

whether to extract and save PFParticle information

Definition at line 1467 of file AnaRootParser_module.cc.

std::string dune::AnaRootParser::fCryGenModuleLabel

private

Definition at line 1417 of file AnaRootParser_module.cc.

std::unique_ptr<AnaRootParserDataStruct> dune::AnaRootParser::fData

private

Definition at line 1404 of file AnaRootParser_module.cc.

std::string dune::AnaRootParser::fElecDriftModuleLabel

private

Definition at line 1414 of file AnaRootParser_module.cc.

short dune::AnaRootParser::fEventsPerSubrun

private

Definition at line 1409 of file AnaRootParser_module.cc.

std::string dune::AnaRootParser::fExternalCounterModuleLabel

private

Definition at line 1424 of file AnaRootParser_module.cc.

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

private

Definition at line 1469 of file AnaRootParser_module.cc.

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

private

Definition at line 1435 of file AnaRootParser_module.cc.

float dune::AnaRootParser::fG4minE

private

Energy threshold to save g4 particle info.

Definition at line 1475 of file AnaRootParser_module.cc.

std::string dune::AnaRootParser::fG4ModuleLabel

private

Definition at line 1419 of file AnaRootParser_module.cc.

std::string dune::AnaRootParser::fGenieGenModuleLabel

private

Definition at line 1416 of file AnaRootParser_module.cc.

std::string dune::AnaRootParser::fHitsModuleLabel

private

Definition at line 1412 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fIsMC

private

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

Definition at line 1439 of file AnaRootParser_module.cc.

int dune::AnaRootParser::fLogLevel

private

Definition at line 1408 of file AnaRootParser_module.cc.

trkf::TrajectoryMCSFitter dune::AnaRootParser::fMCSFitter

private

Definition at line 1477 of file AnaRootParser_module.cc.

std::string dune::AnaRootParser::fMCShowerModuleLabel

private

Definition at line 1425 of file AnaRootParser_module.cc.

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

private

Definition at line 1436 of file AnaRootParser_module.cc.

std::string dune::AnaRootParser::fMCTrackModuleLabel

private

Definition at line 1426 of file AnaRootParser_module.cc.

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

private

Definition at line 1433 of file AnaRootParser_module.cc.

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

private

Definition at line 1434 of file AnaRootParser_module.cc.

std::string dune::AnaRootParser::fOpFlashModuleLabel

private

Definition at line 1423 of file AnaRootParser_module.cc.

std::string dune::AnaRootParser::fPandoraNuVertexModuleLabel

private

Definition at line 1422 of file AnaRootParser_module.cc.

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

private

Definition at line 1432 of file AnaRootParser_module.cc.

std::string dune::AnaRootParser::fPFParticleModuleLabel

private

Definition at line 1428 of file AnaRootParser_module.cc.

std::string dune::AnaRootParser::fPhotonPropS1ModuleLabel

private

Definition at line 1413 of file AnaRootParser_module.cc.

std::string dune::AnaRootParser::fPOTModuleLabel

private

Definition at line 1437 of file AnaRootParser_module.cc.

std::string dune::AnaRootParser::fProtoGenModuleLabel

private

Definition at line 1418 of file AnaRootParser_module.cc.

std::string dune::AnaRootParser::fRawDigitModuleLabel

private

Definition at line 1411 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSaveAuxDetInfo

private

whether to extract and save auxiliary detector data

Definition at line 1443 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSaveCaloCosmics

private

save calorimetry information for cosmics

Definition at line 1474 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSaveClusterInfo

private

whether to extract and save Vertex information

Definition at line 1460 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSaveCryInfo

private

Definition at line 1444 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSaveExternCounterInfo

private

whether to extract and save Flash information

Definition at line 1463 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSaveFlashInfo

private

whether to extract and save nu vertex information from Pandora

Definition at line 1462 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSaveGeantInAVInfo

private

whether to extract and save Geant information

Definition at line 1450 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSaveGeantInfo

private

whether to extract and save Geant information

Definition at line 1449 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSaveGeantTrajectoryInfo

private

whether to extract and save Geant information

Definition at line 1451 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSaveGeneratorInfo

private

whether to extract and save collected photons

Definition at line 1448 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSaveGenieInfo

private

whether to extract and save CRY particle data

Definition at line 1445 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSaveHitInfo

private

whether to extract and save MC Track information

Definition at line 1455 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSaveMCShowerInfo

private

whether to extract and save Cluster information

Definition at line 1453 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSaveMCTrackInfo

private

whether to extract and save MC Shower information

Definition at line 1454 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSavePandoraNuVertexInfo

private

whether to extract and save Cluster information

Definition at line 1461 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSavePFParticleInfo

private

whether to extract and save Shower information

Definition at line 1465 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSavePhotonInfo

private

whether to extract and save ProtDUNE beam simulation information

Definition at line 1447 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSaveProtoInfo

private

whether to extract and save Genie information

Definition at line 1446 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSaveRawDigitInfo

private

whether to extract and save Hit information

Definition at line 1456 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSaveRecobWireInfo

private

whether to extract and save raw digit information

Definition at line 1457 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSaveShowerInfo

private

whether to extract and save External Counter information

Definition at line 1464 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSaveSimEnergyDepositTPCActiveInfo

private

whether to extract and save Geant information

Definition at line 1452 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSaveTrackInfo

private

whether to extract and save recob wire information

Definition at line 1458 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fSaveVertexInfo

private

whether to extract and save Track information

Definition at line 1459 of file AnaRootParser_module.cc.

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

private

Definition at line 1430 of file AnaRootParser_module.cc.

std::string dune::AnaRootParser::fSimEnergyDepositTPCActiveLabel

private

Definition at line 1420 of file AnaRootParser_module.cc.

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

private

Definition at line 1427 of file AnaRootParser_module.cc.

TTree* dune::AnaRootParser::fTree

private

Definition at line 1399 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::fUseBuffer

private

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

Definition at line 1440 of file AnaRootParser_module.cc.

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

private

Definition at line 1429 of file AnaRootParser_module.cc.

bool dune::AnaRootParser::isCosmics

private

if it contains cosmics

Definition at line 1473 of file AnaRootParser_module.cc.

AnaRootParserDataStruct::SubRunData_t dune::AnaRootParser::SubRunData

private

Definition at line 1406 of file AnaRootParser_module.cc.

---

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

• duneprototypes/duneprototypes/Protodune/dualphase/AnaRootParser_module.cc