Public Member Functions | Private Member Functions | Private Attributes | List of all members
dune::AnalysisTree Class Reference

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

Inheritance diagram for dune::AnalysisTree:
art::EDAnalyzer art::detail::Analyzer art::detail::LegacyModule art::Observer art::ModuleBase

Public Member Functions

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

Private Member Functions

void HitsPurity (detinfo::DetectorClocksData const &clockData, std::vector< art::Ptr< recob::Hit > > const &hits, Int_t &trackid, Float_t &purity, double &maxe)
 
double length (const recob::Track &track)
 
double driftedLength (detinfo::DetectorPropertiesData const &detProp, const simb::MCParticle &part, TLorentzVector &start, TLorentzVector &end, unsigned int &starti, unsigned int &endi)
 
double driftedLength (detinfo::DetectorPropertiesData const &detProp, const sim::MCTrack &mctrack, TLorentzVector &tpcstart, TLorentzVector &tpcend, TLorentzVector &tpcmom)
 
double length (const simb::MCParticle &part, TLorentzVector &start, TLorentzVector &end, unsigned int &starti, unsigned int &endi)
 
double bdist (const TVector3 &pos)
 
size_t GetNTrackers () const
 Returns the number of trackers configured. More...
 
size_t GetNVertexAlgos () const
 
size_t GetNShowerAlgos () const
 Returns the number of shower algorithms configured. More...
 
std::vector< std::stringGetShowerAlgos () const
 Returns the name of configured shower algorithms (converted to string) More...
 
void CreateData (bool bClearData=false)
 Creates the structure for the tree data; optionally initializes it. More...
 
void SetAddresses ()
 Sets the addresses of all the tree branches, creating the missing ones. More...
 
void SetTrackerAddresses (size_t iTracker)
 
void SetVertexAddresses (size_t iVertexAlg)
 
void SetShowerAddresses (size_t iShower)
 
void SetPFParticleAddress ()
 
void CreateTree (bool bClearData=false)
 Create the output tree and the data structures, if needed. More...
 
void DestroyData ()
 Destroy the local buffers (existing branches will point to invalid address!) More...
 
void CheckData (std::string caller) const
 Helper function: throws if no data structure is available. More...
 
void CheckTree (std::string caller) const
 Helper function: throws if no tree is available. More...
 
void FillShower (AnalysisTreeDataStruct::ShowerDataStruct &showerData, size_t iShower, recob::Shower const &showers, const bool fSavePFParticleInfo, const std::map< Short_t, Short_t > &showerIDtoPFParticleIDMap) const
 Stores the information of shower in slot iShower of showerData. More...
 
void FillShowers (AnalysisTreeDataStruct::ShowerDataStruct &showerData, std::vector< recob::Shower > const &showers, const bool fSavePFParticleInfo, const std::map< Short_t, Short_t > &showerIDtoPFParticleIDMap) const
 Stores the information of all showers into showerData. More...
 

Private Attributes

TTree * fTree
 
TTree * fPOT
 
std::unique_ptr< AnalysisTreeDataStructfData
 
AnalysisTreeDataStruct::SubRunData_t SubRunData
 
std::string fDigitModuleLabel
 
std::string fHitsModuleLabel
 
std::string fLArG4ModuleLabel
 
std::string fCalDataModuleLabel
 
std::string fGenieGenModuleLabel
 
std::string fCryGenModuleLabel
 
std::string fProtoGenModuleLabel
 
std::string fG4ModuleLabel
 
std::string fClusterModuleLabel
 
std::string fPandoraNuVertexModuleLabel
 
std::string fOpFlashModuleLabel
 
std::string fExternalCounterModuleLabel
 
std::string fMCShowerModuleLabel
 
std::string fMCTrackModuleLabel
 
std::string fSpacePointSolverModuleLabel
 
std::string fCnnModuleLabel
 
std::vector< std::stringfTrackModuleLabel
 
std::string fPFParticleModuleLabel
 
std::vector< std::stringfVertexModuleLabel
 
std::vector< std::stringfShowerModuleLabel
 
std::vector< std::stringfCalorimetryModuleLabel
 
std::vector< std::stringfParticleIDModuleLabel
 
std::vector< std::stringfMVAPIDShowerModuleLabel
 
std::vector< std::stringfMVAPIDTrackModuleLabel
 
std::vector< std::stringfFlashT0FinderLabel
 
std::vector< std::stringfMCT0FinderLabel
 
std::string fPOTModuleLabel
 
std::string fCosmicClusterTaggerAssocLabel
 
bool fUseBuffer
 whether to use a permanent buffer (faster, huge memory) More...
 
bool fSaveAuxDetInfo
 whether to extract and save auxiliary detector data More...
 
bool fSaveCryInfo
 
bool fSaveGenieInfo
 whether to extract and save CRY particle data More...
 
bool fSaveProtoInfo
 whether to extract and save Genie information More...
 
bool fSaveGeantInfo
 whether to extract and save ProtDUNE beam simulation information More...
 
bool fSaveMCShowerInfo
 whether to extract and save Geant information More...
 
bool fSaveMCTrackInfo
 whether to extract and save MC Shower information More...
 
bool fSaveHitInfo
 whether to extract and save MC Track information More...
 
bool fSaveRawDigitInfo
 whether to extract and save Hit information More...
 
bool fSaveTrackInfo
 whether to extract and save Raw Digit information More...
 
bool fSaveVertexInfo
 whether to extract and save Track information More...
 
bool fSaveClusterInfo
 whether to extract and save Vertex information More...
 
bool fSavePandoraNuVertexInfo
 whether to extract and save Cluster information More...
 
bool fSaveFlashInfo
 whether to extract and save nu vertex information from Pandora More...
 
bool fSaveExternCounterInfo
 whether to extract and save Flash information More...
 
bool fSaveShowerInfo
 whether to extract and save External Counter information More...
 
bool fSavePFParticleInfo
 whether to extract and save Shower information More...
 
bool fSaveSpacePointSolverInfo
 whether to extract and save PFParticle information More...
 
bool fSaveCnnInfo
 whether to extract and save SpacePointSolver information More...
 
std::vector< std::stringfCosmicTaggerAssocLabel
 whether to extract and save CNN information More...
 
std::vector< std::stringfContainmentTaggerAssocLabel
 
std::vector< std::stringfFlashMatchAssocLabel
 
bool bIgnoreMissingShowers
 whether to ignore missing shower information More...
 
bool isCosmics
 if it contains cosmics More...
 
bool fSaveCaloCosmics
 save calorimetry information for cosmics More...
 
float fG4minE
 Energy threshold to save g4 particle info. More...
 
double ActiveBounds [6]
 

Additional Inherited Members

- Public Types inherited from art::EDAnalyzer
using WorkerType = WorkerT< EDAnalyzer >
 
using ModuleType = EDAnalyzer
 
- Protected Member Functions inherited from art::Observer
std::string const & processName () const
 
bool wantAllEvents () const noexcept
 
bool wantEvent (ScheduleID id, Event const &e) const
 
