PathSegmentList.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Robert Hatcher <rhatcher@fnal.gov>
8  FNAL - May 26, 2009
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ May 26, 2009 - RWH
14  First added in v2.5.1
15  @ July 16, 2009 - RWH
16  Reworked ROOTGeomAnalyzer code so that each neutrino from the flux is only
17  swum through the geometry once (not once per tgt PDG + 2 more to pick vertex
18  if interaction is selected). The PathSegmentList class holds all the individual
19  steps taken and is amenable to the next step which is trimming the list to
20  restrict the interactions to a non-physically realized volume.
21  @ July 27, 2009 - RWH
22  Print() also shows start pos and dir
23  @ August 3, 2009 - RWH
24  Individual path segments can be further trimmed to be less than the full
25  step (which corresponds to a volume).
26  The list is now internally a STL list (no longer a STL vector).
27  This is in anticipation of a possible need for segment splitting if
28  trimming causes disjoint segments.
29  Both classes have provisions for cross checking ray distance and step size
30  values after completion of the list. Generally I've found ray distance
31  errors of no more than ~8e-9 cm (80000fm) in a 42000cm swim through a
32  geometry. This is on order an atomic size so one shouldn't worry about it.
33  This precision will allow trimming to be based on the distance travelled
34  along the ray.
35  @ August 8, 2009 - RWH
36  Push PathSegmentList down a namespace so it's in genie::geometry.
37  Fix const'ness of some functions. Add SetPath() function and fPathString
38  data member for recording geometry path info in case users want to cut on
39  something there.
40  @ February 4, 2010 - RWH
41  Correct statement about how overhead much fetching geometry path adds
42  (significant, not negligable).
43  Generalize from (lo,hi) pair to vector of (lo,hi) pairs to allow the
44  geometry step to be split and not just squeezed. This might not be
45  necessary and how much this adds to overhead isn't quite clear
46  (needs more testing).
47  New methods IsTrimmedEmpty() and GetSummedStepRange().
48  Also GetPosition() to pick position within vector of (lo,hi) pairs based on
49  fraction of total. Needs more testing for case of split segments.
50 
51 */
52 //____________________________________________________________________________
53 
54 #include <fstream>
55 #include <iomanip>
56 #include <sstream>
57 #include <cassert>
58 
59 //#include "libxml/parser.h"
60 //#include "libxml/xmlmemory.h"
61 
62 #include <TLorentzVector.h>
63 #include <TGeoVolume.h>
64 #include <TGeoMaterial.h>
65 
72 //#include "Framework/Utils/XmlParserUtils.h"
73 
74 using std::ofstream;
75 using std::setw;
76 using std::setfill;
77 using std::endl;
78 
79 using namespace genie;
80 using namespace genie::geometry;
81 
82 //#define RWH_DEBUG
83 
84 //____________________________________________________________________________
85 namespace genie {
86 namespace geometry {
87 
88 ostream & operator << (ostream & stream, const geometry::PathSegment & ps)
89  {
90  ps.Print(stream);
91  return stream;
92  }
93 
94 ostream & operator << (ostream & stream, const geometry::PathSegmentList & list)
95  {
96  list.Print(stream);
97  return stream;
98  }
99 
100 } // namespace geometry
101 } // namespace genie
102 
103 //____________________________________________________________________________
104 namespace genie {
105  namespace pathsegutils {
106  string Vec3AsString (const TVector3 * vec)
107  {
108  int w=17, p=10; // precision setting only affects ostringstream
109  std::ostringstream fmt;
110  fmt << "(" << std::setw(w) << std::setprecision(p) << vec->x()
111  << "," << std::setw(w) << std::setprecision(p) << vec->y()
112  << "," << std::setw(w+1) << std::setprecision(p) << vec->z() << ")";
113  return fmt.str();
114  }
115  }
116 }
117 
118 //___________________________________________________________________________
120  fRayDist(0), fStepLength(0),
121  fVolume(0), fMedium(0), fMaterial(0),
122  fEnter(), fExit()
123 {
124 }
125 
126 //___________________________________________________________________________
127 void PathSegment::DoCrossCheck(const TVector3& startpos,
128  double& ddist, double& dstep) const
129 {
130  double dist_recalc = (fEnter-startpos).Mag();
131  ddist = dist_recalc - fRayDist;
132 
133  double step_recalc = (fExit-fEnter).Mag();
134  dstep = step_recalc - fStepLength;
135 }
136 
137 //___________________________________________________________________________
138 void PathSegment::Print(ostream & stream) const
139 {
140  const char* vname = (fVolume) ? fVolume->GetName() : "no volume";
141  const char* mname = (fMaterial) ? fMaterial->GetName() : "no material";
142  stream << genie::pathsegutils::Vec3AsString(&fEnter) << " "
143  //<< genie::pathsegutils::Vec3AsString(&fExit)
144  << " " // "raydist "
145  << std::setw(12) << fRayDist
146  << " " // "step "
147  << std::setw(12) << fStepLength << " "
148  << std::left
149  << std::setw(16) << vname << " '"
150  << std::setw(18) << mname << "' ";
151  size_t n = fStepRangeSet.size();
152  const int rngw = 24;
153  if ( n == 0 ) {
154  stream << std::setw(rngw) << "[ ]";
155  } else {
156  std::ostringstream rngset;
157  for ( size_t i = 0 ; i < n; ++i ) {
158  const StepRange& sr = fStepRangeSet[i];
159  rngset << "[" << sr.first << ":" << sr.second << "]";
160  }
161  stream << std::setw(rngw) << rngset.str();
162  }
163  stream << std::right;
164 #ifdef PATHSEG_KEEP_PATH
165  stream << " " << fPathString;
166 #endif
167 
168 }
169 
170 //___________________________________________________________________________
171 void PathSegment::SetStep(Double_t step, bool setlimits)
172 {
173  fStepLength = step;
174  if (setlimits) {
175  fStepRangeSet.clear();
176  fStepRangeSet.push_back(StepRange(0,step));
177  }
178 }
179 
180 //___________________________________________________________________________
182 {
183  Double_t sum = 0;
184  for ( size_t i = 0; i < fStepRangeSet.size(); ++i ) {
185  const StepRange& sr = fStepRangeSet[i];
186  sum += ( sr.second - sr.first );
187  }
188  return sum;
189 }
190 
191 TVector3 PathSegment::GetPosition(Double_t fractrim) const
192 {
193  /// calculate position within allowed ranges passed as
194  /// fraction of trimmed segment
195  /// seg.fEnter + fractotal * ( seg.fExit - seg.fEnter );
196  Double_t sumrange = GetSummedStepRange();
197  if ( sumrange < 0.0 ) {
198  LOG("PathS", pFATAL) << "GetPosition failed fractrim=" << fractrim
199  << " because sumrange = " << sumrange;
200  return TVector3(0,0,0);
201  }
202  Double_t target = fractrim * sumrange;
203  Double_t sum = 0;
204  for ( size_t i = 0; i < fStepRangeSet.size(); ++i ) {
205  const StepRange& sr = fStepRangeSet[i];
206  Double_t ds = ( sr.second - sr.first );
207  sum += ds;
208 #ifdef RWH_DEBUG
209  LOG("PathS", pINFO) << "GetPosition fractrim=" << fractrim
210  << " target " << target << " [" << i << "] "
211  << " ds " << ds << " sum " << sum;
212 #endif
213  if ( sum >= target ) {
214  Double_t overstep = sum - target;
215  Double_t fractotal = (sr.second - overstep)/fStepLength;
216 #ifdef RWH_DEBUG
217  LOG("PathS", pINFO) << "GetPosition fractrim=" << fractrim
218  << " overstep " << overstep
219  << " fractotal " << fractotal;
220 #endif
221  return fEnter + fractotal * ( fExit - fEnter );
222  }
223  }
224  LOG("PathS", pFATAL) << "GetPosition failed fractrim=" << fractrim;
225  return TVector3(0,0,0);
226 }
227 
228 
229 //===========================================================================
230 //___________________________________________________________________________
232  : fDoCrossCheck(false), fPrintVerbose(false)
233 {
234 
235 }
236 //___________________________________________________________________________
238 {
239  this->Copy(plist);
240 }
241 //__________________________________________________________________________
243 {
244 
245 }
246 
247 //___________________________________________________________________________
249 {
250  LOG("PathS", pDEBUG) << "SetAllToZero called";
251 
252  this->fStartPos.SetXYZ(0,0,1e37); // clear cache of position/direction
253  this->fDirection.SetXYZ(0,0,0); //
254  this->fSegmentList.clear(); // clear the vector
255  this->fMatStepSum.clear(); // clear the re-factorized info
256 }
257 
258 //___________________________________________________________________________
259 void PathSegmentList::SetStartInfo(const TVector3& pos, const TVector3& dir)
260 {
261  this->fStartPos = pos;
262  this->fDirection = dir;
263 }
264 
265 //___________________________________________________________________________
266 bool PathSegmentList::IsSameStart(const TVector3& pos, const TVector3& dir) const
267 {
268  return ( this->fStartPos == pos && this->fDirection == dir );
269 }
270 
271 //___________________________________________________________________________
273 {
274  fMatStepSum.clear();
275 
278  for ( ; sitr != sitr_end ; ++sitr ) {
279  const PathSegment& ps = *sitr;
280  const TGeoMaterial* mat = ps.fMaterial;
281  // use the post-trim limits on how much material is stepped through
283  }
284 
285 }
286 
287 //___________________________________________________________________________
289 {
290  fSegmentList.clear();
291  fMatStepSum.clear();
292 
293  // copy the segments
294  //vector<PathSegment>::const_iterator pl_iter;
295  //for (pl_iter = plist.fSegmentList.begin(); pl_iter != plist.fSegmentList.end(); ++pl_iter) {
296  // this->fSegmentList.push_back( *pl_iter );
297  //}
298 
299  // other elements
300  fStartPos = plist.fStartPos;
301  fDirection = plist.fDirection;
302  fSegmentList = plist.fSegmentList;
303  fMatStepSum = plist.fMatStepSum;
306 }
307 
308 //___________________________________________________________________________
309 void PathSegmentList::CrossCheck(double& mxddist, double& mxdstep) const
310 {
311 
312  double dstep, ddist;
313  mxdstep = 0;
314  mxddist = 0;
317  for ( ; sitr != sitr_end ; ++sitr ) {
318  const PathSegment& ps = *sitr;
319  ps.DoCrossCheck(fStartPos,ddist,dstep);
320  double addist = TMath::Abs(ddist);
321  double adstep = TMath::Abs(dstep);
322  if ( addist > mxddist ) mxddist = addist;
323  if ( adstep > mxdstep ) mxdstep = adstep;
324  }
325 
326 }
327 
328 //___________________________________________________________________________
329 void PathSegmentList::Print(ostream & stream) const
330 {
331  stream << "\nPathSegmentList [-]" << endl;
332  stream << " start " << pathsegutils::Vec3AsString(&fStartPos)
333  << " dir " << pathsegutils::Vec3AsString(&fDirection) << endl;
334 
335  double dstep, ddist, mxdstep = 0, mxddist = 0;
336  int k = 0, nseg = 0;
339  for ( ; sitr != sitr_end ; ++sitr, ++k ) {
340  const PathSegment& ps = *sitr;
341  ++nseg;
342  stream << " [" << setw(4) << k << "] " << ps;
343  if ( fDoCrossCheck ) {
344  ps.DoCrossCheck(fStartPos,ddist,dstep);
345  double addist = TMath::Abs(ddist);
346  double adstep = TMath::Abs(dstep);
347  if ( addist > mxddist ) mxddist = addist;
348  if ( adstep > mxdstep ) mxdstep = adstep;
349  stream << " recalc diff"
350  << " dist " << std::setw(12) << ddist
351  << " step " << std::setw(12) << dstep;
352  }
353  stream << std::endl;
354  }
355  if ( nseg == 0 ) stream << " holds no segments." << std::endl;
356 
357  if ( fDoCrossCheck )
358  stream << "PathSegmentList "
359  << " mxddist " << mxddist
360  << " mxdstep " << mxdstep
361  << std::endl;
362 
363  if ( fPrintVerbose ) {
366  // loop over map to get tgt weight for each material (once)
367  // steps outside the geometry may have no assigned material
368  for ( ; mitr != mitr_end; ++mitr ) {
369  const TGeoMaterial* mat = mitr->first;
370  double sumsteps = mitr->second;
371  stream << " fMatStepSum[" << mat->GetName() << "] = " << sumsteps << std::endl;
372  }
373  }
374 
375 }
376 //___________________________________________________________________________
377 #ifdef PATH_SEGMENT_SUPPORT_XML
378 XmlParserStatus_t PathSegmentList::LoadFromXml(string filename)
379 {
380  this->clear();
381  PDGLibrary * pdglib = PDGLibrary::Instance();
382 
383  LOG("PathS", pINFO)
384  << "Loading PathSegmentList from XML file: " << filename;
385 
386  xmlDocPtr xml_doc = xmlParseFile(filename.c_str() );
387 
388  if(xml_doc==NULL) {
389  LOG("PathS", pERROR)
390  << "XML file could not be parsed! [filename: " << filename << "]";
391  return kXmlNotParsed;
392  }
393 
394  xmlNodePtr xmlCur = xmlDocGetRootElement(xml_doc);
395 
396  if(xmlCur==NULL) {
397  LOG("PathS", pERROR)
398  << "XML doc. has null root element! [filename: " << filename << "]";
399  return kXmlEmpty;
400  }
401 
402  if( xmlStrcmp(xmlCur->name, (const xmlChar *) "path_length_list") ) {
403  LOG("PathS", pERROR)
404  << "XML doc. has invalid root element! [filename: " << filename << "]";
405  return kXmlInvalidRoot;
406  }
407 
408  LOG("PathS", pINFO) << "XML file was successfully parsed";
409 
410  xmlCur = xmlCur->xmlChildrenNode; // <path_length>'s
411 
412  // loop over all xml tree nodes that are children of the root node
413  while (xmlCur != NULL) {
414 
415  // enter everytime you find a <path_length> tag
416  if( (!xmlStrcmp(xmlCur->name, (const xmlChar *) "path_length")) ) {
417 
418  xmlNodePtr xmlPlVal = xmlCur->xmlChildrenNode;
419 
420  string spdgc = utils::str::TrimSpaces(
421  XmlParserUtils::GetAttribute(xmlCur, "pdgc"));
422 
423  string spl = XmlParserUtils::TrimSpaces(
424  xmlNodeListGetString(xml_doc, xmlPlVal, 1));
425 
426  LOG("PathS", pDEBUG) << "pdgc = " << spdgc << " --> pl = " << spl;
427 
428  int pdgc = atoi( spdgc.c_str() );
429  double pl = atof( spl.c_str() );
430 
431  TParticlePDG * p = pdglib->Find(pdgc);
432  if(!p) {
433  LOG("PathS", pERROR)
434  << "No particle with pdgc " << pdgc
435  << " found. Will not load its path length";
436  } else
437  this->insert( map<int, double>::value_type(pdgc, pl) );
438 
439  xmlFree(xmlPlVal);
440  }
441  xmlCur = xmlCur->next;
442  } // [end of] loop over tags within root elements
443 
444  xmlFree(xmlCur);
445  return kXmlOK;
446 }
447 //___________________________________________________________________________
448 void PathSegmentList::SaveAsXml(string filename) const
449 {
450 //! Save path length list to XML file
451 
452  LOG("PathS", pINFO)
453  << "Saving PathSegmentList as XML in file: " << filename;
454 
455  ofstream outxml(filename.c_str());
456  if(!outxml.is_open()) {
457  LOG("PathS", pERROR) << "Couldn't create file = " << filename;
458  return;
459  }
460  outxml << "<?xml version=\"1.0\" encoding=\"ISO-8859-1\"?>";
461  outxml << endl << endl;
462  outxml << "<!-- generated by PathSegmentList::SaveAsXml() -->";
463  outxml << endl << endl;
464 
465  outxml << "<path_length_list>" << endl << endl;
466 
467  PathSegmentList::const_iterator pl_iter;
468 
469  for(pl_iter = this->begin(); pl_iter != this->end(); ++pl_iter) {
470 
471  int pdgc = pl_iter->first;
472  double pl = pl_iter->second; // path length
473 
474  outxml << " <path_length pdgc=\"" << pdgc << "\">"
475  << pl << "</path_length>" << endl;
476  }
477  outxml << "</path_length_list>";
478  outxml << endl;
479 
480  outxml.close();
481 
482 }
483 #endif
484 //___________________________________________________________________________
486 {
487  this->Copy(list);
488  return (*this);
489 }
490 //___________________________________________________________________________
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:112
PathSegmentV_t fSegmentList
Actual list of segments.
const XML_Char * target
Definition: expat.h:268
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
void CrossCheck(double &mxddist, double &mxdstep) const
MaterialMap_t fMatStepSum
Segment list re-evaluated by material for fast lookup of path lengths.
string Vec3AsString(const TVector3 *vec)
StepRangeSet fStepRangeSet
collection of {steplo,stephi} pairs
TVector3 fEnter
top vol coordinates and units
const char * p
Definition: xmltok.h:285
#define pFATAL
Definition: Messenger.h:57
TVector3 fStartPos
Record, for future comparison, the path taken.
vector< vector< double > > clear
GENIE geometry drivers.
Object to be filled with the neutrino path-segments representing geometry volume steps (generally bou...
string filename
Definition: shutoffs.py:106
bool IsSameStart(const TVector3 &pos, const TVector3 &dir) const
void SetStep(Double_t step, bool setlimits=true)
step taken in the geometry element
const TGeoMaterial * fMaterial
ref only ptr to TGeoMaterial
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
PathSegmentV_t::const_iterator PathSegVCItr_t
void Print(ostream &stream) const
TVector3 fDirection
direction (in top vol coords)
const TGeoVolume * fVolume
ref only ptr to TGeoVolume
TVector3 fExit
top vol coordinates and units
void Copy(const PathSegmentList &plist)
caf::StandardRecord * sr
void SetStartInfo(const TVector3 &pos=TVector3(0, 0, 1e37), const TVector3 &dir=TVector3(0, 0, 0))
#define pINFO
Definition: Messenger.h:63
Eigen::VectorXd vec
Double_t GetSummedStepRange() const
get the sum of all the step range (in case step has been trimmed or split)
std::pair< Double_t, Double_t > StepRange
A very simple service to remember what detector we&#39;re working in.
string TrimSpaces(string input)
Definition: StringUtils.cxx:24
MaterialMap_t::const_iterator MaterialMapCItr_t
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:104
Double_t fRayDist
distance from start of ray
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
Singleton class to load & serve a TDatabasePDG.
Definition: PDGLibrary.h:30
TDirectory * dir
Definition: macro.C:5
TVector3 GetPosition(Double_t frac) const
calculate position within allowed ranges passed on fraction of total
void Print(ostream &stream) const
const MaterialMap_t & GetMatStepSumMap(void) const
PathSegmentList & operator=(const PathSegmentList &list)
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
float Mag() const
Double_t sum
Definition: plot.C:31
void DoCrossCheck(const TVector3 &startpos, double &ddist, double &dstep) const
perform cross check on segment, return differences
enum genie::EXmlParseStatus XmlParserStatus_t
string GetAttribute(xmlNodePtr xml_cur, string attr_name)
Float_t w
Definition: plot.C:20
std::ostream & operator<<(std::ostream &stream, const genie::geometry::PlaneParam &pparam)
Definition: FidShape.cxx:26
Double_t fStepLength
total step size in volume
Eigen::MatrixXd mat
#define pDEBUG
Definition: Messenger.h:64
std::string fPathString
full path names