MCTruthToDk2NuHackItr.cxx
Go to the documentation of this file.
2 
3 #include <iostream>
4 #include <iomanip>
5 
6 #ifndef ART_V1
9 #else
12 #endif
13 
14 #include <float.h> // FLT_EPSILON and DBL_EPSILON definitions used for
15  // floating point comparisons
16 
17 #include <cmath> // for frexp(), ldexp()
18 
20  std::vector<std::string> const & labels,
21  int verboseIn)
22  : evt(evtIn)
23  , fInputModuleLabels(labels)
24  , indx_itr(indices.begin())
25  , nmctruth(0)
26  , imctruth(-1)
27  , thisMCTruth(0)
28  , thisGTruth(0)
29  , thisMCFlux(0)
30  , thisDk2Nu(0)
31  , thisNuChoice(0)
32  , thisLabel("")
33  , verbose(verboseIn)
34 {
35 
36  // look for any existing MCTruth info in this event
37  mclists.clear();
38 
39  if ( fInputModuleLabels.size()==0) {
41  // std::cout << "evt.getManyByType" << std::endl;
42  } else {
43  mclists.resize(fInputModuleLabels.size());
44  for (size_t i=0; i<fInputModuleLabels.size(); ++i) {
46  //std::cout << "evt.getByLabel " << fInputModuleLabels[i] << std::endl;
47  }
48  }
49 
50  for (size_t mcl = 0; mcl < mclists.size(); ++mcl) {
52  if ( ! mclistHandle.isValid() ) {
53  // std::cerr << "not valid mcl=" << mcl << "---------------- " << std::endl;
54  continue;
55  }
56  // processName, moduleLabel, instance (productInstanceName?)
57  /*
58  std::cout << " mclistHandle processName '"
59  << mclistHandle.provenance()->processName() // top of fcl file
60  << "' moduleLabel '"
61  << mclistHandle.provenance()->moduleLabel()
62  << "' productInstanceName '"
63  << mclistHandle.provenance()->productInstanceName()
64  << "'"
65  << std::endl;
66  */
67 
68  std::string handleLabel = mclistHandle.provenance()->moduleLabel();
69  outlabels.push_back(handleLabel);
70  /*
71  std::cerr << "mcl=" << mcl << " '" << handleLabel << "' ---------------- " << std::endl;
72  */
73 
74  // loop over mctruths in a list
75  for(size_t nmc = 0; nmc < mclistHandle->size(); ++nmc) {
76  art::Ptr<simb::MCTruth> mct(mclistHandle, nmc);
77 
78  std::pair<int,int> ipair(mcl,nmc);
79  indices.insert(ipair);
80 
81  /**
82  std::cout << "+++ mcl " << mcl << "[" << mclists.size() << "] "
83  << "nmc " << nmc << "[" << mclistHandle->size() << "] "
84  << std::endl;
85  std::cout << *(mct.get()) << std::endl;
86  ***/
87  }
88 
89  }
90 
91  indx_itr = indices.begin();
92  nmctruth = (int)(indices.size());
93  if ( verbose > 0 ) {
94  std::cout << ".... found nmctruth " << nmctruth
95  << " nlabels " << outlabels.size()
96  << std::endl;
97  }
98 }
99 
101 {
102  ++imctruth; // started at -1, so first call to Next() prepares us for indx=0
103  thisMCTruth = 0;
104  thisGTruth = 0;
105  thisMCFlux = 0;
106  thisDk2Nu = 0;
107  thisNuChoice = 0;
108  //std::cerr << "Next() called ... moved to imctruth " << imctruth << std::endl;
109  if ( imctruth >= nmctruth ) return false;
110 
111  size_t indx_handle = (*indx_itr).first;
112  size_t indx_within = (*indx_itr).second;
113 
114  thisLabel = outlabels[indx_handle];
115 
116  art::Handle< std::vector<simb::MCTruth> > hvMCTruth = mclists[indx_handle];
117 
118  /**
119  std::cout << "imctruth " << std::setw(3) << imctruth
120  << " [" << indx_handle << "," << indx_within << "]"
121  << " hvMCTruth.isValid() " << hvMCTruth.isValid()
122  << " '" << outlabels[indx_handle] << "' "
123  << std::endl;
124  **/
125 
126  thisMCTruth = &(hvMCTruth->at(indx_within));
127  // also need art::Ptr<simb::MCTruth> ... for the association
128  // this should be checked ...
129  thisMCTruthPtr = art::Ptr<simb::MCTruth>(hvMCTruth,indx_within);
130 
131  // inefficient ... should only need to do this bit for every new
132  // Handle ...
133 
134  //art::FindOneP<recob::Hit> findSPToHits(fSpacePoints,evt,fSpacePointLabel);
135  //const art::Ptr<recob::Hit> hit = findSPToHits.at(iSP);
136 
137  try {
138  art::FindOneP<simb::GTruth> qgtruth(hvMCTruth,evt,outlabels[indx_handle]);
139  const art::Ptr<simb::GTruth> gtruthptr = qgtruth.at(indx_within);
140  thisGTruth = gtruthptr.get();
141  }
142  catch (...) {
143  // std::cerr << "no GTruth for this handle" << std::endl;
144  }
145 
146  try {
147  art::FindOneP<simb::MCFlux> qmcflux(hvMCTruth,evt,outlabels[indx_handle]);
148  const art::Ptr<simb::MCFlux> mcfluxptr = qmcflux.at(indx_within);
149  thisMCFlux = mcfluxptr.get();
150  }
151  catch (...) {
152  // std::cerr << "no MCFlux for this handle" << std::endl;
153  }
154 
155  int nuchoice_indx1 = -1;
156  int nuchoice_indx2 = -1;
157  try {
158  art::FindOneP<bsim::NuChoice> qnuchoice(hvMCTruth,evt,outlabels[indx_handle]);
159  const art::Ptr<bsim::NuChoice> nuchoiceptr = qnuchoice.at(indx_within);
160  thisNuChoice = nuchoiceptr.get();
161  }
162  catch (...) {
163  // except this is the hack ... no art::association
164  // assumes NuChoices put in w/ same label as MCTruth's they're matching to
165  //RWH//art::Handle<std::vector<bsim::NuChoice>> nuchoices;
166  //RWH//evt.getByLabel(thisLabel,nuchoices);
167  // allow for label skew induced by mixing
168  std::vector<art::Handle<std::vector<bsim::NuChoice>>> nuchoices2;
169  evt.getManyByType(nuchoices2);
170  if ( verbose > 2 ) {
171  std::cout << ".... getManyByType vector of " << nuchoices2.size()
172  << std::endl;
173  }
174 
175  for (size_t nuc2=0; nuc2 < nuchoices2.size(); ++nuc2) {
176  art::Handle<std::vector<bsim::NuChoice>> nuchoices = nuchoices2[nuc2];
177  // okay look through these NuChoices for a match
178  // our criteria will be a good match on the 4-vector
179 
180  size_t nnuchoices = nuchoices->size();
181  if ( verbose > 2 ) {
182  std::cout << ".... found n(NuChoice)s " << nnuchoices
183  << " for [" << nuc2 << "]" << std::endl;
184  }
185  for (size_t inuchoices=0; inuchoices<nnuchoices; ++inuchoices) {
186  // Handle->at() gives Ptr, & gets the object pointer?
187  thisNuChoice = &(nuchoices->at(inuchoices));
189  nuchoice_indx1 = (int)nuc2;
190  nuchoice_indx2 = (int)inuchoices;
191  if ( verbose > 2 ) {
192  std::cout << ".... found NuChoice match [" << nuchoice_indx1
193  << "," << nuchoice_indx2 << "] for MCTruth @ "
194  << thisMCTruth << std::endl;
195  }
196  break;
197  }
198  thisNuChoice = 0; // wasn't last one checked
199  }
200  } // end-of-loop nuc2
201  if ( nuchoice_indx1 < 0 || nuchoice_indx2 < 0 ) {
202  if ( verbose > 1 ) {
203  std::cerr << "no bsim::NuChoice for this handle via art::association"
204  << std::endl
205  << "also no good p4 match to MCTruth @ " << thisMCTruth
206  << std::endl;
207  }
208  }
209  }
210 
211  int dk2nu_indx1 = -1;
212  int dk2nu_indx2 = -1;
213  try {
214  art::FindOneP<bsim::Dk2Nu> qdk2nu(hvMCTruth,evt,outlabels[indx_handle]);
215  const art::Ptr<bsim::Dk2Nu> dk2nuptr = qdk2nu.at(indx_within);
216  thisDk2Nu = dk2nuptr.get();
217  }
218  catch (...) {
219  // except this is the hack ... no art::association
220  // assumes Dk2Nu put in w/ same label as MCTruth's they're matching to
221  //RWH//art::Handle<std::vector<bsim::Dk2Nu>> dk2nus;
222  //RWH//evt.getByLabel(thisLabel,dk2nus);
223  // allow for label skew induced by mixing
224  std::vector<art::Handle<std::vector<bsim::Dk2Nu>>> dk2nus2;
225  evt.getManyByType(dk2nus2);
226 
227  for (size_t dk2=0; dk2 < dk2nus2.size(); ++dk2) {
228  art::Handle<std::vector<bsim::Dk2Nu>> dk2nus = dk2nus2[dk2];
229  // okay look through these Dk2Nus for a match
230  // our criteria will be a match to MCFlux
231 
232  // could be more than one dk2nu entry for same job + proton #
233  // e.g. muon decay gives two nu's, or multiple nu out of shower
234  // job("run") & proton # & ndecay & ptype & nimpwt
235  // (can't use ntype as "oscillation" or flavor swap may have occurred)
236  // or fall back to the above index for nuchoice
237 
238  if ( dk2nu_indx1 >= 0 && dk2nu_indx2 >= 0 ) break;
239  if ( thisMCFlux ) {
240  // there is a MCFlux associated with the MCTruth
241  // this should have been filled via a simple copy of elements
242  size_t ndk2nus = dk2nus->size();
243  for (size_t idk2nus=0; idk2nus<ndk2nus; ++idk2nus) {
244  // Handle->at() gives Ptr, & gets the object pointer?
245  thisDk2Nu = &(dk2nus->at(idk2nus));
246  bool matched_dk2nu = true;
247  // our critera ... fail any and we're done
248  matched_dk2nu &= ( thisMCFlux->frun == thisDk2Nu->job );
249  matched_dk2nu &= ( thisMCFlux->fevtno == thisDk2Nu->potnum );
250  matched_dk2nu &= ( thisMCFlux->fndecay == thisDk2Nu->decay.ndecay );
251  matched_dk2nu &= ( thisMCFlux->fptype == thisDk2Nu->decay.ptype );
252  matched_dk2nu &=
254 
255  if ( matched_dk2nu ) {
256  dk2nu_indx1 = (int)dk2;
257  dk2nu_indx2 = (int)idk2nus;
258  break;
259  }
260  thisDk2Nu = 0; // wasn't the last one checked
261  }
262  }
263  } // end-of-loop nuc2
264 
265  if ( dk2nu_indx1 < 0 && dk2nu_indx2 < 0 ) {
266  if ( nuchoice_indx1 >= 0 && nuchoice_indx2 >= 0 ) {
267  // couldn't match Dk2Nu to MCFlux above, but do have a NuChoice match
268  // use the same index we found for NuChoice
269  dk2nu_indx1 = nuchoice_indx1;
270  dk2nu_indx2 = nuchoice_indx2;
271  thisDk2Nu = &((dk2nus2[dk2nu_indx1])->at(dk2nu_indx2));
272  }
273  }
274 
275  if ( dk2nu_indx1 < 0 || dk2nu_indx2 < 0 ) {
276  std::cerr << "no bsim::Dk2Nu for this handle via art::association"
277  << std::endl
278  << "also no good match to MCTruth/MCFlux/NuChoice"
279  << std::endl;
280  }
281  }
282 
283  //std::cerr << "Next() called ... seeing " << thisMCTruth
284  // << " " << thisGTruth << " " << thisMCFlux << std::endl;
285 
286  // so user code looks like
287  // evg::MCTruthToDk2NuHackItr mcitr(evt,labels);
288  // while ( mcitr.Next() ) {
289  // const simb::MCTruth* amctruth = mcitr.GetMCTruth();
290  // const simb::GTruth* agtruth = mcitr.GetGTruth();
291  //...
292 
293  /*
294  // loop over lists
295  try {
296  //art::FindOneP<simb::GTruth> QueryG(mclistHandle,evt,handleLabel);
297  }
298  catch (art::Exception) {
299  //std::cerr << "no GTruth for " << handleLabel << std::endl;
300  }
301  */
302 
303  // move the iterator on
304  ++indx_itr;
305  return true;
306 }
307 
308 bool evg::MCTruthToDk2NuHackItr::pretty_darn_close(double a, double b, double eps_precision)
309 {
310  // based in part on Sue Kasahara's implementaion in MINOS compareNtuple.C
311  // floating point comparisons are done with an adaptation of T. Belding's
312  // fcmp implementation of a D. Knuth algorithm.
313 
314  // For float or double, need to determine equality while
315  // accounting for some degree of imprecision
316  // The method applied here is adopted from Theodore C. Belding's
317  // fcmp implementation of Donald Knuth's algorithm.
318  double value[2] = {a,b};
319  double epsilon = eps_precision*FLT_EPSILON;
320  int exponent;
321  std::frexp(std::fabs(value[0]) > std::fabs(value[1]) ? value[0] : value[1], &exponent);
322  double delta = std::ldexp(epsilon,exponent);
323  double difference = value[0] - value[1];
324  if ( difference > delta || difference < -delta ) return false;
325  return true;
326 }
327 
328 
330  const bsim::NuChoice* pnuchoice)
331 {
332 
333  const TLorentzVector& p4NuMCTruth = pmctruth->GetNeutrino().Nu().Momentum();
334  const TLorentzVector& p4NuNuChoice = pnuchoice->p4NuUser;
335 
336  for (int indx=0; indx<4; ++indx) {
337  if ( ! pretty_darn_close(p4NuMCTruth[indx],p4NuNuChoice[indx]) ) return false;
338  }
339  return true; // all 4 components are the same
340 
341 }
int frun
Definition: MCFlux.h:35
Int_t job
identifying job #
Definition: dk2nu.h:323
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
std::vector< std::string > const & fInputModuleLabels
double delta
Definition: runWimpSim.h:98
std::set< std::pair< int, int > >::const_iterator indx_itr
Int_t potnum
proton # processed by simulation
Definition: dk2nu.h:324
const simb::MCTruth * thisMCTruth
std::vector< art::Handle< std::vector< simb::MCTruth > > > mclists
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
T ldexp(const T &a, int b)
Definition: ldexp.hpp:19
Double_t nimpwt
% cumulative importance weight prod to decay
Definition: dk2nu.h:166
OStream cerr
Definition: OStream.cxx:7
bsim::Decay decay
basic decay information
Definition: dk2nu.h:327
std::vector< std::string > outlabels
bool isValid() const
Definition: Handle.h:189
const bsim::NuChoice * thisNuChoice
int fptype
Definition: MCFlux.h:63
const XML_Char int const XML_Char * value
Definition: expat.h:331
int fndecay
Definition: MCFlux.h:50
Provenance const * provenance() const
Definition: Handle.h:203
const double a
void getManyByType(std::vector< Handle< PROD >> &results) const
Definition: DataViewImpl.h:446
bool match_mctruth_nuchoice(const simb::MCTruth *pmctruth, const bsim::NuChoice *pnchoice)
int fevtno
Definition: MCFlux.h:36
bool pretty_darn_close(double a, double b, double eps_precision=10.0)
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
art::Ptr< simb::MCTruth > thisMCTruthPtr
T const * get() const
Definition: Ptr.h:321
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
TLorentzVector p4NuUser
generated nu 4-momentum in user/det coord
Definition: NuChoice.h:57
std::set< std::pair< int, int > > indices
double epsilon
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
const hit & b
Definition: hits.cxx:21
Int_t ndecay
decay process (see dkproc_t)
Definition: dk2nu.h:127
const simb::MCFlux * thisMCFlux
Event generator information.
Definition: MCTruth.h:32
const simb::GTruth * thisGTruth
Int_t ptype
% nu parent species (PDG? code)
Definition: dk2nu.h:144
MCTruthToDk2NuHackItr(art::Event const &evtIn, std::vector< std::string > const &labels, int verbose=0)
double fnimpwt
Definition: MCFlux.h:72