Handle< TriggerResultsgetTriggerResults (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
ConsumesCollectorconsumesCollector ()
 
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

Definition at line 1213 of file AnalysisTree_module.cc.

Constructor & Destructor Documentation

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

Definition at line 3326 of file AnalysisTree_module.cc.

3326  :
3327  EDAnalyzer(pset),
3328  fTree(nullptr), fPOT(nullptr),
3329  fDigitModuleLabel (pset.get< std::string >("DigitModuleLabel") ),
3330  fHitsModuleLabel (pset.get< std::string >("HitsModuleLabel") ),
3331  fLArG4ModuleLabel (pset.get< std::string >("LArGeantModuleLabel") ),
3332  fCalDataModuleLabel (pset.get< std::string >("CalDataModuleLabel") ),
3333  fGenieGenModuleLabel (pset.get< std::string >("GenieGenModuleLabel") ),
3334  fCryGenModuleLabel (pset.get< std::string >("CryGenModuleLabel") ),
3335  fProtoGenModuleLabel (pset.get< std::string >("ProtoGenModuleLabel") ),
3336  fG4ModuleLabel (pset.get< std::string >("G4ModuleLabel") ),
3337  fClusterModuleLabel (pset.get< std::string >("ClusterModuleLabel") ),
3338  fPandoraNuVertexModuleLabel (pset.get< std::string >("PandoraNuVertexModuleLabel") ),
3339  fOpFlashModuleLabel (pset.get< std::string >("OpFlashModuleLabel") ),
3340  fExternalCounterModuleLabel (pset.get< std::string >("ExternalCounterModuleLabel") ),
3341  fMCShowerModuleLabel (pset.get< std::string >("MCShowerModuleLabel") ),
3342  fMCTrackModuleLabel (pset.get< std::string >("MCTrackModuleLabel") ),
3343  fSpacePointSolverModuleLabel (pset.get< std::string >("SpacePointSolverModuleLabel")),
3344  fCnnModuleLabel (pset.get< std::string >("CnnModuleLabel")),
3345  fTrackModuleLabel (pset.get< std::vector<std::string> >("TrackModuleLabel")),
3346  fVertexModuleLabel (pset.get< std::vector<std::string> >("VertexModuleLabel")),
3347  fShowerModuleLabel (pset.get< std::vector<std::string> >("ShowerModuleLabel")),
3348  fCalorimetryModuleLabel (pset.get< std::vector<std::string> >("CalorimetryModuleLabel")),
3349  fParticleIDModuleLabel (pset.get< std::vector<std::string> >("ParticleIDModuleLabel") ),
3350  fMVAPIDShowerModuleLabel (pset.get< std::vector<std::string> >("MVAPIDShowerModuleLabel") ),
3351  fMVAPIDTrackModuleLabel (pset.get< std::vector<std::string> >("MVAPIDTrackModuleLabel") ),
3352  fFlashT0FinderLabel (pset.get< std::vector<std::string> >("FlashT0FinderLabel") ),
3353  fMCT0FinderLabel (pset.get< std::vector<std::string> >("MCT0FinderLabel") ),
3354  fPOTModuleLabel (pset.get< std::string >("POTModuleLabel")),
3355  fCosmicClusterTaggerAssocLabel (pset.get< std::string >("CosmicClusterTaggerAssocLabel")),
3356  fUseBuffer (pset.get< bool >("UseBuffers", false)),
3357  fSaveAuxDetInfo (pset.get< bool >("SaveAuxDetInfo", false)),
3358  fSaveCryInfo (pset.get< bool >("SaveCryInfo", false)),
3359  fSaveGenieInfo (pset.get< bool >("SaveGenieInfo", false)),
3360  fSaveProtoInfo (pset.get< bool >("SaveProtoInfo", false)),
3361  fSaveGeantInfo (pset.get< bool >("SaveGeantInfo", false)),
3362  fSaveMCShowerInfo (pset.get< bool >("SaveMCShowerInfo", false)),
3363  fSaveMCTrackInfo (pset.get< bool >("SaveMCTrackInfo", false)),
3364  fSaveHitInfo (pset.get< bool >("SaveHitInfo", false)),
3365  fSaveRawDigitInfo (pset.get< bool >("SaveRawDigitInfo", false)),
3366  fSaveTrackInfo (pset.get< bool >("SaveTrackInfo", false)),
3367  fSaveVertexInfo (pset.get< bool >("SaveVertexInfo", false)),
3368  fSaveClusterInfo (pset.get< bool >("SaveClusterInfo", false)),
3369  fSavePandoraNuVertexInfo (pset.get< bool >("SavePandoraNuVertexInfo", false)),
3370  fSaveFlashInfo (pset.get< bool >("SaveFlashInfo", false)),
3371  fSaveExternCounterInfo (pset.get< bool >("SaveExternCounterInfo", false)),
3372  fSaveShowerInfo (pset.get< bool >("SaveShowerInfo", false)),
3373  fSavePFParticleInfo (pset.get< bool >("SavePFParticleInfo", false)),
3374  fSaveSpacePointSolverInfo (pset.get< bool >("SaveSpacePointSolverInfo", false)),
3375  fSaveCnnInfo (pset.get< bool >("SaveCnnInfo", false)),
3376  fCosmicTaggerAssocLabel (pset.get<std::vector< std::string > >("CosmicTaggerAssocLabel") ),
3377  fContainmentTaggerAssocLabel (pset.get<std::vector< std::string > >("ContainmentTaggerAssocLabel") ),
3378  fFlashMatchAssocLabel (pset.get<std::vector< std::string > >("FlashMatchAssocLabel") ),
3379  bIgnoreMissingShowers (pset.get< bool >("IgnoreMissingShowers", false)),
3380  isCosmics(false),
3381  fSaveCaloCosmics (pset.get< bool >("SaveCaloCosmics",false)),
3382  fG4minE (pset.get< float>("G4minE",0.01))
3383 {
3384 
3385  if (fSavePFParticleInfo) fPFParticleModuleLabel = pset.get<std::string>("PFParticleModuleLabel");
3386 
3387  if (fSaveAuxDetInfo == true) fSaveGeantInfo = true;
3388  if (fSaveRawDigitInfo == true) fSaveHitInfo = true;
3389  mf::LogInfo("AnalysisTree") << "Configuration:"
3390  << "\n UseBuffers: " << std::boolalpha << fUseBuffer
3391  ;
3392  if (GetNTrackers() > kMaxTrackers) {
3394  << "AnalysisTree currently supports only up to " << kMaxTrackers
3395  << " tracking algorithms, but " << GetNTrackers() << " are specified."
3396  << "\nYou can increase kMaxTrackers and recompile.";
3397  } // if too many trackers
3398  if (fTrackModuleLabel.size() != fCalorimetryModuleLabel.size()){
3400  << "fTrackModuleLabel.size() = "<<fTrackModuleLabel.size()<<" does not match "
3401  << "fCalorimetryModuleLabel.size() = "<<fCalorimetryModuleLabel.size();
3402  }
3403  if (fTrackModuleLabel.size() != fParticleIDModuleLabel.size()){
3405  << "fTrackModuleLabel.size() = "<<fTrackModuleLabel.size()<<" does not match "
3406  << "fParticleIDModuleLabel.size() = "<<fParticleIDModuleLabel.size();
3407  }
3408  if (fTrackModuleLabel.size() != fFlashT0FinderLabel.size()){
3410  << "fTrackModuleLabel.size() = "<<fTrackModuleLabel.size()<<" does not match "
3411  << "fFlashT0FinderLabel.size() = "<<fFlashT0FinderLabel.size();
3412  }
3413  if (fTrackModuleLabel.size() != fMCT0FinderLabel.size()){
3415  << "fTrackModuleLabel.size() = "<<fTrackModuleLabel.size()<<" does not match "
3416  << "fMCT0FinderLabel.size() = "<<fMCT0FinderLabel.size();
3417  }
3418  if (fTrackModuleLabel.size() != fCosmicTaggerAssocLabel.size()) {
3420  << "fTrackModuleLabel.size() = "<<fTrackModuleLabel.size()<<" does not match "
3421  << "fCosmicTaggerAssocLabel.size() = "<<fCosmicTaggerAssocLabel.size();
3422  }
3423  if (fTrackModuleLabel.size() != fContainmentTaggerAssocLabel.size()) {
3425  << "fTrackModuleLabel.size() = "<<fTrackModuleLabel.size()<<" does not match "
3426  << "fCosmicTaggerAssocLabel.size() = "<<fContainmentTaggerAssocLabel.size();
3427  }
3428  if (fTrackModuleLabel.size() != fFlashMatchAssocLabel.size()) {
3430  << "fTrackModuleLabel.size() = "<<fTrackModuleLabel.size()<<" does not match "
3431  << "fCosmicTaggerAssocLabel.size() = "<<fFlashMatchAssocLabel.size();
3432  }
3433  if (GetNVertexAlgos() > kMaxVertexAlgos) {
3435  << "AnalysisTree currently supports only up to " << kMaxVertexAlgos
3436  << " tracking algorithms, but " << GetNVertexAlgos() << " are specified."
3437  << "\nYou can increase kMaxVertexAlgos and recompile.";
3438  } // if too many trackers
3439 
3440  // Build my Cryostat boundaries array...Taken from Tyler Alion in Geometry Core. Should still return the same values for uBoone.
3441  ActiveBounds[0] = ActiveBounds[2] = ActiveBounds[4] = DBL_MAX;
3442  ActiveBounds[1] = ActiveBounds[3] = ActiveBounds[5] = -DBL_MAX;
3443  // assume single cryostats
3444  auto const* geom = lar::providerFrom<geo::Geometry>();
3445  for (geo::TPCGeo const& TPC: geom->IterateTPCs()) {
3446  // get center in world coordinates
3447  double origin[3] = {0.};
3448  double center[3] = {0.};
3449  TPC.LocalToWorld(origin, center);
3450  double tpcDim[3] = {TPC.HalfWidth(), TPC.HalfHeight(), 0.5*TPC.Length() };
3451 
3452  if( center[0] - tpcDim[0] < ActiveBounds[0] ) ActiveBounds[0] = center[0] - tpcDim[0];
3453  if( center[0] + tpcDim[0] > ActiveBounds[1] ) ActiveBounds[1] = center[0] + tpcDim[0];
3454  if( center[1] - tpcDim[1] < ActiveBounds[2] ) ActiveBounds[2] = center[1] - tpcDim[1];
3455  if( center[1] + tpcDim[1] > ActiveBounds[3] ) ActiveBounds[3] = center[1] + tpcDim[1];
3456  if( center[2] - tpcDim[2] < ActiveBounds[4] ) ActiveBounds[4] = center[2] - tpcDim[2];
3457  if( center[2] + tpcDim[2] > ActiveBounds[5] ) ActiveBounds[5] = center[2] + tpcDim[2];
3458  } // for all TPC
3459  std::cout << "Active Boundaries: "
3460  << "\n\tx: " << ActiveBounds[0] << " to " << ActiveBounds[1]
3461  << "\n\ty: " << ActiveBounds[2] << " to " << ActiveBounds[3]
3462  << "\n\tz: " << ActiveBounds[4] << " to " << ActiveBounds[5]
3463  << std::endl;
3464 } // dune::AnalysisTree::AnalysisTree()
bool fSaveRawDigitInfo
whether to extract and save Hit information
bool fSavePFParticleInfo
whether to extract and save Shower information
std::vector< std::string > fCosmicTaggerAssocLabel
whether to extract and save CNN information
bool fSaveAuxDetInfo
whether to extract and save auxiliary detector data
std::string string
Definition: nybbler.cc:12
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Geometry information for a single TPC.
Definition: TPCGeo.h:38
bool fSaveCaloCosmics
save calorimetry information for cosmics
std::string fSpacePointSolverModuleLabel
std::vector< std::string > fMVAPIDShowerModuleLabel
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
std::vector< std::string > fContainmentTaggerAssocLabel
std::vector< std::string > fMCT0FinderLabel
std::vector< std::string > fVertexModuleLabel
constexpr int kMaxTrackers
constexpr int kMaxVertexAlgos
std::vector< std::string > fParticleIDModuleLabel
std::vector< std::string > fFlashT0FinderLabel
bool fSaveExternCounterInfo
whether to extract and save Flash information
std::string fCosmicClusterTaggerAssocLabel
std::vector< std::string > fMVAPIDTrackModuleLabel
bool fSaveVertexInfo
whether to extract and save Track information
std::string fPandoraNuVertexModuleLabel
std::vector< std::string > fCalorimetryModuleLabel
bool fUseBuffer
whether to use a permanent buffer (faster, huge memory)
size_t GetNVertexAlgos() const
std::string fPFParticleModuleLabel
bool fSaveProtoInfo
whether to extract and save Genie information
bool fSaveTrackInfo
whether to extract and save Raw Digit information
bool fSaveClusterInfo
whether to extract and save Vertex information
bool isCosmics
if it contains cosmics
std::vector< std::string > fTrackModuleLabel
bool fSavePandoraNuVertexInfo
whether to extract and save Cluster information
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
bool fSaveFlashInfo
whether to extract and save nu vertex information from Pandora
std::string fExternalCounterModuleLabel
bool fSaveMCShowerInfo
whether to extract and save Geant information
bool fSaveMCTrackInfo
whether to extract and save MC Shower information
def center(depos, point)
Definition: depos.py:117
bool fSaveGeantInfo
whether to extract and save ProtDUNE beam simulation information
float fG4minE
Energy threshold to save g4 particle info.
size_t GetNTrackers() const
Returns the number of trackers configured.
bool bIgnoreMissingShowers
whether to ignore missing shower information
bool fSaveCnnInfo
whether to extract and save SpacePointSolver information
std::vector< std::string > fFlashMatchAssocLabel
std::vector< std::string > fShowerModuleLabel
bool fSaveHitInfo
whether to extract and save MC Track information
bool fSaveGenieInfo
whether to extract and save CRY particle data
constexpr Point origin()
Returns a origin position with a point of the specified type.
Definition: geo_vectors.h:227
QTextStream & endl(QTextStream &s)
bool fSaveSpacePointSolverInfo
whether to extract and save PFParticle information
bool fSaveShowerInfo
whether to extract and save External Counter information
AnalysisTree::~AnalysisTree ( )
virtual

Definition at line 3467 of file AnalysisTree_module.cc.

3468 {
3469  DestroyData();
3470 }
void DestroyData()
Destroy the local buffers (existing branches will point to invalid address!)

Member Function Documentation

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

read access to event

transfer the run and subrun data to the tree data object

Definition at line 3538 of file AnalysisTree_module.cc.

3539 {
3540  std::cout << "Analysing.\n\n";
3541  //services
3544 
3545  // collect the sizes which might me needed to resize the tree data structure:
3546  bool isMC = !evt.isRealData();
3547 
3548  // * hits
3549  std::vector<art::Ptr<recob::Hit> > hitlist;
3550  auto hitListHandle = evt.getHandle< std::vector<recob::Hit> >(fHitsModuleLabel);
3551  if (hitListHandle)
3552  art::fill_ptr_vector(hitlist, hitListHandle);
3553 
3554  // * clusters
3555  art::Handle< std::vector<recob::Cluster> > clusterListHandle;
3556  std::vector<art::Ptr<recob::Cluster> > clusterlist;
3557  if (fSaveClusterInfo){
3558  clusterListHandle = evt.getHandle< std::vector<recob::Cluster> >(fClusterModuleLabel);
3559  if (clusterListHandle)
3560  art::fill_ptr_vector(clusterlist, clusterListHandle);
3561  }
3562 
3563  // * spacepoints
3564  art::Handle<std::vector<recob::SpacePoint>> spacepointListHandle;
3565  art::Handle<std::vector<recob::PointCharge>> pointchargeListHandle;
3567  std::cout << "Saving SpacepointSolver info.";
3568  spacepointListHandle = evt.getHandle<std::vector<recob::SpacePoint>>(fSpacePointSolverModuleLabel);
3569  pointchargeListHandle = evt.getHandle<std::vector<recob::PointCharge>>(fSpacePointSolverModuleLabel);
3570  if (spacepointListHandle->size() != pointchargeListHandle->size()) {
3571  throw cet::exception("tutorial::ReadSpacePointAndCnn")
3572  << "size of point and charge containers must be equal" << std::endl;
3573  }
3574  }
3575 
3576  // * flashes
3577  std::vector<art::Ptr<recob::OpFlash> > flashlist;
3578  auto flashListHandle = evt.getHandle< std::vector<recob::OpFlash> >(fOpFlashModuleLabel);
3579  if (flashListHandle)
3580  art::fill_ptr_vector(flashlist, flashListHandle);
3581 
3582  // * External Counters
3583  std::vector<art::Ptr<raw::ExternalTrigger> > countlist;
3584  auto countListHandle = evt.getHandle< std::vector<raw::ExternalTrigger> >(fExternalCounterModuleLabel);
3585  if (countListHandle)
3586  art::fill_ptr_vector(countlist, countListHandle);
3587 
3588  // * MC truth information
3589  art::Handle< std::vector<simb::MCTruth> > mctruthListHandle;
3590  std::vector<art::Ptr<simb::MCTruth> > mclist;
3591  if (isMC){
3592  mctruthListHandle = evt.getHandle< std::vector<simb::MCTruth> >(fGenieGenModuleLabel);
3593  if (mctruthListHandle)
3594  art::fill_ptr_vector(mclist, mctruthListHandle);
3595  }
3596 
3597  // *MC truth cosmic generator information
3598  art::Handle< std::vector<simb::MCTruth> > mctruthcryListHandle;
3599  std::vector<art::Ptr<simb::MCTruth> > mclistcry;
3600  if (isMC && fSaveCryInfo){
3601  mctruthcryListHandle = evt.getHandle< std::vector<simb::MCTruth> >(fCryGenModuleLabel);
3602  if (mctruthcryListHandle) {
3603  art::fill_ptr_vector(mclistcry, mctruthcryListHandle);
3604  }
3605  else{
3606  // If we requested this info but it doesn't exist, don't use it.
3607  fSaveCryInfo = false;
3608  mf::LogError("AnalysisTree") << "Requested CRY information but none exists, hence not saving CRY information.";
3609  }
3610  }
3611 
3612  // ProtoDUNE beam generator information
3613  art::Handle< std::vector<simb::MCTruth> > mctruthprotoListHandle;
3614  std::vector<art::Ptr<simb::MCTruth> > mclistproto;
3615  if (isMC && fSaveProtoInfo){
3616  mctruthprotoListHandle = evt.getHandle< std::vector<simb::MCTruth> >(fProtoGenModuleLabel);
3617  if (mctruthprotoListHandle) {
3618  art::fill_ptr_vector(mclistproto,mctruthprotoListHandle);
3619  }
3620  else{
3621  // If we requested this info but it doesn't exist, don't use it.
3622  fSaveProtoInfo = false;
3623  mf::LogError("AnalysisTree") << "Requested protoDUNE beam generator information but none exists, hence not saving protoDUNE beam generator information.";
3624  }
3625  }
3626 
3627  // *MC Shower information
3629  if (isMC)
3630  mcshowerh = evt.getHandle< std::vector<sim::MCShower> >(fMCShowerModuleLabel);
3631 
3632  int nMCShowers = 0;
3633  if (fSaveMCShowerInfo && mcshowerh.isValid())
3634  nMCShowers = mcshowerh->size();
3635 
3636  // *MC Track information
3638  if (isMC)
3639  mctrackh = evt.getHandle< std::vector<sim::MCTrack> >(fMCTrackModuleLabel);
3640 
3641  int nMCTracks = 0;
3642  if (fSaveMCTrackInfo && mctrackh.isValid())
3643  nMCTracks = mctrackh->size();
3644 
3645  art::Ptr<simb::MCTruth> mctruthcry;
3646  int nCryPrimaries = 0;
3647 
3648  if (fSaveCryInfo){
3649  mctruthcry = mclistcry[0];
3650  nCryPrimaries = mctruthcry->NParticles();
3651  }
3652 
3653  art::Ptr<simb::MCTruth> mctruthproto;
3654  int nProtoPrimaries = 0;
3655  if( fSaveProtoInfo){
3656  mctruthproto = mclistproto[0];
3657  nProtoPrimaries = mctruthproto->NParticles();
3658  }
3659 
3660  int nGeniePrimaries = 0, nGEANTparticles = 0;
3661 
3662  art::Ptr<simb::MCTruth> mctruth;
3663 
3664 
3665  if (isMC) { //is MC
3666 
3667  //find origin
3668  auto const allmclists = evt.getMany<std::vector<simb::MCTruth>>();
3669  for(size_t mcl = 0; mcl < allmclists.size(); ++mcl){
3670  art::Handle< std::vector<simb::MCTruth> > mclistHandle = allmclists[mcl];
3671  for(size_t m = 0; m < mclistHandle->size(); ++m){
3672  art::Ptr<simb::MCTruth> mct(mclistHandle, m);
3673  if (mct->Origin() == simb::kCosmicRay) isCosmics = true;
3674  }
3675  }
3676  if (fSaveCaloCosmics) isCosmics = false; //override to save calo info
3677 
3678  // GENIE
3679  if (!mclist.empty()){//at least one mc record
3680 
3681  // double maxenergy = -1;
3682  // int imc0 = 0;
3683  // for (std::map<art::Ptr<simb::MCTruth>,double>::iterator ii=mctruthemap.begin(); ii!=mctruthemap.end(); ++ii){
3684  // if ((ii->second)>maxenergy){
3685  // maxenergy = ii->second;
3686  // mctruth = ii->first;
3687  // imc = imc0;
3688  // }
3689  // imc0++;
3690  // }
3691 
3692  //imc = 0; //set imc to 0 to solve a confusion for BNB+cosmic files where there are two MCTruth
3693  mctruth = mclist[0];
3694 
3695  if (mctruth->NeutrinoSet()) nGeniePrimaries = mctruth->NParticles();
3696  //} //end (fSaveGenieInfo)
3697 
3698  const sim::ParticleList& plist = pi_serv->ParticleList();
3699  nGEANTparticles = plist.size();
3700 
3701  // to know the number of particles in AV would require
3702  // looking at all of them; so we waste some memory here
3703  } // if have MC truth
3704  MF_LOG_DEBUG("AnalysisTree") << "Expected "
3705  << nGEANTparticles << " GEANT particles, "
3706  << nGeniePrimaries << " GENIE particles";
3707  } // if MC
3708 
3709  // SpacePointSolver information
3710  int nSpacePoints = 0;
3711  if (fSaveSpacePointSolverInfo && spacepointListHandle.isValid())
3712  nSpacePoints = spacepointListHandle->size();
3713 
3714  CreateData(); // tracker data is created with default constructor
3715  if (fSaveGenieInfo)
3716  fData->ResizeGenie(nGeniePrimaries);
3717  if (fSaveCryInfo)
3718  fData->ResizeCry(nCryPrimaries);
3719  if (fSaveProtoInfo)
3720  fData->ResizeProto(nProtoPrimaries);
3721  if (fSaveGeantInfo)
3722  fData->ResizeGEANT(nGEANTparticles);
3723  if (fSaveMCShowerInfo)
3724  fData->ResizeMCShower(nMCShowers);
3725  if (fSaveMCTrackInfo)
3726  fData->ResizeMCTrack(nMCTracks);
3728  fData->ResizeSpacePointSolver(nSpacePoints);
3729 
3730  fData->ClearLocalData(); // don't bother clearing tracker data yet
3731 
3732  const size_t NTrackers = GetNTrackers(); // number of trackers passed into fTrackModuleLabel
3733  const size_t NShowerAlgos = GetNShowerAlgos(); // number of shower algorithms into fShowerModuleLabel
3734  const size_t NHits = hitlist.size(); // number of hits
3735  const size_t NVertexAlgos = GetNVertexAlgos(); // number of vertex algos
3736  const size_t NClusters = clusterlist.size(); //number of clusters
3737  const size_t NFlashes = flashlist.size(); // number of flashes
3738  const size_t NExternCounts = countlist.size(); // number of External Counters
3739  // make sure there is the data, the tree and everything;
3740  CreateTree();
3741 
3742  /// transfer the run and subrun data to the tree data object
3743  // fData->RunData = RunData;
3744  fData->SubRunData = SubRunData;
3745 
3746  fData->isdata = int(!isMC);
3747 
3748  // * raw trigger
3749  std::vector<art::Ptr<raw::Trigger>> triggerlist;
3750  auto triggerListHandle = evt.getHandle< std::vector<raw::Trigger>>(fDigitModuleLabel);
3751  if (triggerListHandle)
3752  art::fill_ptr_vector(triggerlist, triggerListHandle);
3753 
3754  if (triggerlist.size()){
3755  fData->triggernumber = triggerlist[0]->TriggerNumber();
3756  fData->triggertime = triggerlist[0]->TriggerTime();
3757  fData->beamgatetime = triggerlist[0]->BeamGateTime();
3758  fData->triggerbits = triggerlist[0]->TriggerBits();
3759  }
3760 
3761  // * vertices
3762  std::vector< art::Handle< std::vector<recob::Vertex> > > vertexListHandle(NVertexAlgos);
3763  std::vector< std::vector<art::Ptr<recob::Vertex> > > vertexlist(NVertexAlgos);
3764  for (unsigned int it = 0; it < NVertexAlgos; ++it){
3765  vertexListHandle[it] = evt.getHandle< std::vector<recob::Vertex> >(fVertexModuleLabel[it]);
3766  if (vertexListHandle[it])
3767  art::fill_ptr_vector(vertexlist[it], vertexListHandle[it]);
3768  }
3769 
3770  // * PFParticles
3771  lar_pandora::PFParticleVector pfparticlelist;
3772  lar_pandora::PFParticlesToClusters pfParticleToClusterMap;
3773  lar_pandora::LArPandoraHelper::CollectPFParticles(evt, fPFParticleModuleLabel, pfparticlelist, pfParticleToClusterMap);
3774 
3775  // * tracks
3776  std::vector< art::Handle< std::vector<recob::Track> > > trackListHandle(NTrackers);
3777  std::vector< std::vector<art::Ptr<recob::Track> > > tracklist(NTrackers);
3778  for (unsigned int it = 0; it < NTrackers; ++it){
3779  trackListHandle[it] = evt.getHandle< std::vector<recob::Track> >(fTrackModuleLabel[it]);
3780  if (trackListHandle[it])
3781  art::fill_ptr_vector(tracklist[it], trackListHandle[it]);
3782  }
3783 
3784  // * showers
3785  // It seems that sometimes the shower data product does not exist;
3786  // in that case, we store a nullptr in place of the pointer to the vector;
3787  // the data structure itself is owned by art::Event and we should not try
3788  // to manage its memory
3789  std::vector<std::vector<recob::Shower> const*> showerList;
3790  std::vector< art::Handle< std::vector<recob::Shower> > > showerListHandle;
3791  showerList.reserve(fShowerModuleLabel.size());
3792  if (fSaveShowerInfo) {
3793  for (art::InputTag ShowerInputTag: fShowerModuleLabel) {
3794  auto ShowerHandle = evt.getHandle<std::vector<recob::Shower>>(ShowerInputTag);
3795  if (!ShowerHandle) {
3796  showerList.push_back(nullptr);
3797  if (!bIgnoreMissingShowers) {
3799  << "Showers with input tag '" << ShowerInputTag.encode()
3800  << "' were not found in the event."
3801  " If you really know what you are doing,"
3802  " set AnalysisTree's configuration parameter IgnoreMissingShowers"
3803  " to \"true\"; the lack of any shower set will be tolerated,"
3804  " and the shower list in the corresponding event will be set to"
3805  " a list of one shower, with an invalid ID.\n";
3806  } // if bIgnoreMissingShowers
3807  else {
3808  // this message is more alarming than strictly necessary; by design.
3809  mf::LogError("AnalysisTree")
3810  << "No showers found for input tag '" << ShowerInputTag.encode()
3811  << "' ; FILLING WITH FAKE DATA AS FOR USER'S REQUEST";
3812  }
3813  }
3814 
3815  else showerList.push_back(ShowerHandle.product());
3816 
3817  showerListHandle.push_back(ShowerHandle); // either way, put it into the handle list
3818 
3819  }
3820 
3821  } // for shower input tag
3822 
3823 
3824 
3825  std::vector<const sim::AuxDetSimChannel*> fAuxDetSimChannels;
3826  if (fSaveAuxDetInfo){
3827  evt.getView(fLArG4ModuleLabel, fAuxDetSimChannels);
3828  }
3829 
3830  std::vector<const sim::SimChannel*> fSimChannels;
3831  if (isMC && fSaveGeantInfo)
3832  evt.getView(fLArG4ModuleLabel, fSimChannels);
3833 
3834  fData->run = evt.run();
3835  fData->subrun = evt.subRun();
3836  fData->event = evt.id().event();
3837 
3838  art::Timestamp ts = evt.time();
3839  TTimeStamp tts(ts.timeHigh(), ts.timeLow());
3840  fData->evttime = tts.AsDouble();
3841 
3842  //copied from MergeDataPaddles.cxx
3843  auto beam = evt.getHandle< raw::BeamInfo >("beamdata");
3844  if (beam){
3845  fData->beamtime = (double)beam->get_t_ms();
3846  fData->beamtime/=1000.; //in second
3847  std::map<std::string, std::vector<double>> datamap = beam->GetDataMap();
3848  if (datamap["E:TOR860"].size()){
3849  fData->potbnb = datamap["E:TOR860"][0];
3850  }
3851  if (datamap["E:TORTGT"].size()){
3852  fData->potnumitgt = datamap["E:TORTGT"][0];
3853  }
3854  if (datamap["E:TOR101"].size()){
3855  fData->potnumi101 = datamap["E:TOR101"][0];
3856  }
3857  }
3858 
3859 
3860  // std::cout<<detProp.NumberTimeSamples()<<" "<<detProp.ReadOutWindowSize()<<std::endl;
3861  // std::cout<<geom->DetHalfHeight()*2<<" "<<geom->DetHalfWidth()*2<<" "<<geom->DetLength()<<std::endl;
3862  // std::cout<<geom->Nwires(0)<<" "<<geom->Nwires(1)<<" "<<geom->Nwires(2)<<std::endl;
3863  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
3864  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clockData);
3865 
3866  //hit information
3867  if (fSaveHitInfo){
3868  fData->no_hits = (int) NHits;
3869  fData->no_hits_stored = TMath::Min( (int) NHits, (int) kMaxHits);
3870  if (NHits > kMaxHits) {
3871  // got this error? consider increasing kMaxHits
3872  // (or ask for a redesign using vectors)
3873  mf::LogError("AnalysisTree:limits") << "event has " << NHits
3874  << " hits, only kMaxHits=" << kMaxHits << " stored in tree";
3875  }
3876  for (size_t i = 0; i < NHits && i < kMaxHits ; ++i){//loop over hits
3877  fData->hit_channel[i] = hitlist[i]->Channel();
3878  fData->hit_tpc[i] = hitlist[i]->WireID().TPC;
3879  fData->hit_plane[i] = hitlist[i]->WireID().Plane;
3880  fData->hit_wire[i] = hitlist[i]->WireID().Wire;
3881  fData->hit_peakT[i] = hitlist[i]->PeakTime();
3882  fData->hit_charge[i] = hitlist[i]->Integral();
3883  fData->hit_ph[i] = hitlist[i]->PeakAmplitude();
3884  fData->hit_startT[i] = hitlist[i]->PeakTimeMinusRMS();
3885  fData->hit_endT[i] = hitlist[i]->PeakTimePlusRMS();
3886  fData->hit_rms[i] = hitlist[i]->RMS();
3887  fData->hit_goodnessOfFit[i] = hitlist[i]->GoodnessOfFit();
3888  fData->hit_multiplicity[i] = hitlist[i]->Multiplicity();
3889  //std::vector<double> xyz = bt_serv->HitToXYZ(hitlist[i]);
3890  //when the size of simIDEs is zero, the above function throws an exception
3891  //and crashes, so check that the simIDEs have non-zero size before
3892  //extracting hit true XYZ from simIDEs
3893  if (isMC){
3894  std::vector<const sim::IDE*> ides;
3895  try{
3896  ides= bt_serv->HitToSimIDEs_Ps(clockData, hitlist[i]);
3897  }
3898  catch(...){}
3899  if (ides.size()>0){
3900  std::vector<double> xyz = bt_serv->SimIDEsToXYZ(ides);
3901  fData->hit_trueX[i] = xyz[0];
3902  }
3903  }
3904 
3905  /*
3906  for (unsigned int it=0; it<fTrackModuleLabel.size();++it){
3907  art::FindManyP<recob::Track> fmtk(hitListHandle,evt,fTrackModuleLabel[it]);
3908  if (fmtk.at(i).size()!=0){
3909  hit_trkid[it][i] = fmtk.at(i)[0]->ID();
3910  }
3911  else
3912  hit_trkid[it][i] = 0;
3913  }
3914  */
3915 
3916  if (fSaveRawDigitInfo){
3917  //Hit to RawDigit information
3918  art::FindManyP<raw::RawDigit> fmrd(hitListHandle,evt,fHitsModuleLabel);
3919  if (hitlist[i]->WireID().Plane==2)
3920  {
3921  int dataSize = fmrd.at(i)[0]->Samples();
3922  short ped = fmrd.at(i)[0]->GetPedestal();
3923 
3924  std::vector<short> rawadc(dataSize);
3925  raw::Uncompress(fmrd.at(i)[0]->ADCs(), rawadc, fmrd.at(i)[0]->Compression());
3926  int t0 = hitlist[i]->PeakTime() - 3*(hitlist[i]->RMS());
3927  if (t0<0) t0 = 0;
3928  int t1 = hitlist[i]->PeakTime() + 3*(hitlist[i]->RMS());
3929  if (t1>=dataSize) t1 = dataSize-1;
3930  fData->rawD_ph[i] = -1;
3931  fData->rawD_peakT[i] = -1;
3932  for (int j = t0; j<=t1; ++j){
3933  if (rawadc[j]-ped>fData->rawD_ph[i]){
3934  fData->rawD_ph[i] = rawadc[j]-ped;
3935  fData->rawD_peakT[i] = j;
3936  }
3937  }
3938  fData->rawD_charge[i] = 0;
3939  fData->rawD_fwhh[i] = 0;
3940  double mean_t = 0.0;
3941  double mean_t2 = 0.0;
3942  for (int j = t0; j<=t1; ++j){
3943  if (rawadc[j]-ped>=0.5*fData->rawD_ph[i]){
3944  ++fData->rawD_fwhh[i];
3945  }
3946  if (rawadc[j]-ped>=0.1*fData->rawD_ph[i]){
3947  fData->rawD_charge[i] += rawadc[j]-ped;
3948  mean_t += (double)j*(rawadc[j]-ped);
3949  mean_t2 += (double)j*(double)j*(rawadc[j]-ped);
3950  }
3951  }
3952  mean_t/=fData->rawD_charge[i];
3953  mean_t2/=fData->rawD_charge[i];
3954  fData->rawD_rms[i] = sqrt(mean_t2-mean_t*mean_t);
3955  }
3956  } // Save RawDigitInfo
3957 
3958  if (!evt.isRealData()&&!isCosmics){
3959  fData -> hit_nelec[i] = 0;
3960  fData -> hit_energy[i] = 0;
3961  const sim::SimChannel* chan = 0;
3962  for(size_t sc = 0; sc < fSimChannels.size(); ++sc){
3963  if(fSimChannels[sc]->Channel() == hitlist[i]->Channel()) chan = fSimChannels[sc];
3964  }
3965  if (chan){
3966  for(auto const& mapitr : chan->TDCIDEMap()){
3967  // loop over the vector of IDE objects.
3968  for(auto const& ide : mapitr.second){
3969  fData -> hit_nelec[i] += ide.numElectrons;
3970  fData -> hit_energy[i] += ide.energy;
3971  }
3972  }
3973  }
3974  }
3975  } // Loop over hits
3976 
3977  hitListHandle = evt.getHandle< std::vector<recob::Hit> >(fHitsModuleLabel);
3978  if (hitListHandle) {
3979  //Find tracks associated with hits
3980  art::FindManyP<recob::Track> fmtk(hitListHandle,evt,fTrackModuleLabel[0]);
3981  for (size_t i = 0; i < NHits && i < kMaxHits ; ++i){//loop over hits
3982  if (fmtk.isValid()){
3983  if (fmtk.at(i).size()!=0){
3984  fData->hit_trkid[i] = fmtk.at(i)[0]->ID();
3985  fData->hit_trkKey[i] = fmtk.at(i)[0].key();
3986 
3987  }
3988  else
3989  fData->hit_trkid[i] = -1;
3990  }
3991  }
3992  }
3993 
3994  //In the case of ClusterCrawler or linecluster, use "linecluster or clustercrawler" as HitModuleLabel.
3995  //using cchit will not make this association. In the case of gaushit, just use gaushit
3996  //Not initializing clusterID to -1 since some clustering algorithms assign negative IDs!
3997 
3998  hitListHandle = evt.getHandle< std::vector<recob::Hit> >(fHitsModuleLabel);
3999  if (hitListHandle) {
4000  //Find clusters and spacepoints associated with hits
4001  art::FindManyP<recob::Cluster> fmcl(hitListHandle,evt,fClusterModuleLabel);
4002  art::FindManyP<recob::SpacePoint> fmsp(hitListHandle,evt,fSpacePointSolverModuleLabel);
4003  for (size_t i = 0; i < NHits && i < kMaxHits ; ++i){//loop over hits
4004  if (fmcl.isValid() && fmcl.at(i).size()!=0){
4005  fData->hit_clusterid[i] = fmcl.at(i)[0]->ID();
4006  fData->hit_clusterKey[i] = fmcl.at(i)[0].key();
4007  // std::cout << "ClusterID " <<
4008  }
4009  if(fmsp.isValid() && fmsp.at(i).size()!=0){
4010  fData->hit_spacepointid[i] = fmsp.at(i)[0]->ID();
4011  fData->hit_spacepointKey[i] = fmsp.at(i)[0].key();
4012  }
4013  }
4014  }
4015  }// end (fSaveHitInfo)
4016 
4017 
4019  lar_pandora::PFParticleVector particleVector;
4021  lar_pandora::VertexVector vertexVector;
4022  lar_pandora::PFParticlesToVertices particlesToVertices;
4023  lar_pandora::LArPandoraHelper::CollectVertices(evt, fPandoraNuVertexModuleLabel, vertexVector, particlesToVertices);
4024 
4025  short nprim = 0;
4026  for (unsigned int n = 0; n < particleVector.size(); ++n) {
4027  const art::Ptr<recob::PFParticle> particle = particleVector.at(n);
4028  if(particle->IsPrimary()) nprim++;
4029  }
4030 
4031  if (nprim > kMaxVertices){
4032  // got this error? consider increasing kMaxClusters
4033  // (or ask for a redesign using vectors)
4034  mf::LogError("AnalysisTree:limits") << "event has " << nprim
4035  << " nu neutrino vertices, only kMaxVertices=" << kMaxVertices << " stored in tree";
4036  }
4037 
4038  fData->nnuvtx = nprim;
4039 
4040  short iv = 0;
4041  for (unsigned int n = 0; n < particleVector.size(); ++n) {
4042  const art::Ptr<recob::PFParticle> particle = particleVector.at(n);
4043  if(particle->IsPrimary()) {
4044  lar_pandora::PFParticlesToVertices::const_iterator vIter = particlesToVertices.find(particle);
4045  if (particlesToVertices.end() != vIter) {
4046  const lar_pandora::VertexVector &vertexVector = vIter->second;
4047  if (!vertexVector.empty()) {
4048  if (vertexVector.size() == 1) {
4049  const art::Ptr<recob::Vertex> vertex = *(vertexVector.begin());
4050  double xyz[3] = {0.0, 0.0, 0.0} ;
4051  vertex->XYZ(xyz);
4052  fData->nuvtxx[iv] = xyz[0];
4053  fData->nuvtxy[iv] = xyz[1];
4054  fData->nuvtxz[iv] = xyz[2];
4055  fData->nuvtxpdg[iv] = particle->PdgCode();
4056  iv++;
4057  }
4058  }
4059  }
4060  }
4061  }
4062  } // save PandoraNuVertexInfo
4063 
4064 
4065  if (fSaveClusterInfo){
4066  fData->nclusters = (int) NClusters;
4067  if (NClusters > kMaxClusters){
4068  // got this error? consider increasing kMaxClusters
4069  // (or ask for a redesign using vectors)
4070  mf::LogError("AnalysisTree:limits") << "event has " << NClusters
4071  << " clusters, only kMaxClusters=" << kMaxClusters << " stored in tree";
4072  }
4073  for(unsigned int ic=0; ic<NClusters;++ic){//loop over clusters
4074  art::Ptr<recob::Cluster> clusterholder(clusterListHandle, ic);
4075  const recob::Cluster& cluster = *clusterholder;
4076  fData->clusterId[ic] = cluster.ID();
4077  fData->clusterView[ic] = cluster.View();
4078  fData->cluster_isValid[ic] = cluster.isValid();
4079  fData->cluster_StartCharge[ic] = cluster.StartCharge();
4080  fData->cluster_StartAngle[ic] = cluster.StartAngle();
4081  fData->cluster_EndCharge[ic] = cluster.EndCharge();
4082  fData->cluster_EndAngle[ic] = cluster.EndAngle();
4083  fData->cluster_Integral[ic] = cluster.Integral();
4084  fData->cluster_IntegralAverage[ic] = cluster.IntegralAverage();
4085  fData->cluster_SummedADC[ic] = cluster.SummedADC();
4086  fData->cluster_SummedADCaverage[ic] = cluster.SummedADCaverage();
4087  fData->cluster_MultipleHitDensity[ic] = cluster.MultipleHitDensity();
4088  fData->cluster_Width[ic] = cluster.Width();
4089  fData->cluster_NHits[ic] = cluster.NHits();
4090  fData->cluster_StartWire[ic] = cluster.StartWire();
4091  fData->cluster_StartTick[ic] = cluster.StartTick();
4092  fData->cluster_EndWire[ic] = cluster.EndWire();
4093  fData->cluster_EndTick[ic] = cluster.EndTick();
4094 
4095  //Cosmic Tagger information for cluster
4096  art::FindManyP<anab::CosmicTag> fmcct(clusterListHandle,evt,fCosmicClusterTaggerAssocLabel);
4097  if (fmcct.isValid()){
4098  fData->cluncosmictags_tagger[ic] = fmcct.at(ic).size();
4099  if (fmcct.at(ic).size()>0){
4100  if(fmcct.at(ic).size()>1)
4101  std::cerr << "\n Warning : more than one cosmic tag per cluster in module! assigning the first tag to the cluster" << fCosmicClusterTaggerAssocLabel;
4102  fData->clucosmicscore_tagger[ic] = fmcct.at(ic).at(0)->CosmicScore();
4103  fData->clucosmictype_tagger[ic] = fmcct.at(ic).at(0)->CosmicType();
4104  }
4105  }
4106  }//end loop over clusters
4107  }//end fSaveClusterInfo
4108 
4110  fData->nspacepoints = (unsigned int) nSpacePoints;
4111 
4112  // Largely copied from Robert Sulej's ReadSpacePointAndCnn_module.cc
4113  if(fSaveCnnInfo) {
4115  if(cluResults) {
4116  size_t emLikeIdx = cluResults->getIndex("em"); // at which index EM-like is stored in CNN output vector
4117 
4118  const art::FindManyP<recob::Hit> hitsFromClusters(cluResults->dataHandle(), evt, cluResults->dataTag());
4119  const art::FindManyP<recob::SpacePoint> spFromHits(hitListHandle, evt, fSpacePointSolverModuleLabel);
4120 
4121  std::vector<size_t> sizeScore(nSpacePoints, 0); // keep track of the max size of a cluster containing hit associated to spacepoint
4122 
4123  for(size_t c = 0; c < cluResults->size(); ++c) {
4124  const std::vector< art::Ptr<recob::Hit> > & hits = hitsFromClusters.at(c);
4125  std::array<float, MVA_LENGTH> cnn_out = cluResults->getOutput(c);
4126 
4127  for (auto& hptr : hits) {
4128  const std::vector< art::Ptr<recob::SpacePoint> > & sp = spFromHits.at(hptr.key());
4129  for(const auto & spptr : sp) { // Should always be just one associated spacepoint
4130  if(hits.size() > sizeScore[spptr.key()]) {
4131  sizeScore[spptr.key()] = hits.size();
4132  fData->SpacePointEmScore[spptr.key()] = cnn_out[emLikeIdx];
4133  }
4134  } // Loop over associated spacepoints
4135  } // Loop over hits
4136  } // Loop over cluResults
4137  } // If cluResults
4138  }
4139 
4140  for (unsigned int is = 0; is < (unsigned int)nSpacePoints;
4141  ++is) { // loop over spacepoints
4142  fData->SpacePointX[is] = (*spacepointListHandle)[is].XYZ()[0];
4143  fData->SpacePointY[is] = (*spacepointListHandle)[is].XYZ()[1];
4144  fData->SpacePointZ[is] = (*spacepointListHandle)[is].XYZ()[2];
4145 
4146  fData->SpacePointQ[is] = (*pointchargeListHandle)[is].charge();
4147 
4148  fData->SpacePointErrX[is] = (*spacepointListHandle)[is].ErrXYZ()[0];
4149  fData->SpacePointErrY[is] = (*spacepointListHandle)[is].ErrXYZ()[1];
4150  fData->SpacePointErrZ[is] = (*spacepointListHandle)[is].ErrXYZ()[2];
4151 
4152  fData->SpacePointID[is] = (*spacepointListHandle)[is].ID();
4153 
4154  fData->SpacePointID[is] = (*spacepointListHandle)[is].Chisq();
4155  }//end loop over spacepoints
4156  }//end fSpacePointSolverInfo
4157 
4158  if (fSaveFlashInfo){
4159  fData->no_flashes = (int) NFlashes;
4160  if (NFlashes > kMaxFlashes) {
4161  // got this error? consider increasing kMaxHits
4162  // (or ask for a redesign using vectors)
4163  mf::LogError("AnalysisTree:limits") << "event has " << NFlashes
4164  << " flashes, only kMaxFlashes=" << kMaxFlashes << " stored in tree";
4165  }
4166 
4167  std::sort(flashlist.begin(), flashlist.end(), recob::OpFlashPtrSortByPE);
4168 
4169  for (size_t i = 0; i < NFlashes && i < kMaxFlashes ; ++i){//loop over hits
4170  fData->flash_time[i] = flashlist[i]->Time();
4171  fData->flash_pe[i] = flashlist[i]->TotalPE();
4172  fData->flash_ycenter[i] = flashlist[i]->YCenter();
4173  fData->flash_zcenter[i] = flashlist[i]->ZCenter();
4174  fData->flash_ywidth[i] = flashlist[i]->YWidth();
4175  fData->flash_zwidth[i] = flashlist[i]->ZWidth();
4176  fData->flash_timewidth[i] = flashlist[i]->TimeWidth();
4177  }
4178  }
4179 
4181  fData->no_ExternCounts = (int) NExternCounts;
4182  if (NExternCounts > kMaxExternCounts) {
4183  // got this error? consider increasing kMaxHits
4184  // (or ask for a redesign using vectors)
4185  mf::LogError("AnalysisTree:limits") << "event has " << NExternCounts
4186  << " External Counters, only kMaxExternCounts=" << kMaxExternCounts << " stored in tree";
4187  }
4188  for (size_t i = 0; i < NExternCounts && i < kMaxExternCounts ; ++i){//loop over hits
4189  fData->externcounts_time[i] = countlist[i]->GetTrigTime();
4190  fData->externcounts_id[i] = countlist[i]->GetTrigID();
4191  }
4192  }
4193 
4194 
4195  // Declare object-ID-to-PFParticleID maps so we can assign hasPFParticle and PFParticleID to the tracks, showers, vertices.
4196  std::map<Short_t, Short_t> trackIDtoPFParticleIDMap, vertexIDtoPFParticleIDMap, showerIDtoPFParticleIDMap;
4197 
4198  //Save PFParticle information
4199  if (fSavePFParticleInfo){
4200  AnalysisTreeDataStruct::PFParticleDataStruct& PFParticleData = fData->GetPFParticleData();
4201  size_t NPFParticles = pfparticlelist.size();
4202 
4203  PFParticleData.SetMaxPFParticles(std::max(NPFParticles, (size_t) 1));
4204  PFParticleData.Clear(); // clear all the data
4205 
4206  PFParticleData.nPFParticles = (short) NPFParticles;
4207 
4208  // now set the tree addresses to the newly allocated memory;
4209  // this creates the tree branches in case they are not there yet
4211 
4212  if (NPFParticles > PFParticleData.GetMaxPFParticles()) {
4213  mf::LogError("AnalysisTree:limits") << "event has " << NPFParticles
4214  << " PFParticles, only "
4215  << PFParticleData.GetMaxPFParticles() << " stored in tree";
4216  }
4217 
4218  lar_pandora::PFParticleVector neutrinoPFParticles;
4219  lar_pandora::LArPandoraHelper::SelectNeutrinoPFParticles(pfparticlelist, neutrinoPFParticles);
4220  PFParticleData.pfp_numNeutrinos = neutrinoPFParticles.size();
4221 
4222  for (size_t i = 0; i < std::min(neutrinoPFParticles.size(), (size_t)kMaxNPFPNeutrinos); ++i) {
4223  PFParticleData.pfp_neutrinoIDs[i] = neutrinoPFParticles[i]->Self();
4224  }
4225 
4226  if (neutrinoPFParticles.size() > kMaxNPFPNeutrinos)
4227  std::cerr << "Warning: there were " << neutrinoPFParticles.size() << " reconstructed PFParticle neutrinos; only the first " << kMaxNPFPNeutrinos << " being stored in tree" << std::endl;
4228 
4229  // Get a PFParticle-to-vertex map.
4230  lar_pandora::VertexVector allPfParticleVertices;
4231  lar_pandora::PFParticlesToVertices pfParticleToVertexMap;
4232  lar_pandora::LArPandoraHelper::CollectVertices(evt, fPFParticleModuleLabel, allPfParticleVertices, pfParticleToVertexMap);
4233 
4234  // Get a PFParticle-to-track map.
4235  lar_pandora::TrackVector allPfParticleTracks;
4236  lar_pandora::PFParticlesToTracks pfParticleToTrackMap;
4237  lar_pandora::LArPandoraHelper::CollectTracks(evt, fPFParticleModuleLabel, allPfParticleTracks, pfParticleToTrackMap);
4238 
4239  // Get a PFParticle-to-shower map.
4240  lar_pandora::ShowerVector allPfParticleShowers;
4241  lar_pandora::PFParticlesToShowers pfParticleToShowerMap;
4242  lar_pandora::LArPandoraHelper::CollectShowers(evt, fPFParticleModuleLabel, allPfParticleShowers, pfParticleToShowerMap);
4243 
4244  for (size_t i = 0; i < NPFParticles && i < PFParticleData.GetMaxPFParticles() ; ++i){
4245  PFParticleData.pfp_selfID[i] = pfparticlelist[i]->Self();
4246  PFParticleData.pfp_isPrimary[i] = (Short_t)pfparticlelist[i]->IsPrimary();
4247  PFParticleData.pfp_numDaughters[i] = pfparticlelist[i]->NumDaughters();
4248  PFParticleData.pfp_parentID[i] = pfparticlelist[i]->Parent();
4249  PFParticleData.pfp_pdgCode[i] = pfparticlelist[i]->PdgCode();
4250  PFParticleData.pfp_isNeutrino[i] = lar_pandora::LArPandoraHelper::IsNeutrino(pfparticlelist[i]);
4251 
4252  // Set the daughter IDs.
4253  std::vector<size_t> daughterIDs = pfparticlelist[i]->Daughters();
4254 
4255  for (size_t j = 0; j < daughterIDs.size(); ++j)
4256  PFParticleData.pfp_daughterIDs[i][j] = daughterIDs[j];
4257 
4258  // Set the vertex ID.
4259  auto vertexMapIter = pfParticleToVertexMap.find(pfparticlelist[i]);
4260  if (vertexMapIter != pfParticleToVertexMap.end()) {
4261  lar_pandora::VertexVector pfParticleVertices = vertexMapIter->second;
4262 
4263  if (pfParticleVertices.size() > 1)
4264  std::cerr << "Warning: there was more than one vertex found for PFParticle with ID " << pfparticlelist[i]->Self() << ", storing only one" << std::endl;
4265 
4266  if (pfParticleVertices.size() > 0) {
4267  PFParticleData.pfp_vertexID[i] = pfParticleVertices.at(0)->ID();
4268  vertexIDtoPFParticleIDMap.insert(std::make_pair(pfParticleVertices.at(0)->ID(), pfparticlelist[i]->Self()));
4269  }
4270  }
4271  else
4272  std::cerr << "Warning: there was no vertex found for PFParticle with ID " << pfparticlelist[i]->Self() << std::endl;
4273 
4274  if (lar_pandora::LArPandoraHelper::IsTrack(pfparticlelist[i])){
4275  PFParticleData.pfp_isTrack[i] = 1;
4276 
4277  // Set the track ID.
4278  auto trackMapIter = pfParticleToTrackMap.find(pfparticlelist[i]);
4279  if (trackMapIter != pfParticleToTrackMap.end()) {
4280  lar_pandora::TrackVector pfParticleTracks = trackMapIter->second;
4281 
4282  if (pfParticleTracks.size() > 1)
4283  std::cerr << "Warning: there was more than one track found for PFParticle with ID " << pfparticlelist[i]->Self() << std::endl;
4284 
4285  if (pfParticleTracks.size() > 0) {
4286  PFParticleData.pfp_trackID[i] = pfParticleTracks.at(0)->ID();
4287  trackIDtoPFParticleIDMap.insert(std::make_pair(pfParticleTracks.at(0)->ID(), pfparticlelist[i]->Self()));
4288  }
4289  }
4290  else
4291  std::cerr << "Warning: there was no track found for track-like PFParticle with ID " << pfparticlelist[i]->Self() << std::endl;
4292  }
4293  else
4294  PFParticleData.pfp_isTrack[i] = 0;
4295 
4296  if (lar_pandora::LArPandoraHelper::IsShower(pfparticlelist[i])) {
4297  PFParticleData.pfp_isShower[i] = 1;
4298  // Set the shower ID.
4299  auto showerMapIter = pfParticleToShowerMap.find(pfparticlelist[i]);
4300  if (showerMapIter != pfParticleToShowerMap.end()) {
4301  lar_pandora::ShowerVector pfParticleShowers = showerMapIter->second;
4302 
4303  if (pfParticleShowers.size() > 1)
4304  std::cerr << "Warning: there was more than one shower found for PFParticle with ID " << pfparticlelist[i]->Self() << std::endl;
4305 
4306  if (pfParticleShowers.size() > 0) {
4307  PFParticleData.pfp_showerID[i] = pfParticleShowers.at(0)->ID();
4308  showerIDtoPFParticleIDMap.insert(std::make_pair(pfParticleShowers.at(0)->ID(), pfparticlelist[i]->Self()));
4309  }
4310  }
4311  else
4312  std::cerr << "Warning: there was no shower found for shower-like PFParticle with ID " << pfparticlelist[i]->Self() << std::endl;
4313  }
4314  else
4315  PFParticleData.pfp_isShower[i] = 0;
4316 
4317  // Set the cluster IDs.
4318  auto clusterMapIter = pfParticleToClusterMap.find(pfparticlelist[i]);
4319  if (clusterMapIter != pfParticleToClusterMap.end()) {
4320  lar_pandora::ClusterVector pfParticleClusters = clusterMapIter->second;
4321  PFParticleData.pfp_numClusters[i] = pfParticleClusters.size();
4322 
4323  for (size_t j = 0; j < pfParticleClusters.size(); ++j)
4324  PFParticleData.pfp_clusterIDs[i][j] = pfParticleClusters[j]->ID();
4325  }
4326  //else
4327  // std::cerr << "Warning: there were no clusters found for PFParticle with ID " << pfparticlelist[i]->Self() << std::endl;
4328  }
4329  } // if fSavePFParticleInfo
4330 
4331  if (fSaveShowerInfo){
4332 
4333  // fill data from all the shower algorithms
4334  for (size_t iShowerAlgo = 0; iShowerAlgo < NShowerAlgos; ++iShowerAlgo) {
4335  AnalysisTreeDataStruct::ShowerDataStruct& ShowerData
4336  = fData->GetShowerData(iShowerAlgo);
4337  std::vector<recob::Shower> const* pShowers = showerList[iShowerAlgo];
4338  art::Handle< std::vector<recob::Shower> > showerHandle = showerListHandle[iShowerAlgo];
4339 
4340  if (pShowers){
4341  FillShowers(ShowerData, *pShowers, fSavePFParticleInfo, showerIDtoPFParticleIDMap);
4342 
4343  if(fMVAPIDShowerModuleLabel[iShowerAlgo].size()){
4344  art::FindOneP<anab::MVAPIDResult> fmvapid(showerHandle, evt, fMVAPIDShowerModuleLabel[iShowerAlgo]);
4345  if(fmvapid.isValid()){
4346  for(unsigned int iShower=0;iShower<showerHandle->size();++iShower){
4347  const art::Ptr<anab::MVAPIDResult> pid = fmvapid.at(iShower);
4348  ShowerData.shwr_pidmvamu[iShower] = pid->mvaOutput.at("muon");
4349  ShowerData.shwr_pidmvae[iShower] = pid->mvaOutput.at("electron");
4350  ShowerData.shwr_pidmvapich[iShower] = pid->mvaOutput.at("pich");
4351  ShowerData.shwr_pidmvaphoton[iShower] = pid->mvaOutput.at("photon");
4352  ShowerData.shwr_pidmvapr[iShower] = pid->mvaOutput.at("proton");
4353  }
4354  } // fmvapid.isValid()
4355  }
4356  }
4357  else ShowerData.MarkMissing(fTree); // tree should reflect lack of data
4358  } // for iShowerAlgo
4359 
4360  } // if fSaveShowerInfo
4361 
4362  //track information for multiple trackers
4363  if (fSaveTrackInfo) {
4364  for (unsigned int iTracker=0; iTracker < NTrackers; ++iTracker){
4365  AnalysisTreeDataStruct::TrackDataStruct& TrackerData = fData->GetTrackerData(iTracker);
4366 
4367  size_t NTracks = tracklist[iTracker].size();
4368  // allocate enough space for this number of tracks (but at least for one of them!)
4369  TrackerData.SetMaxTracks(std::max(NTracks, (size_t) 1));
4370  TrackerData.Clear(); // clear all the data
4371 
4372  TrackerData.ntracks = (int) NTracks;
4373 
4374  // now set the tree addresses to the newly allocated memory;
4375  // this creates the tree branches in case they are not there yet
4376  SetTrackerAddresses(iTracker);
4377  if (NTracks > TrackerData.GetMaxTracks()) {
4378  // got this error? it might be a bug,
4379  // since we are supposed to have allocated enough space to fit all tracks
4380  mf::LogError("AnalysisTree:limits") << "event has " << NTracks
4381  << " " << fTrackModuleLabel[iTracker] << " tracks, only "
4382  << TrackerData.GetMaxTracks() << " stored in tree";
4383  }
4384 
4385  //call the track momentum algorithm that gives you momentum based on track range
4386  // - Should the minimal track length be 50 cm? The default of 100 has been used.
4387  trkf::TrackMomentumCalculator trkm{/*100.*/};
4388 
4389  for(size_t iTrk=0; iTrk < NTracks; ++iTrk){//loop over tracks
4390 
4391  //save t0 from reconstructed flash track matching for every track
4392  art::FindManyP<anab::T0> fmt0(trackListHandle[iTracker],evt,fFlashT0FinderLabel[iTracker]);
4393  if (fmt0.isValid()){
4394  if(fmt0.at(iTrk).size()>0){
4395  if(fmt0.at(iTrk).size()>1)
4396  std::cerr << "\n Warning : more than one cosmic tag per track in module! assigning the first tag to the track" << fFlashT0FinderLabel[iTracker];
4397  TrackerData.trkflashT0[iTrk] = fmt0.at(iTrk).at(0)->Time();
4398  }
4399  }
4400 
4401  //save t0 from reconstructed flash track matching for every track
4402  art::FindManyP<anab::T0> fmmct0(trackListHandle[iTracker],evt,fMCT0FinderLabel[iTracker]);
4403  if (fmmct0.isValid()){
4404  if(fmmct0.at(iTrk).size()>0){
4405  if(fmmct0.at(iTrk).size()>1)
4406  std::cerr << "\n Warning : more than one cosmic tag per track in module! assigning the first tag to the cluster" << fMCT0FinderLabel[iTracker];
4407  TrackerData.trktrueT0[iTrk] = fmmct0.at(iTrk).at(0)->Time();
4408  }
4409  }
4410 
4411  //Cosmic Tagger information
4412  art::FindManyP<anab::CosmicTag> fmct(trackListHandle[iTracker],evt,fCosmicTaggerAssocLabel[iTracker]);
4413  if (fmct.isValid()){
4414  TrackerData.trkncosmictags_tagger[iTrk] = fmct.at(iTrk).size();
4415  if (fmct.at(iTrk).size()>0){
4416  if(fmct.at(iTrk).size()>1)
4417  std::cerr << "\n Warning : more than one cosmic tag per track in module! assigning the first tag to the track" << fCosmicTaggerAssocLabel[iTracker];
4418  TrackerData.trkcosmicscore_tagger[iTrk] = fmct.at(iTrk).at(0)->CosmicScore();
4419  TrackerData.trkcosmictype_tagger[iTrk] = fmct.at(iTrk).at(0)->CosmicType();
4420  }
4421  }
4422 
4423  //Containment Tagger information
4424  art::FindManyP<anab::CosmicTag> fmcnt(trackListHandle[iTracker],evt,fContainmentTaggerAssocLabel[iTracker]);
4425  if (fmcnt.isValid()){
4426  TrackerData.trkncosmictags_containmenttagger[iTrk] = fmcnt.at(iTrk).size();
4427  if (fmcnt.at(iTrk).size()>0){
4428  if(fmcnt.at(iTrk).size()>1)
4429  std::cerr << "\n Warning : more than one containment tag per track in module! assigning the first tag to the track" << fContainmentTaggerAssocLabel[iTracker];
4430  TrackerData.trkcosmicscore_containmenttagger[iTrk] = fmcnt.at(iTrk).at(0)->CosmicScore();
4431  TrackerData.trkcosmictype_containmenttagger[iTrk] = fmcnt.at(iTrk).at(0)->CosmicType();
4432  }
4433  }
4434 
4435  //Flash match compatibility information
4436  //Unlike CosmicTagger, Flash match doesn't assign a cosmic tag for every track. For those tracks, AnalysisTree initializes them with -9999 or -99999
4437  art::FindManyP<anab::CosmicTag> fmbfm(trackListHandle[iTracker],evt,fFlashMatchAssocLabel[iTracker]);
4438  if (fmbfm.isValid()){
4439  TrackerData.trkncosmictags_flashmatch[iTrk] = fmbfm.at(iTrk).size();
4440  if (fmbfm.at(iTrk).size()>0){
4441  if(fmbfm.at(iTrk).size()>1)
4442  std::cerr << "\n Warning : more than one cosmic tag per track in module! assigning the first tag to the track" << fFlashMatchAssocLabel[iTracker];
4443  TrackerData.trkcosmicscore_flashmatch[iTrk] = fmbfm.at(iTrk).at(0)->CosmicScore();
4444  TrackerData.trkcosmictype_flashmatch[iTrk] = fmbfm.at(iTrk).at(0)->CosmicType();
4445  //std::cout<<"\n"<<evt.event()<<"\t"<<iTrk<<"\t"<<fmbfm.at(iTrk).at(0)->CosmicScore()<<"\t"<<fmbfm.at(iTrk).at(0)->CosmicType();
4446  }
4447  }
4448 
4449  art::Ptr<recob::Track> ptrack(trackListHandle[iTracker], iTrk);
4450  const recob::Track& track = *ptrack;
4451 
4452  TVector3 pos, dir_start, dir_end, end;
4453 
4454  double tlen = 0., mom = 0.;
4455  int TrackID = -1;
4456 
4457  int ntraj = track.NumberTrajectoryPoints();
4458  if (ntraj > 0) {
4459  pos = track.Vertex<TVector3>();
4460  dir_start = track.VertexDirection<TVector3>();
4461  dir_end = track.EndDirection<TVector3>();
4462  end = track.End<TVector3>();
4463 
4464  tlen = track.Length();
4465  if(track.NumberTrajectoryPoints() > 0)
4466  mom = track.VertexMomentum();
4467  // fill non-bezier-track reco branches
4468  TrackID = track.ID();
4469 
4470  double theta_xz = std::atan2(dir_start.X(), dir_start.Z());
4471  double theta_yz = std::atan2(dir_start.Y(), dir_start.Z());
4472  double dpos = bdist(pos); // FIXME - Passing an uncorrected position....
4473  double dend = bdist(end); // FIXME - Passing an uncorrected position....
4474 
4475  TrackerData.trkId[iTrk] = TrackID;
4476  TrackerData.trkstartx[iTrk] = pos.X();
4477  TrackerData.trkstarty[iTrk] = pos.Y();
4478  TrackerData.trkstartz[iTrk] = pos.Z();
4479  TrackerData.trkstartd[iTrk] = dpos;
4480  TrackerData.trkendx[iTrk] = end.X();
4481  TrackerData.trkendy[iTrk] = end.Y();
4482  TrackerData.trkendz[iTrk] = end.Z();
4483  TrackerData.trkendd[iTrk] = dend;
4484  TrackerData.trktheta[iTrk] = dir_start.Theta();
4485  TrackerData.trkphi[iTrk] = dir_start.Phi();
4486  TrackerData.trkstartdcosx[iTrk] = dir_start.X();
4487  TrackerData.trkstartdcosy[iTrk] = dir_start.Y();
4488  TrackerData.trkstartdcosz[iTrk] = dir_start.Z();
4489  TrackerData.trkenddcosx[iTrk] = dir_end.X();
4490  TrackerData.trkenddcosy[iTrk] = dir_end.Y();
4491  TrackerData.trkenddcosz[iTrk] = dir_end.Z();
4492  TrackerData.trkthetaxz[iTrk] = theta_xz;
4493  TrackerData.trkthetayz[iTrk] = theta_yz;
4494  TrackerData.trkmom[iTrk] = mom;
4495  TrackerData.trklen[iTrk] = tlen;
4496  TrackerData.trkmomrange[iTrk] = trkm.GetTrackMomentum(tlen,13);
4497  //TrackerData.trkmommschi2[iTrk] = trkm.GetMomentumMultiScatterChi2(ptrack);
4498  //TrackerData.trkmommsllhd[iTrk] = trkm.GetMomentumMultiScatterLLHD(ptrack);
4499 
4500  if (fSavePFParticleInfo) {
4501  auto mapIter = trackIDtoPFParticleIDMap.find(TrackID);
4502  if (mapIter != trackIDtoPFParticleIDMap.end()) {
4503  // This track has a corresponding PFParticle.
4504  TrackerData.trkhasPFParticle[iTrk] = 1;
4505  TrackerData.trkPFParticleID[iTrk] = mapIter->second;
4506  }
4507  else
4508  TrackerData.trkhasPFParticle[iTrk] = 0;
4509  }
4510 
4511  } // if we have trajectory
4512 
4513  // find vertices associated with this track
4514  /*
4515  art::FindMany<recob::Vertex> fmvtx(trackListHandle[iTracker], evt, fVertexModuleLabel[iTracker]);
4516  if(fmvtx.isValid()) {
4517  std::vector<const recob::Vertex*> verts = fmvtx.at(iTrk);
4518  // should have two at most
4519  for(size_t ivx = 0; ivx < verts.size(); ++ivx) {
4520  verts[ivx]->XYZ(xyz);
4521  // find the vertex in TrackerData to get the index
4522  short theVtx = -1;
4523  for(short jvx = 0; jvx < TrackerData.nvtx; ++jvx) {
4524  if(TrackerData.vtx[jvx][2] == xyz[2]) {
4525  theVtx = jvx;
4526  break;
4527  }
4528  } // jvx
4529  // decide if it should be assigned to the track Start or End.
4530  // A simple dz test should suffice
4531  if(fabs(xyz[2] - TrackerData.trkstartz[iTrk]) <
4532  fabs(xyz[2] - TrackerData.trkendz[iTrk])) {
4533  TrackerData.trksvtxid[iTrk] = theVtx;
4534  } else {
4535  TrackerData.trkevtxid[iTrk] = theVtx;
4536  }
4537  } // vertices
4538  } // fmvtx.isValid()
4539  */
4540 
4541 
4542  /* //commented out because now have several Vertices
4543  Float_t minsdist = 10000;
4544  Float_t minedist = 10000;
4545  for (int ivx = 0; ivx < NVertices && ivx < kMaxVertices; ++ivx){
4546  Float_t sdist = sqrt(pow(TrackerData.trkstartx[iTrk]-VertexData.vtxx[ivx],2)+
4547  pow(TrackerData.trkstarty[iTrk]-VertexData.vtxy[ivx],2)+
4548  pow(TrackerData.trkstartz[iTrk]-VertexData.vtxz[ivx],2));
4549  Float_t edist = sqrt(pow(TrackerData.trkendx[iTrk]-VertexData.vtxx[ivx],2)+
4550  pow(TrackerData.trkendy[iTrk]-VertexData.vtxy[ivx],2)+
4551  pow(TrackerData.trkendz[iTrk]-VertexData.vtxz[ivx],2));
4552  if (sdist<minsdist){
4553  minsdist = sdist;
4554  if (minsdist<10) TrackerData.trksvtxid[iTrk] = ivx;
4555  }
4556  if (edist<minedist){
4557  minedist = edist;
4558  if (minedist<10) TrackerData.trkevtxid[iTrk] = ivx;
4559  }
4560  }*/
4561 
4562  // find particle ID info
4563  /*
4564  // Note from Jake Calcutt: There has been a breaking change in the definition
4565  // of the anab::ParticleID class. If you are using this class and want this info
4566  // in this tree, these must be updated
4567  art::FindMany<anab::ParticleID> fmpid(trackListHandle[iTracker], evt, fParticleIDModuleLabel[iTracker]);
4568  if(fmpid.isValid()) {
4569  std::vector<const anab::ParticleID*> pids = fmpid.at(iTrk);
4570  //if(pids.size() > 1) {
4571  //mf::LogError("AnalysisTree:limits")
4572  //<< "the " << fTrackModuleLabel[iTracker] << " track #" << iTrk
4573  //<< " has " << pids.size()
4574  //<< " set of ParticleID variables. Only one stored in the tree";
4575  //}
4576  for (size_t ipid = 0; ipid < pids.size(); ++ipid){
4577  if (!pids[ipid]->PlaneID().isValid) continue;
4578  int planenum = pids[ipid]->PlaneID().Plane;
4579  if (planenum<0||planenum>2) continue;
4580  TrackerData.trkpidpdg[iTrk][planenum] = pids[ipid]->Pdg();
4581  TrackerData.trkpidchi[iTrk][planenum] = pids[ipid]->MinChi2();
4582  TrackerData.trkpidchipr[iTrk][planenum] = pids[ipid]->Chi2Proton();
4583  TrackerData.trkpidchika[iTrk][planenum] = pids[ipid]->Chi2Kaon();
4584  TrackerData.trkpidchipi[iTrk][planenum] = pids[ipid]->Chi2Pion();
4585  TrackerData.trkpidchimu[iTrk][planenum] = pids[ipid]->Chi2Muon();
4586  TrackerData.trkpidpida[iTrk][planenum] = pids[ipid]->PIDA();
4587  }
4588  } // fmpid.isValid()
4589  */
4590  if(fMVAPIDTrackModuleLabel[iTracker].size()){
4591  art::FindOneP<anab::MVAPIDResult> fmvapid(trackListHandle[iTracker], evt, fMVAPIDTrackModuleLabel[iTracker]);
4592  if(fmvapid.isValid()) {
4593  const art::Ptr<anab::MVAPIDResult> pid = fmvapid.at(iTrk);
4594  TrackerData.trkpidmvamu[iTrk] = pid->mvaOutput.at("muon");
4595  TrackerData.trkpidmvae[iTrk] = pid->mvaOutput.at("electron");
4596  TrackerData.trkpidmvapich[iTrk] = pid->mvaOutput.at("pich");
4597  TrackerData.trkpidmvaphoton[iTrk] = pid->mvaOutput.at("photon");
4598  TrackerData.trkpidmvapr[iTrk] = pid->mvaOutput.at("proton");
4599  } // fmvapid.isValid()
4600  }
4601  art::FindMany<anab::Calorimetry> fmcal(trackListHandle[iTracker], evt, fCalorimetryModuleLabel[iTracker]);
4602  if (fmcal.isValid()){
4603  std::vector<const anab::Calorimetry*> calos = fmcal.at(iTrk);
4604  if (calos.size() > TrackerData.GetMaxPlanesPerTrack(iTrk)) {
4605  // if you get this message, there is probably a bug somewhere since
4606  // the calorimetry planes should be 3.
4607  mf::LogError("AnalysisTree:limits")
4608  << "the " << fTrackModuleLabel[iTracker] << " track #" << iTrk
4609  << " has " << calos.size() << " planes for calorimetry , only "
4610  << TrackerData.GetMaxPlanesPerTrack(iTrk) << " stored in tree";
4611  }
4612  for (size_t ical = 0; ical<calos.size(); ++ical){
4613  if (!calos[ical]) continue;
4614  if (!calos[ical]->PlaneID().isValid) continue;
4615  int planenum = calos[ical]->PlaneID().Plane;
4616  if (planenum<0||planenum>2) continue;
4617  TrackerData.trkke[iTrk][planenum] = calos[ical]->KineticEnergy();
4618  TrackerData.trkrange[iTrk][planenum] = calos[ical]->Range();
4619  //For now make the second argument as 13 for muons.
4620  TrackerData.trkpitchc[iTrk][planenum]= calos[ical] -> TrkPitchC();
4621  const size_t NHits = calos[ical] -> dEdx().size();
4622  TrackerData.ntrkhits[iTrk][planenum] = (int) NHits;
4623  if (NHits > TrackerData.GetMaxHitsPerTrack(iTrk, planenum)) {
4624  // if you get this error, you'll have to increase kMaxTrackHits
4625  mf::LogError("AnalysisTree:limits")
4626  << "the " << fTrackModuleLabel[iTracker] << " track #" << iTrk
4627  << " has " << NHits << " hits on calorimetry plane #" << planenum
4628  <<", only "
4629  << TrackerData.GetMaxHitsPerTrack(iTrk, planenum) << " stored in tree";
4630  }
4631  if (!isCosmics){
4632  for(size_t iTrkHit = 0; iTrkHit < NHits && iTrkHit < TrackerData.GetMaxHitsPerTrack(iTrk, planenum); ++iTrkHit) {
4633  TrackerData.trkdedx[iTrk][planenum][iTrkHit] = (calos[ical] -> dEdx())[iTrkHit];
4634  TrackerData.trkdqdx[iTrk][planenum][iTrkHit] = (calos[ical] -> dQdx())[iTrkHit];
4635  TrackerData.trkresrg[iTrk][planenum][iTrkHit] = (calos[ical] -> ResidualRange())[iTrkHit];
4636  TrackerData.trktpc[iTrk][planenum][iTrkHit] = (calos[ical] -> PlaneID()).TPC;
4637  const auto& TrkPos = (calos[ical] -> XYZ())[iTrkHit];
4638  auto& TrkXYZ = TrackerData.trkxyz[iTrk][planenum][iTrkHit];
4639  TrkXYZ[0] = TrkPos.X();
4640  TrkXYZ[1] = TrkPos.Y();
4641  TrkXYZ[2] = TrkPos.Z();
4642  } // for track hits
4643  }
4644  } // for calorimetry info
4645  if(TrackerData.ntrkhits[iTrk][0] > TrackerData.ntrkhits[iTrk][1] && TrackerData.ntrkhits[iTrk][0] > TrackerData.ntrkhits[iTrk][2]) TrackerData.trkpidbestplane[iTrk] = 0;
4646  else if(TrackerData.ntrkhits[iTrk][1] > TrackerData.ntrkhits[iTrk][0] && TrackerData.ntrkhits[iTrk][1] > TrackerData.ntrkhits[iTrk][2]) TrackerData.trkpidbestplane[iTrk] = 1;
4647  else if(TrackerData.ntrkhits[iTrk][2] > TrackerData.ntrkhits[iTrk][0] && TrackerData.ntrkhits[iTrk][2] > TrackerData.ntrkhits[iTrk][1]) TrackerData.trkpidbestplane[iTrk] = 2;
4648  else if(TrackerData.ntrkhits[iTrk][2] == TrackerData.ntrkhits[iTrk][0] && TrackerData.ntrkhits[iTrk][2] > TrackerData.ntrkhits[iTrk][1]) TrackerData.trkpidbestplane[iTrk] = 2;
4649  else if(TrackerData.ntrkhits[iTrk][2] == TrackerData.ntrkhits[iTrk][1] && TrackerData.ntrkhits[iTrk][2] > TrackerData.ntrkhits[iTrk][0]) TrackerData.trkpidbestplane[iTrk] = 2;
4650  else if(TrackerData.ntrkhits[iTrk][1] == TrackerData.ntrkhits[iTrk][0] && TrackerData.ntrkhits[iTrk][1] > TrackerData.ntrkhits[iTrk][2]) TrackerData.trkpidbestplane[iTrk] = 0;
4651  else if(TrackerData.ntrkhits[iTrk][1] == TrackerData.ntrkhits[iTrk][0] && TrackerData.ntrkhits[iTrk][1] == TrackerData.ntrkhits[iTrk][2]) TrackerData.trkpidbestplane[iTrk] = 2;
4652 
4653  // FIXME - Do i want to add someway to work out the best TPC???....
4654  } // if has calorimetry info
4655 
4656  //track truth information
4657  if (isMC){
4658  //get the hits on each plane
4659  art::FindManyP<recob::Hit> fmht(trackListHandle[iTracker], evt, fTrackModuleLabel[iTracker]);
4660  std::vector< art::Ptr<recob::Hit> > allHits = fmht.at(iTrk);
4661  std::vector< art::Ptr<recob::Hit> > hits[kNplanes];
4662 
4663  for(size_t ah = 0; ah < allHits.size(); ++ah){
4664  if (/* allHits[ah]->WireID().Plane >= 0 && */ // always true
4665  allHits[ah]->WireID().Plane < 3){
4666  hits[allHits[ah]->WireID().Plane].push_back(allHits[ah]);
4667  }
4668  }
4669  for (size_t ipl = 0; ipl < 3; ++ipl){
4670  double maxe = 0;
4671  HitsPurity(clockData, hits[ipl],TrackerData.trkidtruth[iTrk][ipl],TrackerData.trkpurtruth[iTrk][ipl],maxe);
4672  //std::cout<<"\n"<<iTracker<<"\t"<<iTrk<<"\t"<<ipl<<"\t"<<trkidtruth[iTracker][iTrk][ipl]<<"\t"<<trkpurtruth[iTracker][iTrk][ipl]<<"\t"<<maxe;
4673  if (TrackerData.trkidtruth[iTrk][ipl]>0){
4674  const art::Ptr<simb::MCTruth> mc = pi_serv->TrackIdToMCTruth_P(TrackerData.trkidtruth[iTrk][ipl]);
4675  TrackerData.trkorigin[iTrk][ipl] = mc->Origin();
4676  const simb::MCParticle *particle = pi_serv->TrackIdToParticle_P(TrackerData.trkidtruth[iTrk][ipl]);
4677  double tote = 0;
4678  const std::vector<const sim::IDE*> vide=bt_serv->TrackIdToSimIDEs_Ps(TrackerData.trkidtruth[iTrk][ipl]);
4679  for (auto ide: vide) {
4680  tote += ide->energy;
4681  }
4682  TrackerData.trkpdgtruth[iTrk][ipl] = particle->PdgCode();
4683  TrackerData.trkefftruth[iTrk][ipl] = maxe/(tote/kNplanes); //tote include both induction and collection energies
4684  //std::cout<<"\n"<<trkpdgtruth[iTracker][iTrk][ipl]<<"\t"<<trkefftruth[iTracker][iTrk][ipl];
4685  }
4686  }
4687 
4688  double maxe = 0;
4689  HitsPurity(clockData, allHits,TrackerData.trkg4id[iTrk],TrackerData.trkpurity[iTrk],maxe);
4690  if (TrackerData.trkg4id[iTrk]>0){
4691  const art::Ptr<simb::MCTruth> mc = pi_serv->TrackIdToMCTruth_P(TrackerData.trkg4id[iTrk]);
4692  TrackerData.trkorig[iTrk] = mc->Origin();
4693  }
4694  if (allHits.size()){
4695  std::vector<art::Ptr<recob::Hit> > all_hits;
4697  float totenergy = 0;
4698  if (evt.get(allHits[0].id(), hithandle)){
4699  art::fill_ptr_vector(all_hits, hithandle);
4700  for(size_t h = 0; h < all_hits.size(); ++h){
4701 
4702  art::Ptr<recob::Hit> hit = all_hits[h];
4703  std::vector<sim::IDE*> ides;
4704  //bt_serv->HitToSimIDEs(hit,ides);
4705  std::vector<sim::TrackIDE> eveIDs = bt_serv->HitToEveTrackIDEs(clockData, hit);
4706 
4707  for(size_t e = 0; e < eveIDs.size(); ++e){
4708  //std::cout<<h<<" "<<e<<" "<<eveIDs[e].trackID<<" "<<eveIDs[e].energy<<" "<<eveIDs[e].energyFrac<<std::endl;
4709  if (eveIDs[e].trackID==TrackerData.trkg4id[iTrk]) totenergy += eveIDs[e].energy;
4710  }
4711  }
4712  }
4713  if (totenergy) TrackerData.trkcompleteness[iTrk] = maxe/totenergy;
4714  }
4715  }//end if (isMC)
4716  }//end loop over track
4717  }//end loop over track module labels
4718  }// end (fSaveTrackInfo)
4719 
4720  /*trkf::TrackMomentumCalculator trkm;
4721  std::cout<<"\t"<<trkm.GetTrackMomentum(200,2212)<<"\t"<<trkm.GetTrackMomentum(-10, 13)<<"\t"<<trkm.GetTrackMomentum(300,-19)<<"\n";
4722  */
4723 
4724  //Save Vertex information for multiple algorithms
4725  if (fSaveVertexInfo){
4726  for (unsigned int iVertexAlg=0; iVertexAlg < NVertexAlgos; ++iVertexAlg){
4727  AnalysisTreeDataStruct::VertexDataStruct& VertexData = fData->GetVertexData(iVertexAlg);
4728 
4729  size_t NVertices = vertexlist[iVertexAlg].size();
4730 
4731  VertexData.SetMaxVertices(std::max(NVertices, (size_t) 1));
4732  VertexData.Clear(); // clear all the data
4733 
4734  VertexData.nvtx = (short) NVertices;
4735 
4736  // now set the tree addresses to the newly allocated memory;
4737  // this creates the tree branches in case they are not there yet
4738  SetVertexAddresses(iVertexAlg);
4739  if (NVertices > VertexData.GetMaxVertices()) {
4740  // got this error? it might be a bug,
4741  // since we are supposed to have allocated enough space to fit all tracks
4742  mf::LogError("AnalysisTree:limits") << "event has " << NVertices
4743  << " " << fVertexModuleLabel[iVertexAlg] << " tracks, only "
4744  << VertexData.GetMaxVertices() << " stored in tree";
4745  }
4746 
4747  for (size_t i = 0; i < NVertices && i < kMaxVertices ; ++i){//loop over hits
4748  VertexData.vtxId[i] = vertexlist[iVertexAlg][i]->ID();
4749  Double_t xyz[3] = {};
4750  vertexlist[iVertexAlg][i] -> XYZ(xyz);
4751  VertexData.vtxx[i] = xyz[0];
4752  VertexData.vtxy[i] = xyz[1];
4753  VertexData.vtxz[i] = xyz[2];
4754 
4755  if (fSavePFParticleInfo) {
4756  auto mapIter = vertexIDtoPFParticleIDMap.find(vertexlist[iVertexAlg][i]->ID());
4757  if (mapIter != vertexIDtoPFParticleIDMap.end()) {
4758  // This vertex has a corresponding PFParticle.
4759  VertexData.vtxhasPFParticle[i] = 1;
4760  VertexData.vtxPFParticleID[i] = mapIter->second;
4761  }
4762  else
4763  VertexData.vtxhasPFParticle[i] = 0;
4764  }
4765 
4766  // find PFParticle ID info
4767  art::FindMany<recob::PFParticle> fmPFParticle(vertexListHandle[iVertexAlg], evt, fPFParticleModuleLabel);
4768  if(fmPFParticle.isValid()) {
4769  std::vector<const recob::PFParticle*> pfparticles = fmPFParticle.at(i);
4770  if(pfparticles.size() > 1)
4771  std::cerr << "Warning: more than one associated PFParticle found for a vertex. Only one stored in tree." << std::endl;
4772  if (pfparticles.size() == 0)
4773  VertexData.vtxhasPFParticle[i] = 0;
4774  else {
4775  VertexData.vtxhasPFParticle[i] = 1;
4776  VertexData.vtxPFParticleID[i] = pfparticles.at(0)->Self();
4777  }
4778  } // fmPFParticle.isValid()
4779  }
4780  }
4781  }
4782 
4783  //mc truth information
4784  if (isMC){
4785 
4786  // Find the simb::MCFlux objects corresponding to
4787  // each simb::MCTruth object made by the generator with
4788  // the label fGenieGenModuleLabel
4789  art::FindOne<simb::MCFlux> find_mcflux(mctruthListHandle,
4790  evt, fGenieGenModuleLabel);
4791 
4792  if (fSaveCryInfo){
4793  //store cry (cosmic generator information)
4794  fData->mcevts_truthcry = mclistcry.size();
4795  fData->cry_no_primaries = nCryPrimaries;
4796  //fData->cry_no_primaries;
4797  for(Int_t iPartc = 0; iPartc < mctruthcry->NParticles(); ++iPartc){
4798  const simb::MCParticle& partc(mctruthcry->GetParticle(iPartc));
4799  fData->cry_primaries_pdg[iPartc]=partc.PdgCode();
4800  fData->cry_Eng[iPartc]=partc.E();
4801  fData->cry_Px[iPartc]=partc.Px();
4802  fData->cry_Py[iPartc]=partc.Py();
4803  fData->cry_Pz[iPartc]=partc.Pz();
4804  fData->cry_P[iPartc]=partc.P();
4805  fData->cry_StartPointx[iPartc] = partc.Vx();
4806  fData->cry_StartPointy[iPartc] = partc.Vy();
4807  fData->cry_StartPointz[iPartc] = partc.Vz();
4808  fData->cry_StartPointt[iPartc] = partc.T();
4809  fData->cry_status_code[iPartc]=partc.StatusCode();
4810  fData->cry_mass[iPartc]=partc.Mass();
4811  fData->cry_trackID[iPartc]=partc.TrackId();
4812  fData->cry_ND[iPartc]=partc.NumberDaughters();
4813  fData->cry_mother[iPartc]=partc.Mother();
4814  } // for cry particles
4815  }// end fSaveCryInfo
4816 
4817  // Save the protoDUNE beam generator information
4818  if(fSaveProtoInfo){
4819  fData->proto_no_primaries = nProtoPrimaries;
4820  for(Int_t iPartp = 0; iPartp < nProtoPrimaries; ++iPartp){
4821  const simb::MCParticle& partp(mctruthproto->GetParticle(iPartp));
4822 
4823  fData->proto_isGoodParticle[iPartp] = (partp.Process() == "primary");
4824  fData->proto_vx[iPartp] = partp.Vx();
4825  fData->proto_vy[iPartp] = partp.Vy();
4826  fData->proto_vz[iPartp] = partp.Vz();
4827  fData->proto_t[iPartp] = partp.T();
4828  fData->proto_px[iPartp] = partp.Px();
4829  fData->proto_py[iPartp] = partp.Py();
4830  fData->proto_pz[iPartp] = partp.Pz();
4831  fData->proto_momentum[iPartp] = partp.P();
4832  fData->proto_energy[iPartp] = partp.E();
4833  fData->proto_pdg[iPartp] = partp.PdgCode();
4834  // We will deal with the matching to GEANT later
4835  }
4836  }
4837 
4838  //save neutrino interaction information
4839  fData->mcevts_truth = mclist.size();
4840  if (fData->mcevts_truth > 0){//at least one mc record
4841  if (fSaveGenieInfo){
4842  int neutrino_i = 0;
4843  for(unsigned int iList = 0; (iList < mclist.size()) && (neutrino_i < kMaxTruth) ; ++iList){
4844  if (mclist[iList]->NeutrinoSet()){
4845  fData->nuPDG_truth[neutrino_i] = mclist[iList]->GetNeutrino().Nu().PdgCode();
4846  fData->ccnc_truth[neutrino_i] = mclist[iList]->GetNeutrino().CCNC();
4847  fData->mode_truth[neutrino_i] = mclist[iList]->GetNeutrino().Mode();
4848  fData->Q2_truth[neutrino_i] = mclist[iList]->GetNeutrino().QSqr();
4849  fData->W_truth[neutrino_i] = mclist[iList]->GetNeutrino().W();
4850  fData->X_truth[neutrino_i] = mclist[iList]->GetNeutrino().X();
4851  fData->Y_truth[neutrino_i] = mclist[iList]->GetNeutrino().Y();
4852  fData->hitnuc_truth[neutrino_i] = mclist[iList]->GetNeutrino().HitNuc();
4853  fData->enu_truth[neutrino_i] = mclist[iList]->GetNeutrino().Nu().E();
4854  fData->nuvtxx_truth[neutrino_i] = mclist[iList]->GetNeutrino().Nu().Vx();
4855  fData->nuvtxy_truth[neutrino_i] = mclist[iList]->GetNeutrino().Nu().Vy();
4856  fData->nuvtxz_truth[neutrino_i] = mclist[iList]->GetNeutrino().Nu().Vz();
4857  if (mclist[iList]->GetNeutrino().Nu().P()){
4858  fData->nu_dcosx_truth[neutrino_i] = mclist[iList]->GetNeutrino().Nu().Px()/mclist[iList]->GetNeutrino().Nu().P();
4859  fData->nu_dcosy_truth[neutrino_i] = mclist[iList]->GetNeutrino().Nu().Py()/mclist[iList]->GetNeutrino().Nu().P();
4860  fData->nu_dcosz_truth[neutrino_i] = mclist[iList]->GetNeutrino().Nu().Pz()/mclist[iList]->GetNeutrino().Nu().P();
4861  }
4862  fData->lep_mom_truth[neutrino_i] = mclist[iList]->GetNeutrino().Lepton().P();
4863  if (mclist[iList]->GetNeutrino().Lepton().P()){
4864  fData->lep_dcosx_truth[neutrino_i] = mclist[iList]->GetNeutrino().Lepton().Px()/mclist[iList]->GetNeutrino().Lepton().P();
4865  fData->lep_dcosy_truth[neutrino_i] = mclist[iList]->GetNeutrino().Lepton().Py()/mclist[iList]->GetNeutrino().Lepton().P();
4866  fData->lep_dcosz_truth[neutrino_i] = mclist[iList]->GetNeutrino().Lepton().Pz()/mclist[iList]->GetNeutrino().Lepton().P();
4867  }
4868 
4869  //flux information
4870  //
4871  // Double-check that a simb::MCFlux object is associated with the
4872  // current simb::MCTruth object. For GENIE events, these should
4873  // always accompany each other. Other generators (e.g., MARLEY) may
4874  // create simb::MCTruth objects without corresponding simb::MCFlux
4875  // objects. -- S. Gardiner
4876  if (find_mcflux.isValid()) {
4877  auto flux_maybe_ref = find_mcflux.at(iList);
4878  if (flux_maybe_ref.isValid()) {
4879  auto flux_ref = flux_maybe_ref.ref();
4880  fData->vx_flux[neutrino_i] = flux_ref.fvx;
4881  fData->vy_flux[neutrino_i] = flux_ref.fvy;
4882  fData->vz_flux[neutrino_i] = flux_ref.fvz;
4883  fData->pdpx_flux[neutrino_i] = flux_ref.fpdpx;
4884  fData->pdpy_flux[neutrino_i] = flux_ref.fpdpy;
4885  fData->pdpz_flux[neutrino_i] = flux_ref.fpdpz;
4886  fData->ppdxdz_flux[neutrino_i] = flux_ref.fppdxdz;
4887  fData->ppdydz_flux[neutrino_i] = flux_ref.fppdydz;
4888  fData->pppz_flux[neutrino_i] = flux_ref.fpppz;
4889 
4890  fData->ptype_flux[neutrino_i] = flux_ref.fptype;
4891  fData->ppvx_flux[neutrino_i] = flux_ref.fppvx;
4892  fData->ppvy_flux[neutrino_i] = flux_ref.fppvy;
4893  fData->ppvz_flux[neutrino_i] = flux_ref.fppvz;
4894  fData->muparpx_flux[neutrino_i] = flux_ref.fmuparpx;
4895  fData->muparpy_flux[neutrino_i] = flux_ref.fmuparpy;
4896  fData->muparpz_flux[neutrino_i] = flux_ref.fmuparpz;
4897  fData->mupare_flux[neutrino_i] = flux_ref.fmupare;
4898 
4899  fData->tgen_flux[neutrino_i] = flux_ref.ftgen;
4900  fData->tgptype_flux[neutrino_i] = flux_ref.ftgptype;
4901  fData->tgppx_flux[neutrino_i] = flux_ref.ftgppx;
4902  fData->tgppy_flux[neutrino_i] = flux_ref.ftgppy;
4903  fData->tgppz_flux[neutrino_i] = flux_ref.ftgppz;
4904  fData->tprivx_flux[neutrino_i] = flux_ref.ftprivx;
4905  fData->tprivy_flux[neutrino_i] = flux_ref.ftprivy;
4906  fData->tprivz_flux[neutrino_i] = flux_ref.ftprivz;
4907 
4908  fData->dk2gen_flux[neutrino_i] = flux_ref.fdk2gen;
4909  fData->gen2vtx_flux[neutrino_i] = flux_ref.fgen2vtx;
4910 
4911  fData->tpx_flux[neutrino_i] = flux_ref.ftpx;
4912  fData->tpy_flux[neutrino_i] = flux_ref.ftpy;
4913  fData->tpz_flux[neutrino_i] = flux_ref.ftpz;
4914  fData->tptype_flux[neutrino_i] = flux_ref.ftptype;
4915  } // flux_maybe_ref.isValid()
4916  } // find_mcflux.isValid()
4917  neutrino_i++;
4918  }//mclist is NeutrinoSet()
4919  }//loop over mclist
4920 
4921  if (mctruth->NeutrinoSet()){
4922  //genie particles information
4923  fData->genie_no_primaries = mctruth->NParticles();
4924 
4925  size_t StoreParticles = std::min((size_t) fData->genie_no_primaries, fData->GetMaxGeniePrimaries());
4926  if (fData->genie_no_primaries > (int) StoreParticles) {
4927  // got this error? it might be a bug,
4928  // since the structure should have enough room for everything
4929  mf::LogError("AnalysisTree:limits") << "event has "
4930  << fData->genie_no_primaries << " MC particles, only "
4931  << StoreParticles << " stored in tree";
4932  }
4933  for(size_t iPart = 0; iPart < StoreParticles; ++iPart){
4934  const simb::MCParticle& part(mctruth->GetParticle(iPart));
4935  fData->genie_primaries_pdg[iPart]=part.PdgCode();
4936  fData->genie_Eng[iPart]=part.E();
4937  fData->genie_Px[iPart]=part.Px();
4938  fData->genie_Py[iPart]=part.Py();
4939  fData->genie_Pz[iPart]=part.Pz();
4940  fData->genie_P[iPart]=part.P();
4941  fData->genie_status_code[iPart]=part.StatusCode();
4942  fData->genie_mass[iPart]=part.Mass();
4943  fData->genie_trackID[iPart]=part.TrackId();
4944  fData->genie_ND[iPart]=part.NumberDaughters();
4945  fData->genie_mother[iPart]=part.Mother();
4946  } // for particle
4947  //const simb::MCNeutrino& nu(mctruth->GetNeutrino());
4948  } //if neutrino set
4949  }// end (fSaveGenieInfo)
4950 
4951  //Extract MC Shower information and fill the Shower branches
4952  if (fSaveMCShowerInfo){
4953  fData->no_mcshowers = nMCShowers;
4954  size_t shwr = 0;
4955  for(std::vector<sim::MCShower>::const_iterator imcshwr = mcshowerh->begin();
4956  imcshwr != mcshowerh->end(); ++imcshwr) {
4957  const sim::MCShower& mcshwr = *imcshwr;
4958  fData->mcshwr_origin[shwr] = mcshwr.Origin();
4959  fData->mcshwr_pdg[shwr] = mcshwr.PdgCode();
4960  fData->mcshwr_TrackId[shwr] = mcshwr.TrackID();
4961  fData->mcshwr_Process[shwr] = mcshwr.Process();
4962  fData->mcshwr_startX[shwr] = mcshwr.Start().X();
4963  fData->mcshwr_startY[shwr] = mcshwr.Start().Y();
4964  fData->mcshwr_startZ[shwr] = mcshwr.Start().Z();
4965  fData->mcshwr_endX[shwr] = mcshwr.End().X();
4966  fData->mcshwr_endY[shwr] = mcshwr.End().Y();
4967  fData->mcshwr_endZ[shwr] = mcshwr.End().Z();
4968  if (mcshwr.DetProfile().E()!= 0){
4969  fData->mcshwr_isEngDeposited[shwr] = 1;
4970  fData->mcshwr_CombEngX[shwr] = mcshwr.DetProfile().X();
4971  fData->mcshwr_CombEngY[shwr] = mcshwr.DetProfile().Y();
4972  fData->mcshwr_CombEngZ[shwr] = mcshwr.DetProfile().Z();
4973  fData->mcshwr_CombEngPx[shwr] = mcshwr.DetProfile().Px();
4974  fData->mcshwr_CombEngPy[shwr] = mcshwr.DetProfile().Py();
4975  fData->mcshwr_CombEngPz[shwr] = mcshwr.DetProfile().Pz();
4976  fData->mcshwr_CombEngE[shwr] = mcshwr.DetProfile().E();
4977  fData->mcshwr_dEdx[shwr] = mcshwr.dEdx();
4978  fData->mcshwr_StartDirX[shwr] = mcshwr.StartDir().X();
4979  fData->mcshwr_StartDirY[shwr] = mcshwr.StartDir().Y();
4980  fData->mcshwr_StartDirZ[shwr] = mcshwr.StartDir().Z();
4981  }
4982  else
4983  fData->mcshwr_isEngDeposited[shwr] = 0;
4984  fData->mcshwr_Motherpdg[shwr] = mcshwr.MotherPdgCode();
4985  fData->mcshwr_MotherTrkId[shwr] = mcshwr.MotherTrackID();
4986  fData->mcshwr_MotherProcess[shwr] = mcshwr.MotherProcess();
4987  fData->mcshwr_MotherstartX[shwr] = mcshwr.MotherStart().X();
4988  fData->mcshwr_MotherstartY[shwr] = mcshwr.MotherStart().Y();
4989  fData->mcshwr_MotherstartZ[shwr] = mcshwr.MotherStart().Z();
4990  fData->mcshwr_MotherendX[shwr] = mcshwr.MotherEnd().X();
4991  fData->mcshwr_MotherendY[shwr] = mcshwr.MotherEnd().Y();
4992  fData->mcshwr_MotherendZ[shwr] = mcshwr.MotherEnd().Z();
4993  fData->mcshwr_Ancestorpdg[shwr] = mcshwr.AncestorPdgCode();
4994  fData->mcshwr_AncestorTrkId[shwr] = mcshwr.AncestorTrackID();
4995  fData->mcshwr_AncestorProcess[shwr] = mcshwr.AncestorProcess();
4996  fData->mcshwr_AncestorstartX[shwr] = mcshwr.AncestorStart().X();
4997  fData->mcshwr_AncestorstartY[shwr] = mcshwr.AncestorStart().Y();
4998  fData->mcshwr_AncestorstartZ[shwr] = mcshwr.AncestorStart().Z();
4999  fData->mcshwr_AncestorendX[shwr] = mcshwr.AncestorEnd().X();
5000  fData->mcshwr_AncestorendY[shwr] = mcshwr.AncestorEnd().Y();
5001  fData->mcshwr_AncestorendZ[shwr] = mcshwr.AncestorEnd().Z();
5002  ++shwr;
5003  }
5004  fData->mcshwr_Process.resize(shwr);
5005  fData->mcshwr_MotherProcess.resize(shwr);
5006  fData->mcshwr_AncestorProcess.resize(shwr);
5007  }//End if (fSaveMCShowerInfo){
5008 
5009  //Extract MC Track information and fill the Shower branches
5010  if (fSaveMCTrackInfo){
5011  fData->no_mctracks = nMCTracks;
5012  size_t trk = 0;
5013  for(std::vector<sim::MCTrack>::const_iterator imctrk = mctrackh->begin();imctrk != mctrackh->end(); ++imctrk) {
5014  const sim::MCTrack& mctrk = *imctrk;
5015  TLorentzVector tpcstart, tpcend, tpcmom;
5016  double plen = driftedLength(detProp, mctrk, tpcstart, tpcend, tpcmom);
5017  fData->mctrk_origin[trk] = mctrk.Origin();
5018  fData->mctrk_pdg[trk] = mctrk.PdgCode();
5019  fData->mctrk_TrackId[trk] = mctrk.TrackID();
5020  fData->mctrk_Process[trk] = mctrk.Process();
5021  fData->mctrk_startX[trk] = mctrk.Start().X();
5022  fData->mctrk_startY[trk] = mctrk.Start().Y();
5023  fData->mctrk_startZ[trk] = mctrk.Start().Z();
5024  fData->mctrk_endX[trk] = mctrk.End().X();
5025  fData->mctrk_endY[trk] = mctrk.End().Y();
5026  fData->mctrk_endZ[trk] = mctrk.End().Z();
5027  fData->mctrk_Motherpdg[trk] = mctrk.MotherPdgCode();
5028  fData->mctrk_MotherTrkId[trk] = mctrk.MotherTrackID();
5029  fData->mctrk_MotherProcess[trk] = mctrk.MotherProcess();
5030  fData->mctrk_MotherstartX[trk] = mctrk.MotherStart().X();
5031  fData->mctrk_MotherstartY[trk] = mctrk.MotherStart().Y();
5032  fData->mctrk_MotherstartZ[trk] = mctrk.MotherStart().Z();
5033  fData->mctrk_MotherendX[trk] = mctrk.MotherEnd().X();
5034  fData->mctrk_MotherendY[trk] = mctrk.MotherEnd().Y();
5035  fData->mctrk_MotherendZ[trk] = mctrk.MotherEnd().Z();
5036  fData->mctrk_Ancestorpdg[trk] = mctrk.AncestorPdgCode();
5037  fData->mctrk_AncestorTrkId[trk] = mctrk.AncestorTrackID();
5038  fData->mctrk_AncestorProcess[trk] = mctrk.AncestorProcess();
5039  fData->mctrk_AncestorstartX[trk] = mctrk.AncestorStart().X();
5040  fData->mctrk_AncestorstartY[trk] = mctrk.AncestorStart().Y();
5041  fData->mctrk_AncestorstartZ[trk] = mctrk.AncestorStart().Z();
5042  fData->mctrk_AncestorendX[trk] = mctrk.AncestorEnd().X();
5043  fData->mctrk_AncestorendY[trk] = mctrk.AncestorEnd().Y();
5044  fData->mctrk_AncestorendZ[trk] = mctrk.AncestorEnd().Z();
5045 
5046  fData->mctrk_len_drifted[trk] = plen;
5047 
5048  if (plen != 0){
5049  fData->mctrk_startX_drifted[trk] = tpcstart.X();
5050  fData->mctrk_startY_drifted[trk] = tpcstart.Y();
5051  fData->mctrk_startZ_drifted[trk] = tpcstart.Z();
5052  fData->mctrk_endX_drifted[trk] = tpcend.X();
5053  fData->mctrk_endY_drifted[trk] = tpcend.Y();
5054  fData->mctrk_endZ_drifted[trk] = tpcend.Z();
5055  fData->mctrk_p_drifted[trk] = tpcmom.Vect().Mag();
5056  fData->mctrk_px_drifted[trk] = tpcmom.X();
5057  fData->mctrk_py_drifted[trk] = tpcmom.Y();
5058  fData->mctrk_pz_drifted[trk] = tpcmom.Z();
5059  }
5060  ++trk;
5061  }
5062 
5063  fData->mctrk_Process.resize(trk);
5064  fData->mctrk_MotherProcess.resize(trk);
5065  fData->mctrk_AncestorProcess.resize(trk);
5066  }//End if (fSaveMCTrackInfo){
5067 
5068 
5069  //GEANT particles information
5070  if (fSaveGeantInfo){
5071 
5072  const sim::ParticleList& plist = pi_serv->ParticleList();
5073 
5074  std::string pri("primary");
5075  int primary=0;
5076  int active = 0;
5077  size_t geant_particle=0;
5078  sim::ParticleList::const_iterator itPart = plist.begin(),
5079  pend = plist.end(); // iterator to pairs (track id, particle)
5080 
5081  // helper map track ID => index
5082  std::map<int, size_t> TrackIDtoIndex;
5083  std::vector<int> gpdg;
5084  std::vector<int> gmother;
5085  for(size_t iPart = 0; (iPart < plist.size()) && (itPart != pend); ++iPart){
5086  const simb::MCParticle* pPart = (itPart++)->second;
5087  if (!pPart) {
5089  << "GEANT particle #" << iPart << " returned a null pointer";
5090  }
5091 
5092  //++geant_particle;
5093  bool isPrimary = pPart->Process() == pri;
5094  int TrackID = pPart->TrackId();
5095  TrackIDtoIndex.emplace(TrackID, iPart);
5096  gpdg.push_back(pPart->PdgCode());
5097  gmother.push_back(pPart->Mother());
5098  if (iPart < fData->GetMaxGEANTparticles()) {
5099  if (pPart->E()<fG4minE&&(!isPrimary)) continue;
5100  if (isPrimary) ++primary;
5101 
5102  TLorentzVector mcstart, mcend, mcstartdrifted, mcenddrifted;
5103  unsigned int pstarti, pendi, pstartdriftedi, penddriftedi; //mcparticle indices for starts and ends in tpc or drifted volumes
5104  double plen = length(*pPart, mcstart, mcend, pstarti, pendi);
5105  double plendrifted = driftedLength(detProp, *pPart, mcstartdrifted, mcenddrifted, pstartdriftedi, penddriftedi);
5106 
5107  bool isActive = plen != 0;
5108  bool isDrifted = plendrifted!= 0;
5109  if (plen) ++active;
5110 
5111  fData->process_primary[geant_particle] = int(isPrimary);
5112  fData->processname[geant_particle]= pPart->Process();
5113  fData->Mother[geant_particle]=pPart->Mother();
5114  fData->TrackId[geant_particle]=TrackID;
5115  fData->pdg[geant_particle]=pPart->PdgCode();
5116  fData->status[geant_particle] = pPart->StatusCode();
5117  fData->Eng[geant_particle]=pPart->E();
5118  fData->EndE[geant_particle]=pPart->EndE();
5119  fData->Mass[geant_particle]=pPart->Mass();
5120  fData->Px[geant_particle]=pPart->Px();
5121  fData->Py[geant_particle]=pPart->Py();
5122  fData->Pz[geant_particle]=pPart->Pz();
5123  fData->P[geant_particle]=pPart->Momentum().Vect().Mag();
5124  fData->StartPointx[geant_particle]=pPart->Vx();
5125  fData->StartPointy[geant_particle]=pPart->Vy();
5126  fData->StartPointz[geant_particle]=pPart->Vz();
5127  fData->StartT[geant_particle] = pPart->T();
5128  fData->EndPointx[geant_particle]=pPart->EndPosition()[0];
5129  fData->EndPointy[geant_particle]=pPart->EndPosition()[1];
5130  fData->EndPointz[geant_particle]=pPart->EndPosition()[2];
5131  fData->EndT[geant_particle] = pPart->EndT();
5132  fData->theta[geant_particle] = pPart->Momentum().Theta();
5133  fData->phi[geant_particle] = pPart->Momentum().Phi();
5134  fData->theta_xz[geant_particle] = std::atan2(pPart->Px(), pPart->Pz());
5135  fData->theta_yz[geant_particle] = std::atan2(pPart->Py(), pPart->Pz());
5136  fData->pathlen[geant_particle] = plen;
5137  fData->pathlen_drifted[geant_particle] = plendrifted;
5138  fData->NumberDaughters[geant_particle]=pPart->NumberDaughters();
5139  fData->inTPCActive[geant_particle] = int(isActive);
5140  fData->inTPCDrifted[geant_particle] = int(isDrifted);
5141  art::Ptr<simb::MCTruth> const& mc_truth = pi_serv->ParticleToMCTruth_P(pPart);
5142  if (mc_truth){
5143  fData->origin[geant_particle] = mc_truth->Origin();
5144  fData->MCTruthIndex[geant_particle] = mc_truth.key();
5145  }
5146  if (isActive){
5147  fData->StartPointx_tpcAV[geant_particle] = mcstart.X();
5148  fData->StartPointy_tpcAV[geant_particle] = mcstart.Y();
5149  fData->StartPointz_tpcAV[geant_particle] = mcstart.Z();
5150  fData->StartT_tpcAV[geant_particle] = mcstart.T();
5151  fData->StartE_tpcAV[geant_particle] = pPart->E(pstarti);
5152  fData->StartP_tpcAV[geant_particle] = pPart->P(pstarti);
5153  fData->StartPx_tpcAV[geant_particle] = pPart->Px(pstarti);
5154  fData->StartPy_tpcAV[geant_particle] = pPart->Py(pstarti);
5155  fData->StartPz_tpcAV[geant_particle] = pPart->Pz(pstarti);
5156  fData->EndPointx_tpcAV[geant_particle] = mcend.X();
5157  fData->EndPointy_tpcAV[geant_particle] = mcend.Y();
5158  fData->EndPointz_tpcAV[geant_particle] = mcend.Z();
5159  fData->EndT_tpcAV[geant_particle] = mcend.T();
5160  fData->EndE_tpcAV[geant_particle] = pPart->E(pendi);
5161  fData->EndP_tpcAV[geant_particle] = pPart->P(pendi);
5162  fData->EndPx_tpcAV[geant_particle] = pPart->Px(pendi);
5163  fData->EndPy_tpcAV[geant_particle] = pPart->Py(pendi);
5164  fData->EndPz_tpcAV[geant_particle] = pPart->Pz(pendi);
5165  }
5166  if (isDrifted){
5167  fData->StartPointx_drifted[geant_particle] = mcstartdrifted.X();
5168  fData->StartPointy_drifted[geant_particle] = mcstartdrifted.Y();
5169  fData->StartPointz_drifted[geant_particle] = mcstartdrifted.Z();
5170  fData->StartT_drifted[geant_particle] = mcstartdrifted.T();
5171  fData->StartE_drifted[geant_particle] = pPart->E(pstartdriftedi);
5172  fData->StartP_drifted[geant_particle] = pPart->P(pstartdriftedi);
5173  fData->StartPx_drifted[geant_particle] = pPart->Px(pstartdriftedi);
5174  fData->StartPy_drifted[geant_particle] = pPart->Py(pstartdriftedi);
5175  fData->StartPz_drifted[geant_particle] = pPart->Pz(pstartdriftedi);
5176  fData->EndPointx_drifted[geant_particle] = mcenddrifted.X();
5177  fData->EndPointy_drifted[geant_particle] = mcenddrifted.Y();
5178  fData->EndPointz_drifted[geant_particle] = mcenddrifted.Z();
5179  fData->EndT_drifted[geant_particle] = mcenddrifted.T();
5180  fData->EndE_drifted[geant_particle] = pPart->E(penddriftedi);
5181  fData->EndP_drifted[geant_particle] = pPart->P(penddriftedi);
5182  fData->EndPx_drifted[geant_particle] = pPart->Px(penddriftedi);
5183  fData->EndPy_drifted[geant_particle] = pPart->Py(penddriftedi);
5184  fData->EndPz_drifted[geant_particle] = pPart->Pz(penddriftedi);
5185  }
5186 
5187  //access auxiliary detector parameters
5188  if (fSaveAuxDetInfo) {
5189  unsigned short nAD = 0; // number of cells that particle hit
5190 
5191  // find deposit of this particle in each of the detector cells
5192  for (const sim::AuxDetSimChannel* c: fAuxDetSimChannels) {
5193 
5194  // find if this cell has a contribution (IDE) from this particle,
5195  // and which one
5196  const std::vector<sim::AuxDetIDE>& setOfIDEs = c->AuxDetIDEs();
5197  // using a C++ "lambda" function here; this one:
5198  // - sees only TrackID from the current scope
5199  // - takes one parameter: the AuxDetIDE to be tested
5200  // - returns if that IDE belongs to the track we are looking for
5202  = std::find_if(
5203  setOfIDEs.begin(), setOfIDEs.end(),
5204  [TrackID](const sim::AuxDetIDE& IDE){ return IDE.trackID == TrackID; }
5205  );
5206  if (iIDE == setOfIDEs.end()) continue;
5207 
5208  // now iIDE points to the energy released by the track #i (TrackID)
5209 
5210  // look for IDE with matching trackID
5211  // find trackIDs stored in setOfIDEs with the same trackID, but negative,
5212  // this is an untracked particle who's energy should be added as deposited by this original trackID
5213  float totalE = 0.; // total energy deposited around by the GEANT particle in this cell
5214  for(const auto& adtracks: setOfIDEs) {
5215  if( fabs(adtracks.trackID) == TrackID )
5216  totalE += adtracks.energyDeposited;
5217  } // for
5218 
5219  // fill the structure
5220  if (nAD < kMaxAuxDets) {
5221  fData->AuxDetID[geant_particle][nAD] = c->AuxDetID();
5222  fData->entryX[geant_particle][nAD] = iIDE->entryX;
5223  fData->entryY[geant_particle][nAD] = iIDE->entryY;
5224  fData->entryZ[geant_particle][nAD] = iIDE->entryZ;
5225  fData->entryT[geant_particle][nAD] = iIDE->entryT;
5226  fData->exitX[geant_particle][nAD] = iIDE->exitX;
5227  fData->exitY[geant_particle][nAD] = iIDE->exitY;
5228  fData->exitZ[geant_particle][nAD] = iIDE->exitZ;
5229  fData->exitT[geant_particle][nAD] = iIDE->exitT;
5230  fData->exitPx[geant_particle][nAD] = iIDE->exitMomentumX;
5231  fData->exitPy[geant_particle][nAD] = iIDE->exitMomentumY;
5232  fData->exitPz[geant_particle][nAD] = iIDE->exitMomentumZ;
5233  fData->CombinedEnergyDep[geant_particle][nAD] = totalE;
5234  }
5235  ++nAD;
5236  } // for aux det sim channels
5237  fData->NAuxDets[geant_particle] = nAD;
5238 
5239  if (nAD > kMaxAuxDets) {
5240  // got this error? consider increasing kMaxAuxDets
5241  mf::LogError("AnalysisTree:limits")
5242  << "particle #" << iPart
5243  << " touches " << nAD << " auxiliary detector cells, only "
5244  << kMaxAuxDets << " of them are saved in the tree";
5245  } // if too many detector cells
5246  } // if (fSaveAuxDetInfo)
5247 
5248  ++geant_particle;
5249  }
5250  else if (iPart == fData->GetMaxGEANTparticles()) {
5251  // got this error? it might be a bug,
5252  // since the structure should have enough room for everything
5253  mf::LogError("AnalysisTree:limits") << "event has "
5254  << plist.size() << " MC particles, only "
5255  << fData->GetMaxGEANTparticles() << " will be stored in tree";
5256  }
5257  } // for particles
5258 
5259  fData->geant_list_size_in_tpcAV = active;
5260  fData->no_primaries = primary;
5261  fData->geant_list_size = geant_particle;
5262  fData->processname.resize(geant_particle);
5263  MF_LOG_DEBUG("AnalysisTree")
5264  << "Counted "
5265  << fData->geant_list_size << " GEANT particles ("
5266  << fData->geant_list_size_in_tpcAV << " in AV), "
5267  << fData->no_primaries << " primaries, "
5268  << fData->genie_no_primaries << " GENIE particles";
5269 
5270  FillWith(fData->MergedId, 0);
5271 
5272  // for each particle, consider all the direct ancestors with the same
5273  // PDG ID, and mark them as belonging to the same "group"
5274  // (having the same MergedId)
5275  /* turn off for now
5276  int currentMergedId = 1;
5277  for(size_t iPart = 0; iPart < geant_particle; ++iPart){
5278  // if the particle already belongs to a group, don't bother
5279  if (fData->MergedId[iPart]) continue;
5280  // the particle starts its own group
5281  fData->MergedId[iPart] = currentMergedId;
5282  int currentMotherTrackId = fData->Mother[iPart];
5283  while (currentMotherTrackId > 0) {
5284  if (TrackIDtoIndex.find(currentMotherTrackId)==TrackIDtoIndex.end()) break;
5285  size_t gindex = TrackIDtoIndex[currentMotherTrackId];
5286  if (gindex<0||gindex>=plist.size()) break;
5287  // if the mother particle is of a different type,
5288  // don't bother with iPart ancestry any further
5289  if (gpdg[gindex]!=fData->pdg[iPart]) break;
5290  if (TrackIDtoIndex.find(currentMotherTrackId)!=TrackIDtoIndex.end()){
5291  size_t igeantMother = TrackIDtoIndex[currentMotherTrackId];
5292  if (igeantMother>=0&&igeantMother<geant_particle){
5293  fData->MergedId[igeantMother] = currentMergedId;
5294  }
5295  }
5296  currentMotherTrackId = gmother[gindex];
5297  }
5298  ++currentMergedId;
5299  }// for merging check
5300  */
5301  } // if (fSaveGeantInfo)
5302 
5303  // Now we have the GEANT info, see if we can match the protoDUNE generator particles
5304  if(fSaveProtoInfo){
5305  for(Int_t prt = 0; prt < nProtoPrimaries; ++prt){
5306  for(Int_t gnt = 0; gnt < fData->geant_list_size; ++gnt){
5307 // if(fData->proto_pdg[prt] == fData->pdg[gnt] && fData->proto_px[prt] == fData->Px[gnt]){
5308  if(fData->proto_pdg[prt] == fData->pdg[gnt] && std::fabs(fData->proto_px[prt] - fData->Px[gnt]) < 0.0001){
5309  fData->proto_geantTrackID[prt] = fData->TrackId[gnt];
5310  fData->proto_geantIndex[prt] = gnt;
5311  break;
5312  }
5313  } // End GEANT loop
5314  } // End protoDUNE generator loop
5315  } // End ProtoDUNE generator if statement
5316 
5317  }//if (mcevts_truth)
5318  }//if (isMC){
5319  fData->taulife = detProp.ElectronLifetime();
5320  fTree->Fill();
5321 
5322  if (mf::isDebugEnabled()) {
5323  // use mf::LogDebug instead of MF_LOG_DEBUG because we reuse it in many lines;
5324  // thus, we protect this part of the code with the line above
5325  mf::LogDebug logStream("AnalysisTreeStructure");
5326  logStream
5327  << "Tree data structure contains:"
5328  << "\n - " << fData->no_hits << " hits (" << fData->GetMaxHits() << ")"
5329  << "\n - " << fData->genie_no_primaries << " genie primaries (" << fData->GetMaxGeniePrimaries() << ")"
5330  << "\n - " << fData->geant_list_size << " GEANT particles (" << fData->GetMaxGEANTparticles() << "), "
5331  << fData->no_primaries << " primaries"
5332  << "\n - " << fData->geant_list_size_in_tpcAV << " GEANT particles in AV "
5333  << "\n - " << ((int) fData->kNTracker) << " trackers:"
5334  ;
5335 
5336  size_t iTracker = 0;
5337  for (auto tracker = fData->TrackData.cbegin();
5338  tracker != fData->TrackData.cend(); ++tracker, ++iTracker
5339  ) {
5340  logStream
5341  << "\n -> " << tracker->ntracks << " " << fTrackModuleLabel[iTracker]
5342  << " tracks (" << tracker->GetMaxTracks() << ")"
5343  ;
5344  for (int iTrk = 0; iTrk < tracker->ntracks; ++iTrk) {
5345  logStream << "\n [" << iTrk << "] "<< tracker->ntrkhits[iTrk][0];
5346  for (size_t ipl = 1; ipl < tracker->GetMaxPlanesPerTrack(iTrk); ++ipl)
5347  logStream << " + " << tracker->ntrkhits[iTrk][ipl];
5348  logStream << " hits (" << tracker->GetMaxHitsPerTrack(iTrk, 0);
5349  for (size_t ipl = 1; ipl < tracker->GetMaxPlanesPerTrack(iTrk); ++ipl)
5350  logStream << " + " << tracker->GetMaxHitsPerTrack(iTrk, ipl);
5351  logStream << ")";
5352  } // for tracks
5353  } // for trackers
5354  } // if logging enabled
5355 
5356 
5357  // if we don't use a permanent buffer (which can be huge),
5358  // delete the current buffer, and we'll create a new one on the next event
5359  if (!fUseBuffer) {
5360  MF_LOG_DEBUG("AnalysisTreeStructure") << "Freeing the tree data structure";
5361  DestroyData();
5362  }
5363 } // dune::AnalysisTree::analyze()
double E(const int i=0) const
Definition: MCParticle.h:233
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
code to link reconstructed objects back to the MC truth information
void DestroyData()
Destroy the local buffers (existing branches will point to invalid address!)
void XYZ(double *xyz) const
Legacy method to access vertex position, preserved to avoid breaking code. Please try to use Vertex::...
Definition: Vertex.cxx:36
bool fSaveRawDigitInfo
whether to extract and save Hit information
void CreateTree(bool bClearData=false)
Create the output tree and the data structures, if needed.
simb::Origin_t Origin() const
Definition: MCTrack.h:40
int PdgCode() const
Definition: MCParticle.h:212
constexpr std::uint32_t timeLow() const
Definition: Timestamp.h:29
double VertexMomentum() const
Definition: Track.h:142
double driftedLength(detinfo::DetectorPropertiesData const &detProp, const simb::MCParticle &part, TLorentzVector &start, TLorentzVector &end, unsigned int &starti, unsigned int &endi)
const MCStep & End() const
Definition: MCShower.h:56
double Py(const int i=0) const
Definition: MCParticle.h:231
float SummedADCaverage() const
Returns the average signal ADC counts of the cluster hits.
Definition: Cluster.h:656
float IntegralAverage() const
Returns the average charge of the cluster hits.
Definition: Cluster.h:622
bool fSavePFParticleInfo
whether to extract and save Shower information
Energy deposited on a readout channel by simulated tracks.
Definition: SimChannel.h:140
double EndE() const
Definition: MCParticle.h:244
constexpr int kMaxNPFPNeutrinos
const TLorentzVector & EndPosition() const
Definition: MCParticle.h:225
std::vector< std::string > fCosmicTaggerAssocLabel
whether to extract and save CNN information
std::map< art::Ptr< recob::PFParticle >, ClusterVector > PFParticlesToClusters
constexpr unsigned short kMaxAuxDets
max number of auxiliary detector cells per MC particle
unsigned int TrackID() const
Definition: MCShower.h:53
static bool IsNeutrino(const art::Ptr< recob::PFParticle > particle)
Determine whether a particle has been reconstructed as a neutrino.
bool fSaveAuxDetInfo
whether to extract and save auxiliary detector data
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
std::string string
Definition: nybbler.cc:12
unsigned int ID
void SetVertexAddresses(size_t iVertexAlg)
const simb::MCParticle * TrackIdToParticle_P(int id) const
const std::string & AncestorProcess() const
Definition: MCTrack.h:57
int Mother() const
Definition: MCParticle.h:213
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:34
simb::Origin_t Origin() const
Definition: MCTruth.h:74
double Mass() const
Definition: MCParticle.h:239
double Px(const int i=0) const
Definition: MCParticle.h:230
bool fSaveCaloCosmics
save calorimetry information for cosmics
const MCStep & MotherEnd() const
Definition: MCTrack.h:53
unsigned int AncestorTrackID() const
Definition: MCTrack.h:56
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:102
float StartWire() const
Returns the wire coordinate of the start of the cluster.
Definition: Cluster.h:286
std::string fSpacePointSolverModuleLabel
int PdgCode() const
Definition: MCShower.h:52
Vector_t VertexDirection() const
Definition: Track.h:132
static bool IsShower(const art::Ptr< recob::PFParticle > particle)
Determine whether a particle has been reconstructed as shower-like.
std::vector< std::string > fMVAPIDShowerModuleLabel
bool get(SelectorBase const &, Handle< PROD > &result) const
Definition: DataViewImpl.h:606
std::vector< const sim::IDE * > TrackIdToSimIDEs_Ps(int const &id) const
int StatusCode() const
Definition: MCParticle.h:211
std::pair< float, std::string > P
int PdgCode() const
Return the type of particle as a PDG ID.
Definition: PFParticle.h:83
Set of hits with a 2D structure.
Definition: Cluster.h:71
std::vector< art::Ptr< recob::Shower > > ShowerVector
intermediate_table::const_iterator const_iterator
float EndTick() const
Returns the tick coordinate of the end of the cluster.
Definition: Cluster.h:342
int NParticles() const
Definition: MCTruth.h:75
float MultipleHitDensity() const
Density of wires in the cluster with more than one hit.
Definition: Cluster.h:723
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Cluster finding and building.
std::string Process() const
Definition: MCParticle.h:215
std::vector< std::string > fContainmentTaggerAssocLabel
int AncestorPdgCode() const
Definition: MCTrack.h:55
constexpr int kMaxExternCounts
bool isValid() const
Returns if the cluster is valid (that is, if its ID is not invalid)
Definition: Cluster.h:753
const MCStep & End() const
Definition: MCTrack.h:45
float StartAngle() const
Returns the starting angle of the cluster.
Definition: Cluster.h:475
std::vector< std::string > fMCT0FinderLabel
std::map< art::Ptr< recob::PFParticle >, TrackVector > PFParticlesToTracks
const art::Ptr< simb::MCTruth > & ParticleToMCTruth_P(const simb::MCParticle *p) const
std::vector< std::string > fVertexModuleLabel
int NumberDaughters() const
Definition: MCParticle.h:217
unsigned int MotherTrackID() const
Definition: MCTrack.h:50
std::map< art::Ptr< recob::PFParticle >, VertexVector > PFParticlesToVertices
int TrackId() const
Definition: MCParticle.h:210
bool isValid() const noexcept
Definition: Handle.h:191
std::vector< std::string > fFlashT0FinderLabel
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
const TVector3 & StartDir() const
Definition: MCShower.h:82
Collection of particles crossing one auxiliary detector cell.
bool isRealData() const
bool fSaveExternCounterInfo
whether to extract and save Flash information
float EndCharge() const
Returns the charge on the last wire of the cluster.
Definition: Cluster.h:498
double dEdx(float dqdx, float Efield)
Definition: doAna.cpp:21
float SummedADC() const
Returns the total charge of the cluster from signal ADC counts.
Definition: Cluster.h:634
double bdist(const TVector3 &pos)
std::string fCosmicClusterTaggerAssocLabel
simb::Origin_t Origin() const
Definition: MCShower.h:50
const double e
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:167
std::map< art::Ptr< recob::PFParticle >, ShowerVector > PFParticlesToShowers
std::vector< std::string > fMVAPIDTrackModuleLabel
double Px() const
Definition: MCStep.h:46
Timestamp time() const
AnalysisTreeDataStruct::SubRunData_t SubRunData
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
int MotherPdgCode() const
Definition: MCShower.h:58
const std::string & AncestorProcess() const
Definition: MCShower.h:66
std::vector< art::Ptr< recob::PFParticle > > PFParticleVector
bool fSaveVertexInfo
whether to extract and save Track information
const art::Ptr< simb::MCTruth > & TrackIdToMCTruth_P(int id) const
std::void_t< T > n
static void CollectShowers(const art::Event &evt, const std::string &label, ShowerVector &showerVector, PFParticlesToShowers &particlesToShowers)
Collect the reconstructed PFParticles and associated Showers from the ART event record.
std::vector< Handle< PROD > > getMany(SelectorBase const &selector=MatchAllSelector{}) const
Definition: DataViewImpl.h:479
double P(const int i=0) const
Definition: MCParticle.h:234
key_type key() const noexcept
Definition: Ptr.h:216
float Width() const
A measure of the cluster width, in homogenized units.
Definition: Cluster.h:727
size_t GetNShowerAlgos() const
Returns the number of shower algorithms configured.
Point_t const & Vertex() const
Definition: Track.h:124
static void CollectVertices(const art::Event &evt, const std::string &label, VertexVector &vertexVector, PFParticlesToVertices &particlesToVertices)
Collect the reconstructed PFParticles and associated Vertices from the ART event record.
std::string fPandoraNuVertexModuleLabel
double T(const int i=0) const
Definition: MCParticle.h:224
double Z() const
Definition: MCStep.h:44
std::vector< std::string > fCalorimetryModuleLabel
constexpr int kMaxClusters
struct recob::OpFlashPtrSortByPE_t OpFlashPtrSortByPE
bool fUseBuffer
whether to use a permanent buffer (faster, huge memory)
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:78
bool IsPrimary() const
Returns whether the particle is the root of the flow.
Definition: PFParticle.h:86
bool isDebugEnabled()
double Py() const
Definition: MCStep.h:47
double Y() const
Definition: MCStep.h:43
std::vector< art::Ptr< recob::Track > > TrackVector
const MCStep & AncestorStart() const
Definition: MCShower.h:67
RunNumber_t run() const
Definition: DataViewImpl.cc:71
static int max(int a, int b)
double EndT() const
Definition: MCParticle.h:229
size_t GetNVertexAlgos() const
constexpr int kNplanes
std::string fPFParticleModuleLabel
const MCStep & AncestorStart() const
Definition: MCTrack.h:58
static void CollectTracks(const art::Event &evt, const std::string &label, TrackVector &trackVector, PFParticlesToTracks &particlesToTracks)
Collect the reconstructed PFParticles and associated Tracks from the ART event record.
bool fSaveProtoInfo
whether to extract and save Genie information
const std::string & MotherProcess() const
Definition: MCShower.h:60
static void CollectPFParticles(const art::Event &evt, const std::string &label, PFParticleVector &particleVector)
Collect the reconstructed PFParticles from the ART event record.
double dEdx() const
Definition: MCShower.h:81
bool fSaveTrackInfo
whether to extract and save Raw Digit information
bool fSaveClusterInfo
whether to extract and save Vertex information
unsigned int AncestorTrackID() const
Definition: MCShower.h:65
const MCStep & AncestorEnd() const
Definition: MCShower.h:68
int PdgCode() const
Definition: MCTrack.h:41
const MCStep & DetProfile() const
Definition: MCShower.h:70
Detector simulation of raw signals on wires.
const sim::ParticleList & ParticleList() const
bool isCosmics
if it contains cosmics
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:76
static void SelectNeutrinoPFParticles(const PFParticleVector &inputParticles, PFParticleVector &outputParticles)
Select reconstructed neutrino particles from a list of all reconstructed particles.
std::vector< std::string > fTrackModuleLabel
bool fSavePandoraNuVertexInfo
whether to extract and save Cluster information
const MCStep & Start() const
Definition: MCShower.h:55
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
constexpr int kMaxVertices
double Vx(const int i=0) const
Definition: MCParticle.h:221
int ID() const
Definition: Track.h:198
constexpr int kMaxHits
geo::View_t View() const
Returns the view for this cluster.
Definition: Cluster.h:741
const MCStep & MotherStart() const
Definition: MCTrack.h:52
bool fSaveFlashInfo
whether to extract and save nu vertex information from Pandora
constexpr int kMaxFlashes
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
int MotherPdgCode() const
Definition: MCTrack.h:49
const MCStep & MotherEnd() const
Definition: MCShower.h:62
std::string fExternalCounterModuleLabel
float StartCharge() const
Returns the charge on the first wire of the cluster.
Definition: Cluster.h:454
bool fSaveMCShowerInfo
whether to extract and save Geant information
bool fSaveMCTrackInfo
whether to extract and save MC Shower information
ID_t ID() const
Identifier of this cluster.
Definition: Cluster.h:738
Vector_t EndDirection() const
Definition: Track.h:133
MC truth information to make RawDigits and do back tracking.
std::vector< sim::TrackIDE > HitToEveTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:220
const std::string & Process() const
Definition: MCTrack.h:43
const std::string & MotherProcess() const
Definition: MCTrack.h:51
bool fSaveGeantInfo
whether to extract and save ProtDUNE beam simulation information
#define MF_LOG_DEBUG(id)
static bool IsTrack(const art::Ptr< recob::PFParticle > particle)
Determine whether a particle has been reconstructed as track-like.
double Pz(const int i=0) const
Definition: MCParticle.h:232
float fG4minE
Energy threshold to save g4 particle info.
size_t GetNTrackers() const
Returns the number of trackers configured.
void SetTrackerAddresses(size_t iTracker)
double E() const
Definition: MCStep.h:49
unsigned int MotherTrackID() const
Definition: MCShower.h:59
double Vz(const int i=0) const
Definition: MCParticle.h:223
std::vector< double > SimIDEsToXYZ(std::vector< sim::IDE > const &ides) const
const std::string & Process() const
Definition: MCShower.h:54
double Pz() const
Definition: MCStep.h:48
std::vector< art::Ptr< recob::Vertex > > VertexVector
std::unique_ptr< AnalysisTreeDataStruct > fData
bool bIgnoreMissingShowers
whether to ignore missing shower information
EventNumber_t event() const
Definition: EventID.h:116
Point_t const & End() const
Definition: Track.h:125
const MCStep & MotherStart() const
Definition: MCShower.h:61
detail::Node< FrameID, bool > PlaneID
Definition: CRTID.h:125
bool fSaveCnnInfo
whether to extract and save SpacePointSolver information
const MCStep & Start() const
Definition: MCTrack.h:44
TDCIDEs_t const & TDCIDEMap() const
Returns all the deposited energy information as stored.
Definition: SimChannel.h:328
double X() const
Definition: MCStep.h:42
std::vector< const sim::IDE * > HitToSimIDEs_Ps(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
std::size_t getView(std::string const &moduleLabel, std::string const &productInstanceName, std::string const &processName, std::vector< ELEMENT const * > &result) const
Definition: DataViewImpl.h:500
std::vector< art::Ptr< recob::Cluster > > ClusterVector
unsigned int NHits() const
Number of hits in the cluster.
Definition: Cluster.h:275
unsigned int TrackID() const
Definition: MCTrack.h:42
bool NeutrinoSet() const
Definition: MCTruth.h:78
TCEvent evt
Definition: DataStructs.cxx:7
int AncestorPdgCode() const
Definition: MCShower.h:64
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
std::map< std::string, double > mvaOutput
Definition: MVAPIDResult.h:27
void CreateData(bool bClearData=false)
Creates the structure for the tree data; optionally initializes it.
float EndAngle() const
Returns the ending angle of the cluster.
Definition: Cluster.h:519
void Uncompress(const std::vector< short > &adc, std::vector< short > &uncompressed, raw::Compress_t compress)
Uncompresses a raw data buffer.
Definition: raw.cxx:776
double length(const recob::Track &track)
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:85
recob::tracking::Plane Plane
Definition: TrackState.h:17
std::vector< std::string > fFlashMatchAssocLabel
float StartTick() const
Returns the tick coordinate of the start of the cluster.
Definition: Cluster.h:297
constexpr int kMaxTruth
void HitsPurity(detinfo::DetectorClocksData const &clockData, std::vector< art::Ptr< recob::Hit > > const &hits, Int_t &trackid, Float_t &purity, double &maxe)
std::vector< std::string > fShowerModuleLabel
void FillShowers(AnalysisTreeDataStruct::ShowerDataStruct &showerData, std::vector< recob::Shower > const &showers, const bool fSavePFParticleInfo, const std::map< Short_t, Short_t > &showerIDtoPFParticleIDMap) const
Stores the information of all showers into showerData.
bool fSaveHitInfo
whether to extract and save MC Track information
static std::unique_ptr< MVAReader > create(const art::Event &evt, const art::InputTag &tag)
Definition: MVAReader.h:110
EventID id() const
Definition: Event.cc:34
double Vy(const int i=0) const
Definition: MCParticle.h:222
const MCStep & AncestorEnd() const
Definition: MCTrack.h:59
bool fSaveGenieInfo
whether to extract and save CRY particle data
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
Definition: Track.h:49
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Cosmic rays.
Definition: MCTruth.h:24
QTextStream & endl(QTextStream &s)
bool fSaveSpacePointSolverInfo
whether to extract and save PFParticle information
float Integral() const
Returns the total charge of the cluster from hit shape.
Definition: Cluster.h:600
float EndWire() const
Returns the wire coordinate of the end of the cluster.
Definition: Cluster.h:329
vertex reconstruction
bool fSaveShowerInfo
whether to extract and save External Counter information
double AnalysisTree::bdist ( const TVector3 &  pos)
private

Definition at line 5491 of file AnalysisTree_module.cc.

5492 {
5493  double d1 = -ActiveBounds[0] + pos.X();
5494  double d2 = ActiveBounds[1] - pos.X();
5495  double d3 = -ActiveBounds[2] + pos.Y();
5496  double d4 = ActiveBounds[3] - pos.Y();
5497  double d5 = -ActiveBounds[4] + pos.Z();
5498  double d6 = ActiveBounds[5] - pos.Z();
5499 
5500  double Xmin = std::min(d1, d2);
5501  double result = 0;
5502  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 )
5503  else result = std::min(std::min(std::min(std::min(Xmin, d3), d4), d5), d6);
5504  if (result<-1) result = -1;
5505  /*
5506  std::cout << "\n" << std::endl;
5507  std::cout << "-"<<ActiveBounds[0]<<" + "<<pos.X()<< " = " << d1 << "\n"
5508  << ActiveBounds[1]<<" - "<<pos.X()<< " = " << d2 << "\n"
5509  << "-"<<ActiveBounds[2]<<" + "<<pos.Y()<< " = " << d3 << "\n"
5510  << ActiveBounds[3]<<" - "<<pos.Y()<< " = " << d4 << "\n"
5511  << "-"<<ActiveBounds[4]<<" + "<<pos.Z()<< " = " << d5 << "\n"
5512  << ActiveBounds[5]<<" - "<<pos.Z()<< " = " << d6 << "\n"
5513  << "And the Minimum is " << result << std::endl;
5514  */
5515  return result;
5516 }
static QCString result
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
void AnalysisTree::beginSubRun ( const art::SubRun sr)

