GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PathSegmentList.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Robert Hatcher <rhatcher@fnal.gov>
7  FNAL
8 */
9 //____________________________________________________________________________
10 
11 #include <fstream>
12 #include <iomanip>
13 #include <sstream>
14 #include <cassert>
15 
16 //#include "libxml/parser.h"
17 //#include "libxml/xmlmemory.h"
18 
19 #include <TLorentzVector.h>
20 #include <TGeoVolume.h>
21 #include <TGeoMaterial.h>
22 
29 //#include "Framework/Utils/XmlParserUtils.h"
30 
31 using std::ofstream;
32 using std::setw;
33 using std::setfill;
34 using std::endl;
35 
36 using namespace genie;
37 using namespace genie::geometry;
38 
39 //#define RWH_DEBUG
40 
41 //____________________________________________________________________________
42 namespace genie {
43 namespace geometry {
44 
45 ostream & operator << (ostream & stream, const geometry::PathSegment & ps)
46  {
47  ps.Print(stream);
48  return stream;
49  }
50 
51 ostream & operator << (ostream & stream, const geometry::PathSegmentList & list)
52  {
53  list.Print(stream);
54  return stream;
55  }
56 
57 } // namespace geometry
58 } // namespace genie
59 
60 //____________________________________________________________________________
61 namespace genie {
62  namespace pathsegutils {
63  string Vec3AsString (const TVector3 * vec)
64  {
65  int w=17, p=10; // precision setting only affects ostringstream
66  std::ostringstream fmt;
67  fmt << "(" << std::setw(w) << std::setprecision(p) << vec->x()
68  << "," << std::setw(w) << std::setprecision(p) << vec->y()
69  << "," << std::setw(w+1) << std::setprecision(p) << vec->z() << ")";
70  return fmt.str();
71  }
72  }
73 }
74 
75 //___________________________________________________________________________
76 PathSegment::PathSegment(void) :
77  fRayDist(0), fStepLength(0),
78  fVolume(0), fMedium(0), fMaterial(0),
79  fEnter(), fExit()
80 {
81 }
82 
83 //___________________________________________________________________________
84 void PathSegment::DoCrossCheck(const TVector3& startpos,
85  double& ddist, double& dstep) const
86 {
87  double dist_recalc = (fEnter-startpos).Mag();
88  ddist = dist_recalc - fRayDist;
89 
90  double step_recalc = (fExit-fEnter).Mag();
91  dstep = step_recalc - fStepLength;
92 }
93 
94 //___________________________________________________________________________
95 void PathSegment::Print(ostream & stream) const
96 {
97  const char* vname = (fVolume) ? fVolume->GetName() : "no volume";
98  const char* mname = (fMaterial) ? fMaterial->GetName() : "no material";
99  stream << genie::pathsegutils::Vec3AsString(&fEnter) << " "
100  //<< genie::pathsegutils::Vec3AsString(&fExit)
101  << " " // "raydist "
102  << std::setw(12) << fRayDist
103  << " " // "step "
104  << std::setw(12) << fStepLength << " "
105  << std::left
106  << std::setw(16) << vname << " '"
107  << std::setw(18) << mname << "' ";
108  size_t n = fStepRangeSet.size();
109  const int rngw = 24;
110  if ( n == 0 ) {
111  stream << std::setw(rngw) << "[ ]";
112  } else {
113  std::ostringstream rngset;
114  for ( size_t i = 0 ; i < n; ++i ) {
115  const StepRange& sr = fStepRangeSet[i];
116  rngset << "[" << sr.first << ":" << sr.second << "]";
117  }
118  stream << std::setw(rngw) << rngset.str();
119  }
120  stream << std::right;
121 #ifdef PATHSEG_KEEP_PATH
122  stream << " " << fPathString;
123 #endif
124 
125 }
126 
127 //___________________________________________________________________________
128 void PathSegment::SetStep(Double_t step, bool setlimits)
129 {
130  fStepLength = step;
131  if (setlimits) {
132  fStepRangeSet.clear();
133  fStepRangeSet.push_back(StepRange(0,step));
134  }
135 }
136 
137 //___________________________________________________________________________
139 {
140  Double_t sum = 0;
141  for ( size_t i = 0; i < fStepRangeSet.size(); ++i ) {
142  const StepRange& sr = fStepRangeSet[i];
143  sum += ( sr.second - sr.first );
144  }
145  return sum;
146 }
147 
148 TVector3 PathSegment::GetPosition(Double_t fractrim) const
149 {
150  /// calculate position within allowed ranges passed as
151  /// fraction of trimmed segment
152  /// seg.fEnter + fractotal * ( seg.fExit - seg.fEnter );
153  Double_t sumrange = GetSummedStepRange();
154  if ( sumrange < 0.0 ) {
155  LOG("PathS", pFATAL) << "GetPosition failed fractrim=" << fractrim
156  << " because sumrange = " << sumrange;
157  return TVector3(0,0,0);
158  }
159  Double_t target = fractrim * sumrange;
160  Double_t sum = 0;
161  for ( size_t i = 0; i < fStepRangeSet.size(); ++i ) {
162  const StepRange& sr = fStepRangeSet[i];
163  Double_t ds = ( sr.second - sr.first );
164  sum += ds;
165 #ifdef RWH_DEBUG
166  LOG("PathS", pINFO) << "GetPosition fractrim=" << fractrim
167  << " target " << target << " [" << i << "] "
168  << " ds " << ds << " sum " << sum;
169 #endif
170  if ( sum >= target ) {
171  Double_t overstep = sum - target;
172  Double_t fractotal = (sr.second - overstep)/fStepLength;
173 #ifdef RWH_DEBUG
174  LOG("PathS", pINFO) << "GetPosition fractrim=" << fractrim
175  << " overstep " << overstep
176  << " fractotal " << fractotal;
177 #endif
178  return fEnter + fractotal * ( fExit - fEnter );
179  }
180  }
181  LOG("PathS", pFATAL) << "GetPosition failed fractrim=" << fractrim;
182  return TVector3(0,0,0);
183 }
184 
185 
186 //===========================================================================
187 //___________________________________________________________________________
189  : fDoCrossCheck(false), fPrintVerbose(false)
190 {
191 
192 }
193 //___________________________________________________________________________
195 {
196  this->Copy(plist);
197 }
198 //__________________________________________________________________________
200 {
201 
202 }
203 
204 //___________________________________________________________________________
206 {
207  LOG("PathS", pDEBUG) << "SetAllToZero called";
208 
209  this->fStartPos.SetXYZ(0,0,1e37); // clear cache of position/direction
210  this->fDirection.SetXYZ(0,0,0); //
211  this->fSegmentList.clear(); // clear the vector
212  this->fMatStepSum.clear(); // clear the re-factorized info
213 }
214 
215 //___________________________________________________________________________
216 void PathSegmentList::SetStartInfo(const TVector3& pos, const TVector3& dir)
217 {
218  this->fStartPos = pos;
219  this->fDirection = dir;
220 }
221 
222 //___________________________________________________________________________
223 bool PathSegmentList::IsSameStart(const TVector3& pos, const TVector3& dir) const
224 {
225  return ( this->fStartPos == pos && this->fDirection == dir );
226 }
227 
228 //___________________________________________________________________________
230 {
231  fMatStepSum.clear();
232 
235  for ( ; sitr != sitr_end ; ++sitr ) {
236  const PathSegment& ps = *sitr;
237  const TGeoMaterial* mat = ps.fMaterial;
238  // use the post-trim limits on how much material is stepped through
239  fMatStepSum[mat] += ps.GetSummedStepRange();
240  }
241 
242 }
243 
244 //___________________________________________________________________________
246 {
247  fSegmentList.clear();
248  fMatStepSum.clear();
249 
250  // copy the segments
251  //vector<PathSegment>::const_iterator pl_iter;
252  //for (pl_iter = plist.fSegmentList.begin(); pl_iter != plist.fSegmentList.end(); ++pl_iter) {
253  // this->fSegmentList.push_back( *pl_iter );
254  //}
255 
256  // other elements
257  fStartPos = plist.fStartPos;
258  fDirection = plist.fDirection;
259  fSegmentList = plist.fSegmentList;
260  fMatStepSum = plist.fMatStepSum;
263 }
264 
265 //___________________________________________________________________________
266 void PathSegmentList::CrossCheck(double& mxddist, double& mxdstep) const
267 {
268 
269  double dstep, ddist;
270  mxdstep = 0;
271  mxddist = 0;
274  for ( ; sitr != sitr_end ; ++sitr ) {
275  const PathSegment& ps = *sitr;
276  ps.DoCrossCheck(fStartPos,ddist,dstep);
277  double addist = TMath::Abs(ddist);
278  double adstep = TMath::Abs(dstep);
279  if ( addist > mxddist ) mxddist = addist;
280  if ( adstep > mxdstep ) mxdstep = adstep;
281  }
282 
283 }
284 
285 //___________________________________________________________________________
286 void PathSegmentList::Print(ostream & stream) const
287 {
288  stream << "\nPathSegmentList [-]" << endl;
289  stream << " start " << pathsegutils::Vec3AsString(&fStartPos)
290  << " dir " << pathsegutils::Vec3AsString(&fDirection) << endl;
291 
292  double dstep, ddist, mxdstep = 0, mxddist = 0;
293  int k = 0, nseg = 0;
296  for ( ; sitr != sitr_end ; ++sitr, ++k ) {
297  const PathSegment& ps = *sitr;
298  ++nseg;
299  stream << " [" << setw(4) << k << "] " << ps;
300  if ( fDoCrossCheck ) {
301  ps.DoCrossCheck(fStartPos,ddist,dstep);
302  double addist = TMath::Abs(ddist);
303  double adstep = TMath::Abs(dstep);
304  if ( addist > mxddist ) mxddist = addist;
305  if ( adstep > mxdstep ) mxdstep = adstep;
306  stream << " recalc diff"
307  << " dist " << std::setw(12) << ddist
308  << " step " << std::setw(12) << dstep;
309  }
310  stream << std::endl;
311  }
312  if ( nseg == 0 ) stream << " holds no segments." << std::endl;
313 
314  if ( fDoCrossCheck )
315  stream << "PathSegmentList "
316  << " mxddist " << mxddist
317  << " mxdstep " << mxdstep
318  << std::endl;
319 
320  if ( fPrintVerbose ) {
323  // loop over map to get tgt weight for each material (once)
324  // steps outside the geometry may have no assigned material
325  for ( ; mitr != mitr_end; ++mitr ) {
326  const TGeoMaterial* mat = mitr->first;
327  double sumsteps = mitr->second;
328  stream << " fMatStepSum[" << mat->GetName() << "] = " << sumsteps << std::endl;
329  }
330  }
331 
332 }
333 //___________________________________________________________________________
334 #ifdef PATH_SEGMENT_SUPPORT_XML
335 XmlParserStatus_t PathSegmentList::LoadFromXml(string filename)
336 {
337  this->clear();
338  PDGLibrary * pdglib = PDGLibrary::Instance();
339 
340  LOG("PathS", pINFO)
341  << "Loading PathSegmentList from XML file: " << filename;
342 
343  xmlDocPtr xml_doc = xmlParseFile(filename.c_str() );
344 
345  if(xml_doc==NULL) {
346  LOG("PathS", pERROR)
347  << "XML file could not be parsed! [filename: " << filename << "]";
348  return kXmlNotParsed;
349  }
350 
351  xmlNodePtr xmlCur = xmlDocGetRootElement(xml_doc);
352 
353  if(xmlCur==NULL) {
354  LOG("PathS", pERROR)
355  << "XML doc. has null root element! [filename: " << filename << "]";
356  return kXmlEmpty;
357  }
358 
359  if( xmlStrcmp(xmlCur->name, (const xmlChar *) "path_length_list") ) {
360  LOG("PathS", pERROR)
361  << "XML doc. has invalid root element! [filename: " << filename << "]";
362  return kXmlInvalidRoot;
363  }
364 
365  LOG("PathS", pINFO) << "XML file was successfully parsed";
366 
367  xmlCur = xmlCur->xmlChildrenNode; // <path_length>'s
368 
369  // loop over all xml tree nodes that are children of the root node
370  while (xmlCur != NULL) {
371 
372  // enter everytime you find a <path_length> tag
373  if( (!xmlStrcmp(xmlCur->name, (const xmlChar *) "path_length")) ) {
374 
375  xmlNodePtr xmlPlVal = xmlCur->xmlChildrenNode;
376 
377  string spdgc = utils::str::TrimSpaces(
378  XmlParserUtils::GetAttribute(xmlCur, "pdgc"));
379 
380  string spl = XmlParserUtils::TrimSpaces(
381  xmlNodeListGetString(xml_doc, xmlPlVal, 1));
382 
383  LOG("PathS", pDEBUG) << "pdgc = " << spdgc << " --> pl = " << spl;
384 
385  int pdgc = atoi( spdgc.c_str() );
386  double pl = atof( spl.c_str() );
387 
388  TParticlePDG * p = pdglib->Find(pdgc);
389  if(!p) {
390  LOG("PathS", pERROR)
391  << "No particle with pdgc " << pdgc
392  << " found. Will not load its path length";
393  } else
394  this->insert( map<int, double>::value_type(pdgc, pl) );
395 
396  xmlFree(xmlPlVal);
397  }
398  xmlCur = xmlCur->next;
399  } // [end of] loop over tags within root elements
400 
401  xmlFree(xmlCur);
402  return kXmlOK;
403 }
404 //___________________________________________________________________________
405 void PathSegmentList::SaveAsXml(string filename) const
406 {
407 //! Save path length list to XML file
408 
409  LOG("PathS", pINFO)
410  << "Saving PathSegmentList as XML in file: " << filename;
411 
412  ofstream outxml(filename.c_str());
413  if(!outxml.is_open()) {
414  LOG("PathS", pERROR) << "Couldn't create file = " << filename;
415  return;
416  }
417  outxml << "<?xml version=\"1.0\" encoding=\"ISO-8859-1\"?>";
418  outxml << endl << endl;
419  outxml << "<!-- generated by PathSegmentList::SaveAsXml() -->";
420  outxml << endl << endl;
421 
422  outxml << "<path_length_list>" << endl << endl;
423 
424  PathSegmentList::const_iterator pl_iter;
425 
426  for(pl_iter = this->begin(); pl_iter != this->end(); ++pl_iter) {
427 
428  int pdgc = pl_iter->first;
429  double pl = pl_iter->second; // path length
430 
431  outxml << " <path_length pdgc=\"" << pdgc << "\">"
432  << pl << "</path_length>" << endl;
433  }
434  outxml << "</path_length_list>";
435  outxml << endl;
436 
437  outxml.close();
438 
439 }
440 #endif
441 //___________________________________________________________________________
443 {
444  this->Copy(list);
445  return (*this);
446 }
447 //___________________________________________________________________________
PathSegmentV_t fSegmentList
Actual list of segments.
#define pERROR
Definition: Messenger.h:59
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
#define pFATAL
Definition: Messenger.h:56
TVector3 fStartPos
Record, for future comparison, the path taken.
Object to be filled with the neutrino path-segments representing geometry volume steps (generally bou...
string dir
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:96
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)
void SetStartInfo(const TVector3 &pos=TVector3(0, 0, 1e37), const TVector3 &dir=TVector3(0, 0, 0))
#define pINFO
Definition: Messenger.h:62
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
static constexpr double ps
Definition: Units.h:99
string TrimSpaces(string input)
Definition: StringUtils.cxx:18
MaterialMap_t::const_iterator MaterialMapCItr_t
Double_t fRayDist
distance from start of ray
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
Singleton class to load &amp; serve a TDatabasePDG.
Definition: PDGLibrary.h:35
ostream & operator<<(ostream &stream, const AlgConfigPool &config_pool)
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)
vector< vector< double > > clear
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
void DoCrossCheck(const TVector3 &startpos, double &ddist, double &dstep) const
perform cross check on segment, return differences
enum genie::EXmlParseStatus XmlParserStatus_t
static constexpr double sr
Definition: Units.h:166
string GetAttribute(xmlNodePtr xml_cur, string attr_name)
string Vec3AsString(const TVector3 *vec)
Definition: PrintUtils.cxx:80
Double_t fStepLength
total step size in volume
#define pDEBUG
Definition: Messenger.h:63
std::string fPathString
full path names