LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
MCTruthAndFriendsItr.cxx
Go to the documentation of this file.
1 #include "MCTruthAndFriendsItr.h"
2 
3 #include <iostream>
4 #include <iomanip>
5 
6 #ifndef ART_V1
9 #else
10  #include "art/Framework/Core/FindOneP.h"
11  #include "art/Utilities/InputTag.h"
12 #endif
13 
14 
16  std::vector<std::string> const & labels)
17  : evt(evtIn)
18  , fInputModuleLabels(labels)
19  , indx_itr(indices.begin())
20  , nmctruth(0)
21  , imctruth(-1)
22  , thisMCTruth(0)
23  , thisGTruth(0)
24  , thisMCFlux(0)
25  , thisDk2Nu(0)
26  , thisNuChoice(0)
27  , thisLabel("")
28 {
29 
30  // look for any existing MCTruth info in this event
31  if ( fInputModuleLabels.size()==0) {
32  mclists = evt.getMany<std::vector<simb::MCTruth>>();
33  // std::cout << "evt.getManyByType" << std::endl;
34  } else {
35  mclists.resize(fInputModuleLabels.size());
36  for (size_t i=0; i<fInputModuleLabels.size(); ++i) {
38  //std::cout << "evt.getByLabel " << fInputModuleLabels[i] << std::endl;
39  }
40  }
41 
42  //std::cerr << "evg::GENIEDumper::analyze got stuff ---------------- "
43  // << evt.id() << std::endl;
44 
45  for (size_t mcl = 0; mcl < mclists.size(); ++mcl) {
47  if ( ! mclistHandle.isValid() ) {
48  // std::cerr << "not valid mcl=" << mcl << "---------------- " << std::endl;
49  continue;
50  }
51  // processName, moduleLabel, instance (productInstanceName?)
52  /*
53  std::cout << " mclistHandle processName '"
54  << mclistHandle.provenance()->processName() // top of fcl file
55  << "' moduleLabel '"
56  << mclistHandle.provenance()->moduleLabel()
57  << "' productInstanceName '"
58  << mclistHandle.provenance()->productInstanceName()
59  << "'"
60  << std::endl;
61  */
62 
63  std::string handleLabel = mclistHandle.provenance()->moduleLabel();
64  outlabels.push_back(handleLabel);
65  /*
66  std::cerr << "mcl=" << mcl << " '" << handleLabel << "' ---------------- " << std::endl;
67  */
68 
69  // loop over mctruths in a list
70  for(size_t nmc = 0; nmc < mclistHandle->size(); ++nmc) {
71  art::Ptr<simb::MCTruth> mct(mclistHandle, nmc);
72 
73  std::pair<int,int> ipair(mcl,nmc);
74  indices.insert(ipair);
75 
82  }
83 
84  }
85 
86  indx_itr = indices.begin();
87  nmctruth = (int)(indices.size());
88  //std::cout << ".... found nmctruth " << nmctruth
89  // << " nlabels " << outlabels.size()
90  // << std::endl;
91 
92 }
93 
95 {
96  ++imctruth; // started at -1, so first call to Next() prepares us for indx=0
97  thisMCTruth = 0;
98  thisGTruth = 0;
99  thisMCFlux = 0;
100  thisDk2Nu = 0;
101  thisNuChoice = 0;
102  //std::cerr << "Next() called ... moved to imctruth " << imctruth << std::endl;
103  if ( imctruth >= nmctruth ) return false;
104 
105  size_t indx_handle = (*indx_itr).first;
106  size_t indx_within = (*indx_itr).second;
107 
108  thisLabel = outlabels[indx_handle];
109 
110  art::Handle< std::vector<simb::MCTruth> > hvMCTruth = mclists[indx_handle];
111 
120  thisMCTruth = &(hvMCTruth->at(indx_within));
121 
122  // inefficient ... should only need to do this bit for every new
123  // Handle ...
124 
125  //art::FindOneP<recob::Hit> findSPToHits(fSpacePoints,evt,fSpacePointLabel);
126  //const art::Ptr<recob::Hit> hit = findSPToHits.at(iSP);
127 
128  try {
129  art::FindOneP<simb::GTruth> qgtruth(hvMCTruth,evt,outlabels[indx_handle]);
130  const art::Ptr<simb::GTruth> gtruthptr = qgtruth.at(indx_within);
131  thisGTruth = gtruthptr.get();
132  }
133  catch (...) {
134  // std::cerr << "no GTruth for this handle" << std::endl;
135  }
136 
137  try {
138  art::FindOneP<simb::MCFlux> qmcflux(hvMCTruth,evt,outlabels[indx_handle]);
139  const art::Ptr<simb::MCFlux> mcfluxptr = qmcflux.at(indx_within);
140  thisMCFlux = mcfluxptr.get();
141  }
142  catch (...) {
143  // std::cerr << "no MCFlux for this handle" << std::endl;
144  }
145 
146  try {
147  art::FindOneP<bsim::Dk2Nu> qdk2nu(hvMCTruth,evt,outlabels[indx_handle]);
148  const art::Ptr<bsim::Dk2Nu> dk2nuptr = qdk2nu.at(indx_within);
149  thisDk2Nu = dk2nuptr.get();
150  }
151  catch (...) {
152  // std::cerr << "no bsim::Dk2NU for this handle" << std::endl;
153  }
154 
155  try {
156  art::FindOneP<bsim::NuChoice> qnuchoice(hvMCTruth,evt,outlabels[indx_handle]);
157  const art::Ptr<bsim::NuChoice> nuchoiceptr = qnuchoice.at(indx_within);
158  thisNuChoice = nuchoiceptr.get();
159  }
160  catch (...) {
161  // std::cerr << "no bsim::NuChoice for this handle" << std::endl;
162  }
163 
164  //std::cerr << "Next() called ... seeing " << thisMCTruth
165  // << " " << thisGTruth << " " << thisMCFlux << std::endl;
166 
167  // so user code looks like
168  // evgb::MCTruthAndFriendsItr mcitr(evt,labels);
169  // while ( mcitr.Next() ) {
170  // const simb::MCTruth* amctruth = mcitr.GetMCTruth();
171  // const simb::GTruth* agtruth = mcitr.GetGTruth();
172  //...
173 
174  /*
175  // loop over lists
176  try {
177  //art::FindOneP<simb::GTruth> QueryG(mclistHandle,evt,handleLabel);
178  }
179  catch (art::Exception) {
180  //std::cerr << "no GTruth for " << handleLabel << std::endl;
181  }
182  */
183 
184  // move the iterator on
185  ++indx_itr;
186  return true;
187 }
std::vector< art::Handle< std::vector< simb::MCTruth > > > mclists
MCTruthAndFriendsItr(art::Event const &evtIn, std::vector< std::string > const &labels)
std::vector< std::string > outlabels
bool isValid() const noexcept
Definition: Handle.h:203
const simb::GTruth * thisGTruth
Provenance const * provenance() const
Definition: Handle.h:217
constexpr std::array< std::size_t, geo::vect::dimension< Vector >)> indices()
Returns a sequence of indices valid for a vector of the specified type.
std::vector< std::string > const & fInputModuleLabels
const bsim::NuChoice * thisNuChoice
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
std::set< std::pair< int, int > > indices
const simb::MCFlux * thisMCFlux
const simb::MCTruth * thisMCTruth
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:69
TCEvent evt
Definition: DataStructs.cxx:8
std::set< std::pair< int, int > >::const_iterator indx_itr
T const * get() const
Definition: Ptr.h:138
std::vector< Handle< PROD > > getMany(SelectorBase const &selector=MatchAllSelector{}) const