Definition at line 3490 of file AnalysisTree_module.cc.

3491 {
3492 
3493 // auto potListHandle = sr.getHandle< sumdata::POTSummary >(fPOTModuleLabel);
3494 // if (potListHandle)
3495 // SubRunData.pot=potListHandle->totpot;
3496 // else
3497 // SubRunData.pot=0.;
3498 
3499 }
void dune::AnalysisTree::CheckData ( std::string  caller) const
inlineprivate

Helper function: throws if no data structure is available.

Definition at line 1418 of file AnalysisTree_module.cc.

1419  {
1420  if (fData) return;
1422  << "AnalysisTree::" << caller << ": no data";
1423  } // CheckData()
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
std::unique_ptr< AnalysisTreeDataStruct > fData
void dune::AnalysisTree::CheckTree ( std::string  caller) const
inlineprivate

Helper function: throws if no tree is available.

Definition at line 1425 of file AnalysisTree_module.cc.

1426  {
1427  if (fTree) return;
1429  << "AnalysisTree::" << caller << ": no tree";
1430  } // CheckData()
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
void dune::AnalysisTree::CreateData ( bool  bClearData = false)
inlineprivate

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

Definition at line 1321 of file AnalysisTree_module.cc.

1322  {
1323  if (!fData) {
1324  fData.reset
1325  (new AnalysisTreeDataStruct(GetNTrackers(), GetNVertexAlgos(), GetShowerAlgos()));
1345  }
1346  else {
1347  fData->SetTrackers(GetNTrackers());
1348  fData->SetVertexAlgos(GetNVertexAlgos());
1349  fData->SetShowerAlgos(GetShowerAlgos());
1350 
1351  if (bClearData) fData->Clear();
1352  }
1353  } // CreateData()
bool fSaveRawDigitInfo
whether to extract and save Hit information
bool fSavePFParticleInfo
whether to extract and save Shower information
bool fSaveAuxDetInfo
whether to extract and save auxiliary detector data
bool fSaveExternCounterInfo
whether to extract and save Flash information
bool fSaveVertexInfo
whether to extract and save Track information
size_t GetNVertexAlgos() const
bool fSaveProtoInfo
whether to extract and save Genie information
bool fSaveTrackInfo
whether to extract and save Raw Digit information
bool fSaveClusterInfo
whether to extract and save Vertex information
bool fSavePandoraNuVertexInfo
whether to extract and save Cluster information
bool fSaveFlashInfo
whether to extract and save nu vertex information from Pandora
bool fSaveMCShowerInfo
whether to extract and save Geant information
bool fSaveMCTrackInfo
whether to extract and save MC Shower information
bool fSaveGeantInfo
whether to extract and save ProtDUNE beam simulation information
size_t GetNTrackers() const
Returns the number of trackers configured.
std::unique_ptr< AnalysisTreeDataStruct > fData
bool fSaveCnnInfo
whether to extract and save SpacePointSolver information
std::vector< std::string > GetShowerAlgos() const
Returns the name of configured shower algorithms (converted to string)
bool fSaveHitInfo
whether to extract and save MC Track information
bool fSaveGenieInfo
whether to extract and save CRY particle data
bool fSaveSpacePointSolverInfo
whether to extract and save PFParticle information
bool fSaveShowerInfo
whether to extract and save External Counter information
void AnalysisTree::CreateTree ( bool  bClearData = false)
private

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

Definition at line 3472 of file AnalysisTree_module.cc.

3472  {
3473  if (!fTree) {
3475  fTree = tfs->make<TTree>("anatree","analysis tree");
3476  }
3477  if (!fPOT) {
3479  fPOT = tfs->make<TTree>("pottree","pot tree");
3480  fPOT->Branch("pot",&SubRunData.pot,"pot/D");
3481  fPOT->Branch("potbnbETOR860",&SubRunData.potbnbETOR860,"potbnbETOR860/D");
3482  fPOT->Branch("potbnbETOR875",&SubRunData.potbnbETOR875,"potbnbETOR875/D");
3483  fPOT->Branch("potnumiETORTGT",&SubRunData.potnumiETORTGT,"potnumiETORTGT/D");
3484  }
3485  CreateData(bClearData);
3486  SetAddresses();
3487 } // dune::AnalysisTree::CreateTree()
AnalysisTreeDataStruct::SubRunData_t SubRunData
void SetAddresses()
Sets the addresses of all the tree branches, creating the missing ones.
void CreateData(bool bClearData=false)
Creates the structure for the tree data; optionally initializes it.
void dune::AnalysisTree::DestroyData ( )
inlineprivate

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

Definition at line 1415 of file AnalysisTree_module.cc.

1415 { fData.reset(); }
std::unique_ptr< AnalysisTreeDataStruct > fData
double AnalysisTree::driftedLength ( detinfo::DetectorPropertiesData const &  detProp,
const simb::MCParticle part,
TLorentzVector &  start,
TLorentzVector &  end,
unsigned int &  starti,
unsigned int &  endi 
)
private

Definition at line 5576 of file AnalysisTree_module.cc.

5578 {
5579  auto const* geom = lar::providerFrom<geo::Geometry>();
5580 
5581  //compute the drift x range
5582  double vDrift = detProp.DriftVelocity()*1e-3; //cm/ns
5583  double xrange[2] = {DBL_MAX, -DBL_MAX };
5584  for (unsigned int c=0; c<geom->Ncryostats(); ++c) {
5585  for (unsigned int t=0; t<geom->NTPC(c); ++t) {
5586  double Xat0 = detProp.ConvertTicksToX(0,0,t,c);
5587  double XatT = detProp.ConvertTicksToX(detProp.NumberTimeSamples(),0,t,c);
5588  xrange[0] = std::min(std::min(Xat0, XatT), xrange[0]);
5589  xrange[1] = std::max(std::max(Xat0, XatT), xrange[1]);
5590  }
5591  }
5592 
5593  double result = 0.;
5594  TVector3 disp;
5595  bool first = true;
5596 
5597  for(unsigned int i = 0; i < p.NumberTrajectoryPoints(); ++i) {
5598  // check if the particle is inside a TPC
5599  if (p.Vx(i) >= ActiveBounds[0] && p.Vx(i) <= ActiveBounds[1] &&
5600  p.Vy(i) >= ActiveBounds[2] && p.Vy(i) <= ActiveBounds[3] &&
5601  p.Vz(i) >= ActiveBounds[4] && p.Vz(i) <= ActiveBounds[5]){
5602  // Doing some manual shifting to account for
5603  // an interaction not occuring with the beam dump
5604  // we will reconstruct an x distance different from
5605  // where the particle actually passed to to the time
5606  // being different from in-spill interactions
5607  double newX = p.Vx(i)+(p.T(i)*vDrift);
5608  if (newX < xrange[0] || newX > xrange[1]) continue;
5609  TLorentzVector pos(newX,p.Vy(i),p.Vz(i),p.T());
5610  if(first){
5611  start = pos;
5612  starti=i;
5613  first = false;
5614  }
5615  else {
5616  disp -= pos.Vect();
5617  result += disp.Mag();
5618  }
5619  disp = pos.Vect();
5620  end = pos;
5621  endi = i;
5622  }
5623  }
5624  return result;
5625 }
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
static QCString result
const double e
p
Definition: test.py:223
static int max(int a, int b)
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
double AnalysisTree::driftedLength ( detinfo::DetectorPropertiesData const &  detProp,
const sim::MCTrack mctrack,
TLorentzVector &  tpcstart,
TLorentzVector &  tpcend,
TLorentzVector &  tpcmom 
)
private

Definition at line 5525 of file AnalysisTree_module.cc.

5526  {
5527  auto const* geom = lar::providerFrom<geo::Geometry>();
5528 
5529  //compute the drift x range
5530  double vDrift = detProp.DriftVelocity()*1e-3; //cm/ns
5531  double xrange[2] = {DBL_MAX, -DBL_MAX };
5532  for (unsigned int c=0; c<geom->Ncryostats(); ++c) {
5533  for (unsigned int t=0; t<geom->NTPC(c); ++t) {
5534  double Xat0 = detProp.ConvertTicksToX(0,0,t,c);
5535  double XatT = detProp.ConvertTicksToX(detProp.NumberTimeSamples(),0,t,c);
5536  xrange[0] = std::min(std::min(Xat0, XatT), xrange[0]);
5537  xrange[1] = std::max(std::max(Xat0, XatT), xrange[1]);
5538  }
5539  }
5540 
5541  double result = 0.;
5542  TVector3 disp;
5543  bool first = true;
5544 
5545  for(auto step: mctrack) {
5546  // check if the particle is inside a TPC
5547  if (step.X() >= ActiveBounds[0] && step.X() <= ActiveBounds[1] &&
5548  step.Y() >= ActiveBounds[2] && step.Y() <= ActiveBounds[3] &&
5549  step.Z() >= ActiveBounds[4] && step.Z() <= ActiveBounds[5] ){
5550  // Doing some manual shifting to account for
5551  // an interaction not occuring with the beam dump
5552  // we will reconstruct an x distance different from
5553  // where the particle actually passed to to the time
5554  // being different from in-spill interactions
5555  double newX = step.X()+(step.T()*vDrift);
5556  if (newX < xrange[0] || newX > xrange[1]) continue;
5557 
5558  TLorentzVector pos(newX,step.Y(),step.Z(),step.T());
5559  if(first){
5560  tpcstart = pos;
5561  tpcmom = step.Momentum();
5562  first = false;
5563  }
5564  else {
5565  disp -= pos.Vect();
5566  result += disp.Mag();
5567  }
5568  disp = pos.Vect();
5569  tpcend = pos;
5570  }
5571  }
5572  return result;
5573 }
static QCString result
const double e
static int max(int a, int b)
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
void AnalysisTree::endSubRun ( const art::SubRun sr)

Definition at line 3501 of file AnalysisTree_module.cc.

3502 {
3503 
3504  auto potListHandle = sr.getHandle< sumdata::POTSummary >(fPOTModuleLabel);
3505  if (potListHandle)
3506  SubRunData.pot=potListHandle->totpot;
3507  else
3508  SubRunData.pot=0.;
3509 
3510  art::InputTag itag1("beamdata","bnbETOR860");
3511  auto potSummaryHandlebnbETOR860 = sr.getHandle<sumdata::POTSummary>(itag1);
3512  if (potSummaryHandlebnbETOR860){
3513  SubRunData.potbnbETOR860 = potSummaryHandlebnbETOR860->totpot;
3514  }
3515  else
3517 
3518  art::InputTag itag2("beamdata","bnbETOR875");
3519  auto potSummaryHandlebnbETOR875 = sr.getHandle<sumdata::POTSummary>(itag2);
3520  if (potSummaryHandlebnbETOR875){
3521  SubRunData.potbnbETOR875 = potSummaryHandlebnbETOR875->totpot;
3522  }
3523  else
3525 
3526  art::InputTag itag3("beamdata","numiETORTGT");
3527  auto potSummaryHandlenumiETORTGT = sr.getHandle<sumdata::POTSummary>(itag3);
3528  if (potSummaryHandlenumiETORTGT){
3529  SubRunData.potnumiETORTGT = potSummaryHandlenumiETORTGT->totpot;
3530  }
3531  else
3533 
3534  if (fPOT) fPOT->Fill();
3535 
3536 }
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
AnalysisTreeDataStruct::SubRunData_t SubRunData
void AnalysisTree::FillShower ( AnalysisTreeDataStruct::ShowerDataStruct showerData,
size_t  iShower,
recob::Shower const &  showers,
const bool  fSavePFParticleInfo,
const std::map< Short_t, Short_t > &  showerIDtoPFParticleIDMap 
) const
private

Stores the information of shower in slot iShower of showerData.

Definition at line 5366 of file AnalysisTree_module.cc.

5369  {
5370 
5371  showerData.showerID[iShower] = shower.ID();
5372  showerData.shwr_bestplane[iShower] = shower.best_plane();
5373  showerData.shwr_length[iShower] = shower.Length();
5374 
5375  TVector3 const& dir_start = shower.Direction();
5376  showerData.shwr_startdcosx[iShower] = dir_start.X();
5377  showerData.shwr_startdcosy[iShower] = dir_start.Y();
5378  showerData.shwr_startdcosz[iShower] = dir_start.Z();
5379 
5380  TVector3 const& pos_start = shower.ShowerStart();
5381  showerData.shwr_startx[iShower] = pos_start.X();
5382  showerData.shwr_starty[iShower] = pos_start.Y();
5383  showerData.shwr_startz[iShower] = pos_start.Z();
5384 
5385  if (fSavePFParticleInfo) {
5386  auto mapIter = showerIDtoPFParticleIDMap.find(shower.ID());
5387  if (mapIter != showerIDtoPFParticleIDMap.end()) {
5388  // This vertex has a corresponding PFParticle.
5389  showerData.shwr_hasPFParticle[iShower] = 1;
5390  showerData.shwr_PFParticleID[iShower] = mapIter->second;
5391  }
5392  else
5393  showerData.shwr_hasPFParticle[iShower] = 0;
5394  }
5395 
5396  if (shower.Energy().size() == kNplanes)
5397  std::copy_n
5398  (shower.Energy().begin(), kNplanes, &showerData.shwr_totEng[iShower][0]);
5399  if (shower.dEdx().size() == kNplanes)
5400  std::copy_n
5401  (shower.dEdx().begin(), kNplanes, &showerData.shwr_dedx[iShower][0]);
5402  if (shower.MIPEnergy().size() == kNplanes)
5403  std::copy_n
5404  (shower.MIPEnergy().begin(), kNplanes, &showerData.shwr_mipEng[iShower][0]);
5405 
5406 } // dune::AnalysisTree::FillShower()
bool fSavePFParticleInfo
whether to extract and save Shower information
constexpr int kNplanes
void AnalysisTree::FillShowers ( AnalysisTreeDataStruct::ShowerDataStruct showerData,
std::vector< recob::Shower > const &  showers,
const bool  fSavePFParticleInfo,
const std::map< Short_t, Short_t > &  showerIDtoPFParticleIDMap 
) const
private

Stores the information of all showers into showerData.

Definition at line 5409 of file AnalysisTree_module.cc.

5412  {
5413 
5414  const size_t NShowers = showers.size();
5415 
5416  //
5417  // prepare the data structures, the tree and the connection between them
5418  //
5419 
5420  // allocate enough space for this number of showers
5421  // (but at least for one of them!)
5422  showerData.SetMaxShowers(std::max(NShowers, (size_t) 1));
5423  showerData.Clear(); // clear all the data
5424 
5425  // now set the tree addresses to the newly allocated memory;
5426  // this creates the tree branches in case they are not there yet
5427  showerData.SetAddresses(fTree);
5428  if (NShowers > showerData.GetMaxShowers()) {
5429  // got this error? it might be a bug,
5430  // since we are supposed to have allocated enough space to fit all showers
5431  mf::LogError("AnalysisTree:limits") << "event has " << NShowers
5432  << " " << showerData.Name() << " showers, only "
5433  << showerData.GetMaxShowers() << " stored in tree";
5434  }
5435 
5436  //
5437  // now set the data
5438  //
5439  // set the record of the number of showers
5440  // (data structures are already properly resized)
5441  showerData.nshowers = (Short_t) NShowers;
5442 
5443  // set all the showers one by one
5444  for (size_t i = 0; i < NShowers; ++i) FillShower(showerData, i, showers[i], fSavePFParticleInfo, showerIDtoPFParticleIDMap);
5445 
5446 } // dune::AnalysisTree::FillShowers()
bool fSavePFParticleInfo
whether to extract and save Shower information
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
static int max(int a, int b)
void FillShower(AnalysisTreeDataStruct::ShowerDataStruct &showerData, size_t iShower, recob::Shower const &showers, const bool fSavePFParticleInfo, const std::map< Short_t, Short_t > &showerIDtoPFParticleIDMap) const
Stores the information of shower in slot iShower of showerData.
size_t dune::AnalysisTree::GetNShowerAlgos ( ) const
inlineprivate

Returns the number of shower algorithms configured.

Definition at line 1314 of file AnalysisTree_module.cc.

1314 { return fShowerModuleLabel.size(); }
std::vector< std::string > fShowerModuleLabel
size_t dune::AnalysisTree::GetNTrackers ( ) const
inlineprivate

Returns the number of trackers configured.

Definition at line 1309 of file AnalysisTree_module.cc.

1309 { return fTrackModuleLabel.size(); }
std::vector< std::string > fTrackModuleLabel
size_t dune::AnalysisTree::GetNVertexAlgos ( ) const
inlineprivate

Definition at line 1311 of file AnalysisTree_module.cc.

1311 { return fVertexModuleLabel.size(); }
std::vector< std::string > fVertexModuleLabel
std::vector<std::string> dune::AnalysisTree::GetShowerAlgos ( ) const
inlineprivate

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

Definition at line 1317 of file AnalysisTree_module.cc.

1318  { return { fShowerModuleLabel.begin(), fShowerModuleLabel.end() }; }
std::vector< std::string > fShowerModuleLabel
void AnalysisTree::HitsPurity ( detinfo::DetectorClocksData const &  clockData,
std::vector< art::Ptr< recob::Hit > > const &  hits,
Int_t &  trackid,
Float_t &  purity,
double &  maxe 
)
private

Definition at line 5450 of file AnalysisTree_module.cc.

5451  {
5452 
5453  trackid = -1;
5454  purity = -1;
5455 
5457 
5458  std::map<int,double> trkide;
5459 
5460  for(size_t h = 0; h < hits.size(); ++h){
5461 
5462  art::Ptr<recob::Hit> hit = hits[h];
5463  std::vector<sim::IDE> ides;
5464  //bt_serv->HitToSimIDEs(hit,ides);
5465  std::vector<sim::TrackIDE> eveIDs = bt_serv->HitToEveTrackIDEs(clockData, hit);
5466 
5467  for(size_t e = 0; e < eveIDs.size(); ++e){
5468  //std::cout<<h<<" "<<e<<" "<<eveIDs[e].trackID<<" "<<eveIDs[e].energy<<" "<<eveIDs[e].energyFrac<<std::endl;
5469  trkide[eveIDs[e].trackID] += eveIDs[e].energy;
5470  }
5471  }
5472 
5473  maxe = -1;
5474  double tote = 0;
5475  for (std::map<int,double>::iterator ii = trkide.begin(); ii!=trkide.end(); ++ii){
5476  tote += ii->second;
5477  if ((ii->second)>maxe){
5478  maxe = ii->second;
5479  trackid = ii->first;
5480  }
5481  }
5482 
5483  //std::cout << "the total energy of this reco track is: " << tote << std::endl;
5484 
5485  if (tote>0){
5486  purity = maxe/tote;
5487  }
5488 }
intermediate_table::iterator iterator
const double e
Detector simulation of raw signals on wires.
std::vector< sim::TrackIDE > HitToEveTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
double AnalysisTree::length ( const recob::Track track)
private

Definition at line 5519 of file AnalysisTree_module.cc.

5520 {
5521  return track.Length();
5522 }
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:167
double AnalysisTree::length ( const simb::MCParticle part,
TLorentzVector &  start,
TLorentzVector &  end,
unsigned int &  starti,
unsigned int &  endi 
)
private

Definition at line 5628 of file AnalysisTree_module.cc.

5629 {
5630  double result = 0.;
5631  TVector3 disp;
5632  bool first = true;
5633 
5634  for(unsigned int i = 0; i < p.NumberTrajectoryPoints(); ++i) {
5635  // check if the particle is inside a TPC
5636  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]){
5637  if(first){
5638  start = p.Position(i);
5639  first = false;
5640  starti = i;
5641  }else{
5642  disp -= p.Position(i).Vect();
5643  result += disp.Mag();
5644  }
5645  disp = p.Position(i).Vect();
5646  end = p.Position(i);
5647  endi = i;
5648  }
5649  }
5650  return result;
5651 }
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
static QCString result
p
Definition: test.py:223
void dune::AnalysisTree::SetAddresses ( )
inlineprivate

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

Definition at line 1356 of file AnalysisTree_module.cc.

1357  {
1358  CheckData(__func__); CheckTree(__func__);
1359  fData->SetAddresses
1361  } // SetAddresses()
std::vector< std::string > fVertexModuleLabel
void CheckData(std::string caller) const
Helper function: throws if no data structure is available.
bool isCosmics
if it contains cosmics
std::vector< std::string > fTrackModuleLabel
std::unique_ptr< AnalysisTreeDataStruct > fData
void CheckTree(std::string caller) const
Helper function: throws if no tree is available.
std::vector< std::string > fShowerModuleLabel
void dune::AnalysisTree::SetPFParticleAddress ( )
inlineprivate

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

Definition at line 1405 of file AnalysisTree_module.cc.

1406  {
1407  CheckData(__func__); CheckTree(__func__);
1408  fData->GetPFParticleData().SetAddresses(fTree);
1409  } // SetPFParticleAddress()
void CheckData(std::string caller) const
Helper function: throws if no data structure is available.
std::unique_ptr< AnalysisTreeDataStruct > fData
void CheckTree(std::string caller) const
Helper function: throws if no tree is available.
void dune::AnalysisTree::SetShowerAddresses ( size_t  iShower)
inlineprivate

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

Definition at line 1392 of file AnalysisTree_module.cc.

1393  {
1394  CheckData(__func__); CheckTree(__func__);
1395  if (iShower >= fData->GetNShowerAlgos()) {
1397  << "AnalysisTree::SetShowerAddresses(): no shower algo #" << iShower
1398  << " (" << fData->GetNShowerAlgos() << " available)";
1399  }
1400  fData->GetShowerData(iShower).SetAddresses(fTree);
1401  } // SetShowerAddresses()
void CheckData(std::string caller) const
Helper function: throws if no data structure is available.
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
std::unique_ptr< AnalysisTreeDataStruct > fData
void CheckTree(std::string caller) const
Helper function: throws if no tree is available.
void dune::AnalysisTree::SetTrackerAddresses ( size_t  iTracker)
inlineprivate

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

Definition at line 1365 of file AnalysisTree_module.cc.

1366  {
1367  CheckData(__func__); CheckTree(__func__);
1368  if (iTracker >= fData->GetNTrackers()) {
1370  << "AnalysisTree::SetTrackerAddresses(): no tracker #" << iTracker
1371  << " (" << fData->GetNTrackers() << " available)";
1372  }
1373  fData->GetTrackerData(iTracker)
1374  .SetAddresses(fTree, fTrackModuleLabel[iTracker], isCosmics);
1375  } // SetTrackerAddresses()
void CheckData(std::string caller) const
Helper function: throws if no data structure is available.
bool isCosmics
if it contains cosmics
std::vector< std::string > fTrackModuleLabel
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
std::unique_ptr< AnalysisTreeDataStruct > fData
void CheckTree(std::string caller) const
Helper function: throws if no tree is available.
void dune::AnalysisTree::SetVertexAddresses ( size_t  iVertexAlg)
inlineprivate

Definition at line 1378 of file AnalysisTree_module.cc.

1379  {
1380  CheckData(__func__); CheckTree(__func__);
1381  if (iVertexAlg >= fData->GetNVertexAlgos()) {
1383  << "AnalysisTree::SetVertexAddresses(): no vertex alg #" << iVertexAlg
1384  << " (" << fData->GetNVertexAlgos() << " available)";
1385  }
1386  fData->GetVertexData(iVertexAlg)
1387  .SetAddresses(fTree, fVertexModuleLabel[iVertexAlg], isCosmics);
1388  } // SetVertexAddresses()
std::vector< std::string > fVertexModuleLabel
void CheckData(std::string caller) const
Helper function: throws if no data structure is available.
bool isCosmics
if it contains cosmics
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
std::unique_ptr< AnalysisTreeDataStruct > fData
void CheckTree(std::string caller) const
Helper function: throws if no tree is available.

Member Data Documentation

double dune::AnalysisTree::ActiveBounds[6]
private

Definition at line 1306 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::bIgnoreMissingShowers
private

whether to ignore missing shower information

Definition at line 1300 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fCalDataModuleLabel
private

Definition at line 1250 of file AnalysisTree_module.cc.

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

Definition at line 1267 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fClusterModuleLabel
private

Definition at line 1255 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fCnnModuleLabel
private

Definition at line 1262 of file AnalysisTree_module.cc.

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

Definition at line 1297 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fCosmicClusterTaggerAssocLabel
private

Definition at line 1274 of file AnalysisTree_module.cc.

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

whether to extract and save CNN information

Definition at line 1296 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fCryGenModuleLabel
private

Definition at line 1252 of file AnalysisTree_module.cc.

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

Definition at line 1243 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fDigitModuleLabel
private

Definition at line 1247 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fExternalCounterModuleLabel
private

Definition at line 1258 of file AnalysisTree_module.cc.

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

Definition at line 1298 of file AnalysisTree_module.cc.

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

Definition at line 1271 of file AnalysisTree_module.cc.

float dune::AnalysisTree::fG4minE
private

Energy threshold to save g4 particle info.

Definition at line 1304 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fG4ModuleLabel
private

Definition at line 1254 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fGenieGenModuleLabel
private

Definition at line 1251 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fHitsModuleLabel
private

Definition at line 1248 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fLArG4ModuleLabel
private

Definition at line 1249 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fMCShowerModuleLabel
private

Definition at line 1259 of file AnalysisTree_module.cc.

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

Definition at line 1272 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fMCTrackModuleLabel
private

Definition at line 1260 of file AnalysisTree_module.cc.

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

Definition at line 1269 of file AnalysisTree_module.cc.

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

Definition at line 1270 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fOpFlashModuleLabel
private

Definition at line 1257 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fPandoraNuVertexModuleLabel
private

Definition at line 1256 of file AnalysisTree_module.cc.

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

Definition at line 1268 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fPFParticleModuleLabel
private

Definition at line 1264 of file AnalysisTree_module.cc.

TTree* dune::AnalysisTree::fPOT
private

Definition at line 1239 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fPOTModuleLabel
private

Definition at line 1273 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fProtoGenModuleLabel
private

Definition at line 1253 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveAuxDetInfo
private

whether to extract and save auxiliary detector data

Definition at line 1276 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveCaloCosmics
private

save calorimetry information for cosmics

Definition at line 1303 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveClusterInfo
private

whether to extract and save Vertex information

Definition at line 1287 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveCnnInfo
private

whether to extract and save SpacePointSolver information

Definition at line 1294 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveCryInfo
private

Definition at line 1277 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveExternCounterInfo
private

whether to extract and save Flash information

Definition at line 1290 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveFlashInfo
private

whether to extract and save nu vertex information from Pandora

Definition at line 1289 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveGeantInfo
private

whether to extract and save ProtDUNE beam simulation information

Definition at line 1280 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveGenieInfo
private

whether to extract and save CRY particle data

Definition at line 1278 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveHitInfo
private

whether to extract and save MC Track information

Definition at line 1283 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveMCShowerInfo
private

whether to extract and save Geant information

Definition at line 1281 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveMCTrackInfo
private

whether to extract and save MC Shower information

Definition at line 1282 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSavePandoraNuVertexInfo
private

whether to extract and save Cluster information

Definition at line 1288 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSavePFParticleInfo
private

whether to extract and save Shower information

Definition at line 1292 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveProtoInfo
private

whether to extract and save Genie information

Definition at line 1279 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveRawDigitInfo
private

whether to extract and save Hit information

Definition at line 1284 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveShowerInfo
private

whether to extract and save External Counter information

Definition at line 1291 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveSpacePointSolverInfo
private

whether to extract and save PFParticle information

Definition at line 1293 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveTrackInfo
private

whether to extract and save Raw Digit information

Definition at line 1285 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fSaveVertexInfo
private

whether to extract and save Track information

Definition at line 1286 of file AnalysisTree_module.cc.

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

Definition at line 1266 of file AnalysisTree_module.cc.

std::string dune::AnalysisTree::fSpacePointSolverModuleLabel
private

Definition at line 1261 of file AnalysisTree_module.cc.

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

Definition at line 1263 of file AnalysisTree_module.cc.

TTree* dune::AnalysisTree::fTree
private

Definition at line 1238 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::fUseBuffer
private

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

Definition at line 1275 of file AnalysisTree_module.cc.

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

Definition at line 1265 of file AnalysisTree_module.cc.

bool dune::AnalysisTree::isCosmics
private

if it contains cosmics

Definition at line 1302 of file AnalysisTree_module.cc.

AnalysisTreeDataStruct::SubRunData_t dune::AnalysisTree::SubRunData
private

Definition at line 1245 of file AnalysisTree_module.cc.


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