ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4GDMLWriteStructure.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4GDMLWriteStructure.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: PHG4GDMLWriteStructure.cc 68053 2013-03-13 14:39:51Z gcosmo $
28 //
29 // class PHG4GDMLWriteStructure Implementation
30 //
31 // Original author: Zoltan Torzsok, November 2007
32 //
33 // --------------------------------------------------------------------
34 
35 #include "PHG4GDMLWriteStructure.hh"
36 #include "PHG4GDMLConfig.hh"
37 
38 #include <Geant4/G4Material.hh>
39 #include <Geant4/G4ReflectedSolid.hh>
40 #include <Geant4/G4DisplacedSolid.hh>
41 #include <Geant4/G4LogicalVolumeStore.hh>
42 #include <Geant4/G4PhysicalVolumeStore.hh>
43 #include <Geant4/G4ReflectionFactory.hh>
44 #include <Geant4/G4PVDivision.hh>
45 #include <Geant4/G4PVReplica.hh>
46 #include <Geant4/G4Region.hh>
47 #include <Geant4/G4OpticalSurface.hh>
48 #include <Geant4/G4LogicalSkinSurface.hh>
49 #include <Geant4/G4LogicalBorderSurface.hh>
50 
51 #include <Geant4/G4ProductionCuts.hh>
52 #include <Geant4/G4ProductionCutsTable.hh>
53 #include <Geant4/G4Gamma.hh>
54 #include <Geant4/G4Electron.hh>
55 #include <Geant4/G4Positron.hh>
56 #include <Geant4/G4Proton.hh>
57 
58 
59 #include <cassert>
60 
61 
62 
63 PHG4GDMLWriteStructure::PHG4GDMLWriteStructure(const PHG4GDMLConfig * config_input)
64  : PHG4GDMLWriteParamvol(),structureElement(nullptr), cexport(false), config(config_input)
65 {
66  assert(config);
67  reflFactory = G4ReflectionFactory::Instance();
68 }
69 
70 PHG4GDMLWriteStructure::~PHG4GDMLWriteStructure()
71 {
72 }
73 
74 void
75 PHG4GDMLWriteStructure::DivisionvolWrite(xercesc::DOMElement* volumeElement,
76  const G4PVDivision* const divisionvol)
77 {
78  EAxis axis = kUndefined;
79  G4int number = 0;
80  G4double width = 0.0;
81  G4double offset = 0.0;
82  G4bool consuming = false;
83 
84  divisionvol->GetReplicationData(axis,number,width,offset,consuming);
85  axis = divisionvol->GetDivisionAxis();
86 
87  G4String unitString("mm");
88  G4String axisString("kUndefined");
89  if (axis==kXAxis) { axisString = "kXAxis"; }
90  else if (axis==kYAxis) { axisString = "kYAxis"; }
91  else if (axis==kZAxis) { axisString = "kZAxis"; }
92  else if (axis==kRho) { axisString = "kRho"; }
93  else if (axis==kPhi) { axisString = "kPhi"; unitString = "rad"; }
94 
95  const G4String name
96  = GenerateName(divisionvol->GetName(),divisionvol);
97  const G4String volumeref
98  = GenerateName(divisionvol->GetLogicalVolume()->GetName(),
99  divisionvol->GetLogicalVolume());
100 
101  xercesc::DOMElement* divisionvolElement = NewElement("divisionvol");
102  divisionvolElement->setAttributeNode(NewAttribute("axis",axisString));
103  divisionvolElement->setAttributeNode(NewAttribute("number",number));
104  divisionvolElement->setAttributeNode(NewAttribute("width",width));
105  divisionvolElement->setAttributeNode(NewAttribute("offset",offset));
106  divisionvolElement->setAttributeNode(NewAttribute("unit",unitString));
107  xercesc::DOMElement* volumerefElement = NewElement("volumeref");
108  volumerefElement->setAttributeNode(NewAttribute("ref",volumeref));
109  divisionvolElement->appendChild(volumerefElement);
110  volumeElement->appendChild(divisionvolElement);
111 }
112 
113 void PHG4GDMLWriteStructure::PhysvolWrite(xercesc::DOMElement* volumeElement,
114  const G4VPhysicalVolume* const physvol,
115  const G4Transform3D& T,
116  const G4String& ModuleName)
117 {
118  HepGeom::Scale3D scale;
119  HepGeom::Rotate3D rotate;
120  HepGeom::Translate3D translate;
121 
122  T.getDecomposition(scale,rotate,translate);
123 
124  const G4ThreeVector scl(scale(0,0),scale(1,1),scale(2,2));
125  const G4ThreeVector rot = GetAngles(rotate.getRotation());
126  const G4ThreeVector pos = T.getTranslation();
127 
128  const G4String name = GenerateName(physvol->GetName(),physvol);
129  const G4int copynumber = physvol->GetCopyNo();
130 
131  xercesc::DOMElement* physvolElement = NewElement("physvol");
132  physvolElement->setAttributeNode(NewAttribute("name",name));
133  if (copynumber) physvolElement->setAttributeNode(NewAttribute("copynumber", copynumber));
134 
135  volumeElement->appendChild(physvolElement);
136 
137  G4LogicalVolume* lv;
138  // Is it reflected?
139  if (reflFactory->IsReflected(physvol->GetLogicalVolume()))
140  {
141  lv = reflFactory->GetConstituentLV(physvol->GetLogicalVolume());
142  }
143  else
144  {
145  lv = physvol->GetLogicalVolume();
146  }
147 
148  const G4String volumeref = GenerateName(lv->GetName(), lv);
149 
150  if (ModuleName.empty())
151  {
152  xercesc::DOMElement* volumerefElement = NewElement("volumeref");
153  volumerefElement->setAttributeNode(NewAttribute("ref",volumeref));
154  physvolElement->appendChild(volumerefElement);
155  }
156  else
157  {
158  xercesc::DOMElement* fileElement = NewElement("file");
159  fileElement->setAttributeNode(NewAttribute("name",ModuleName));
160  fileElement->setAttributeNode(NewAttribute("volname",volumeref));
161  physvolElement->appendChild(fileElement);
162  }
163 
164  if (std::fabs(pos.x()) > kLinearPrecision
165  || std::fabs(pos.y()) > kLinearPrecision
166  || std::fabs(pos.z()) > kLinearPrecision)
167  {
168  PositionWrite(physvolElement,name+"_pos",pos);
169  }
170  if (std::fabs(rot.x()) > kAngularPrecision
171  || std::fabs(rot.y()) > kAngularPrecision
172  || std::fabs(rot.z()) > kAngularPrecision)
173  {
174  RotationWrite(physvolElement,name+"_rot",rot);
175  }
176  if (std::fabs(scl.x()-1.0) > kRelativePrecision
177  || std::fabs(scl.y()-1.0) > kRelativePrecision
178  || std::fabs(scl.z()-1.0) > kRelativePrecision)
179  {
180  ScaleWrite(physvolElement,name+"_scl",scl);
181  }
182 }
183 
184 void PHG4GDMLWriteStructure::ReplicavolWrite(xercesc::DOMElement* volumeElement,
185  const G4VPhysicalVolume* const replicavol)
186 {
187  EAxis axis = kUndefined;
188  G4int number = 0;
189  G4double width = 0.0;
190  G4double offset = 0.0;
191  G4bool consuming = false;
192  G4String unitString("mm");
193 
194  replicavol->GetReplicationData(axis,number,width,offset,consuming);
195 
196  const G4String volumeref
197  = GenerateName(replicavol->GetLogicalVolume()->GetName(),
198  replicavol->GetLogicalVolume());
199 
200  xercesc::DOMElement* replicavolElement = NewElement("replicavol");
201  replicavolElement->setAttributeNode(NewAttribute("number",number));
202  xercesc::DOMElement* volumerefElement = NewElement("volumeref");
203  volumerefElement->setAttributeNode(NewAttribute("ref",volumeref));
204  replicavolElement->appendChild(volumerefElement);
205  xercesc::DOMElement* replicateElement = NewElement("replicate_along_axis");
206  replicavolElement->appendChild(replicateElement);
207 
208  xercesc::DOMElement* dirElement = NewElement("direction");
209  if(axis==kXAxis)
210  { dirElement->setAttributeNode(NewAttribute("x","1")); }
211  else if(axis==kYAxis)
212  { dirElement->setAttributeNode(NewAttribute("y","1")); }
213  else if(axis==kZAxis)
214  { dirElement->setAttributeNode(NewAttribute("z","1")); }
215  else if(axis==kRho)
216  { dirElement->setAttributeNode(NewAttribute("rho","1")); }
217  else if(axis==kPhi)
218  { dirElement->setAttributeNode(NewAttribute("phi","1"));
219  unitString="rad"; }
220  replicateElement->appendChild(dirElement);
221 
222  xercesc::DOMElement* widthElement = NewElement("width");
223  widthElement->setAttributeNode(NewAttribute("value",width));
224  widthElement->setAttributeNode(NewAttribute("unit",unitString));
225  replicateElement->appendChild(widthElement);
226 
227  xercesc::DOMElement* offsetElement = NewElement("offset");
228  offsetElement->setAttributeNode(NewAttribute("value",offset));
229  offsetElement->setAttributeNode(NewAttribute("unit",unitString));
230  replicateElement->appendChild(offsetElement);
231 
232  volumeElement->appendChild(replicavolElement);
233 }
234 
235 void PHG4GDMLWriteStructure::
236 BorderSurfaceCache(const G4LogicalBorderSurface* const bsurf)
237 {
238  if (!bsurf) { return; }
239 
240  const G4SurfaceProperty* psurf = bsurf->GetSurfaceProperty();
241 
242  // Generate the new element for border-surface
243  //
244  xercesc::DOMElement* borderElement = NewElement("bordersurface");
245  borderElement->setAttributeNode(NewAttribute("name", bsurf->GetName()));
246  borderElement->setAttributeNode(NewAttribute("surfaceproperty",
247  psurf->GetName()));
248 
249  const G4String volumeref1 = GenerateName(bsurf->GetVolume1()->GetName(),
250  bsurf->GetVolume1());
251  const G4String volumeref2 = GenerateName(bsurf->GetVolume2()->GetName(),
252  bsurf->GetVolume2());
253  xercesc::DOMElement* volumerefElement1 = NewElement("physvolref");
254  xercesc::DOMElement* volumerefElement2 = NewElement("physvolref");
255  volumerefElement1->setAttributeNode(NewAttribute("ref",volumeref1));
256  volumerefElement2->setAttributeNode(NewAttribute("ref",volumeref2));
257  borderElement->appendChild(volumerefElement1);
258  borderElement->appendChild(volumerefElement2);
259 
260  if (FindOpticalSurface(psurf))
261  {
262  const G4OpticalSurface* opsurf =
263  dynamic_cast<const G4OpticalSurface*>(psurf);
264  if (!opsurf)
265  {
266  G4Exception("PHG4GDMLWriteStructure::BorderSurfaceCache()",
267  "InvalidSetup", FatalException, "No optical surface found!");
268  return;
269  }
270  OpticalSurfaceWrite(solidsElement, opsurf);
271  }
272 
273  borderElementVec.push_back(borderElement);
274 }
275 
276 void PHG4GDMLWriteStructure::
277 SkinSurfaceCache(const G4LogicalSkinSurface* const ssurf)
278 {
279  if (!ssurf) { return; }
280 
281  const G4SurfaceProperty* psurf = ssurf->GetSurfaceProperty();
282 
283  // Generate the new element for border-surface
284  //
285  xercesc::DOMElement* skinElement = NewElement("skinsurface");
286  skinElement->setAttributeNode(NewAttribute("name", ssurf->GetName()));
287  skinElement->setAttributeNode(NewAttribute("surfaceproperty",
288  psurf->GetName()));
289 
290  const G4String volumeref = GenerateName(ssurf->GetLogicalVolume()->GetName(),
291  ssurf->GetLogicalVolume());
292  xercesc::DOMElement* volumerefElement = NewElement("volumeref");
293  volumerefElement->setAttributeNode(NewAttribute("ref",volumeref));
294  skinElement->appendChild(volumerefElement);
295 
296  if (FindOpticalSurface(psurf))
297  {
298  const G4OpticalSurface* opsurf =
299  dynamic_cast<const G4OpticalSurface*>(psurf);
300  if (!opsurf)
301  {
302  G4Exception("PHG4GDMLWriteStructure::SkinSurfaceCache()",
303  "InvalidSetup", FatalException, "No optical surface found!");
304  return;
305  }
306  OpticalSurfaceWrite(solidsElement, opsurf);
307  }
308 
309  skinElementVec.push_back(skinElement);
310 }
311 
312 G4bool PHG4GDMLWriteStructure::FindOpticalSurface(const G4SurfaceProperty* psurf)
313 {
314  const G4OpticalSurface* osurf = dynamic_cast<const G4OpticalSurface*>(psurf);
315  std::vector<const G4OpticalSurface*>::const_iterator pos;
316  pos = std::find(opt_vec.begin(), opt_vec.end(), osurf);
317  if (pos != opt_vec.end()) { return false; } // item already created!
318 
319  opt_vec.push_back(osurf); // cache it for future reference
320  return true;
321 }
322 
323 const G4LogicalSkinSurface*
324 PHG4GDMLWriteStructure::GetSkinSurface(const G4LogicalVolume* const lvol)
325 {
326  G4LogicalSkinSurface* surf = 0;
327  G4int nsurf = G4LogicalSkinSurface::GetNumberOfSkinSurfaces();
328  if (nsurf)
329  {
330  const G4LogicalSkinSurfaceTable* stable =
331  G4LogicalSkinSurface::GetSurfaceTable();
332  std::vector<G4LogicalSkinSurface*>::const_iterator pos;
333  for (pos = stable->begin(); pos != stable->end(); ++pos)
334  {
335  if (lvol == (*pos)->GetLogicalVolume())
336  {
337  surf = *pos; break;
338  }
339  }
340  }
341  return surf;
342 }
343 
344 const G4LogicalBorderSurface*
345 PHG4GDMLWriteStructure::GetBorderSurface(const G4VPhysicalVolume* const pvol)
346 {
347  G4LogicalBorderSurface* surf = 0;
348  G4int nsurf = G4LogicalBorderSurface::GetNumberOfBorderSurfaces();
349  if (nsurf)
350  {
351  const G4LogicalBorderSurfaceTable* btable =
352  G4LogicalBorderSurface::GetSurfaceTable();
353  std::vector<G4LogicalBorderSurface*>::const_iterator pos;
354  for (pos = btable->begin(); pos != btable->end(); ++pos)
355  {
356  if (pvol == (*pos)->GetVolume1()) // just the first in the couple
357  { // is enough
358  surf = *pos; break;
359  }
360  }
361  }
362  return surf;
363 }
364 
365 void PHG4GDMLWriteStructure::SurfacesWrite()
366 {
367  std::cout << "PHG4GDML: Writing surfaces..." << std::endl;
368 
369  std::vector<xercesc::DOMElement*>::const_iterator pos;
370  for (pos = skinElementVec.begin(); pos != skinElementVec.end(); ++pos)
371  {
372  structureElement->appendChild(*pos);
373  }
374  for (pos = borderElementVec.begin(); pos != borderElementVec.end(); ++pos)
375  {
376  structureElement->appendChild(*pos);
377  }
378 }
379 
380 void PHG4GDMLWriteStructure::StructureWrite(xercesc::DOMElement* gdmlElement)
381 {
382  std::cout << "PHG4GDML: Writing structure..." << std::endl;
383 
384  structureElement = NewElement("structure");
385  gdmlElement->appendChild(structureElement);
386 }
387 
388 G4Transform3D PHG4GDMLWriteStructure::
389 TraverseVolumeTree(const G4LogicalVolume* const volumePtr, const G4int depth)
390 {
391  if (VolumeMap().find(volumePtr) != VolumeMap().end())
392  {
393  return VolumeMap()[volumePtr]; // Volume is already processed
394  }
395 
396  //jump over the exclusions
397  assert(config);
398  if (config->get_excluded_logical_vol().find(volumePtr) != config->get_excluded_logical_vol().end())
399  {
400  return G4Transform3D::Identity;
401  }
402 
403 
404  G4VSolid* solidPtr = volumePtr->GetSolid();
405  G4Transform3D R,invR;
406  G4int trans=0;
407 
408  std::map<const G4LogicalVolume*, PHG4GDMLAuxListType>::iterator auxiter;
409 
410  while (true) // Solve possible displacement/reflection
411  { // of the referenced solid!
412  if (trans>maxTransforms)
413  {
414  G4String ErrorMessage = "Referenced solid in volume '"
415  + volumePtr->GetName()
416  + "' was displaced/reflected too many times!";
417  G4Exception("PHG4GDMLWriteStructure::TraverseVolumeTree()",
418  "InvalidSetup", FatalException, ErrorMessage);
419  }
420 
421  if (G4ReflectedSolid* refl = dynamic_cast<G4ReflectedSolid*>(solidPtr))
422  {
423  R = R*refl->GetTransform3D();
424  solidPtr = refl->GetConstituentMovedSolid();
425  trans++;
426  continue;
427  }
428 
429  if (G4DisplacedSolid* disp = dynamic_cast<G4DisplacedSolid*>(solidPtr))
430  {
431  R = R*G4Transform3D(disp->GetObjectRotation(),
432  disp->GetObjectTranslation());
433  solidPtr = disp->GetConstituentMovedSolid();
434  trans++;
435  continue;
436  }
437 
438  break;
439  }
440 
441  //check if it is a reflected volume
442  G4LogicalVolume* tmplv = const_cast<G4LogicalVolume*>(volumePtr);
443 
444  if (reflFactory->IsReflected(tmplv))
445  {
446  tmplv = reflFactory->GetConstituentLV(tmplv);
447  if (VolumeMap().find(tmplv) != VolumeMap().end())
448  {
449  return R; // Volume is already processed
450  }
451  }
452 
453  // Only compute the inverse when necessary!
454  //
455  if (trans>0) { invR = R.inverse(); }
456 
457  const G4String name
458  = GenerateName(tmplv->GetName(), tmplv);
459  const G4String materialref
460  = GenerateName(volumePtr->GetMaterial()->GetName(),
461  volumePtr->GetMaterial());
462  const G4String solidref
463  = GenerateName(solidPtr->GetName(),solidPtr);
464 
465  xercesc::DOMElement* volumeElement = NewElement("volume");
466  volumeElement->setAttributeNode(NewAttribute("name",name));
467  xercesc::DOMElement* materialrefElement = NewElement("materialref");
468  materialrefElement->setAttributeNode(NewAttribute("ref",materialref));
469  volumeElement->appendChild(materialrefElement);
470  xercesc::DOMElement* solidrefElement = NewElement("solidref");
471  solidrefElement->setAttributeNode(NewAttribute("ref",solidref));
472  volumeElement->appendChild(solidrefElement);
473 
474  const G4int daughterCount = volumePtr->GetNoDaughters();
475 
476  for (G4int i=0;i<daughterCount;i++) // Traverse all the children!
477  {
478  const G4VPhysicalVolume* const physvol = volumePtr->GetDaughter(i);
479 
480 
481  //jump over the exclusions
482  assert(config);
483  if (config->get_excluded_physical_vol().find(physvol) != config->get_excluded_physical_vol().end())
484  continue;
485 
486 
487  const G4String ModuleName = Modularize(physvol,depth);
488 
489  G4Transform3D daughterR;
490 
491  if (ModuleName.empty()) // Check if subtree requested to be
492  { // a separate module!
493  daughterR = TraverseVolumeTree(physvol->GetLogicalVolume(),depth+1);
494  }
495  else
496  {
497  PHG4GDMLWriteStructure writer(config);
498  daughterR = writer.Write(ModuleName,physvol->GetLogicalVolume(),
499  SchemaLocation,depth+1);
500  }
501 
502  if (const G4PVDivision* const divisionvol
503  = dynamic_cast<const G4PVDivision*>(physvol)) // Is it division?
504  {
505  if (!G4Transform3D::Identity.isNear(invR*daughterR,kRelativePrecision))
506  {
507  G4String ErrorMessage = "Division volume in '" + name
508  + "' can not be related to reflected solid!";
509  G4Exception("PHG4GDMLWriteStructure::TraverseVolumeTree()",
510  "InvalidSetup", FatalException, ErrorMessage);
511  }
512  DivisionvolWrite(volumeElement,divisionvol);
513  } else
514  if (physvol->IsParameterised()) // Is it a paramvol?
515  {
516  if (!G4Transform3D::Identity.isNear(invR*daughterR,kRelativePrecision))
517  {
518  G4String ErrorMessage = "Parameterised volume in '" + name
519  + "' can not be related to reflected solid!";
520  G4Exception("PHG4GDMLWriteStructure::TraverseVolumeTree()",
521  "InvalidSetup", FatalException, ErrorMessage);
522  }
523  ParamvolWrite(volumeElement,physvol);
524  } else
525  if (physvol->IsReplicated()) // Is it a replicavol?
526  {
527  if (!G4Transform3D::Identity.isNear(invR*daughterR,kRelativePrecision))
528  {
529  G4String ErrorMessage = "Replica volume in '" + name
530  + "' can not be related to reflected solid!";
531  G4Exception("PHG4GDMLWriteStructure::TraverseVolumeTree()",
532  "InvalidSetup", FatalException, ErrorMessage);
533  }
534  ReplicavolWrite(volumeElement,physvol);
535  }
536  else // Is it a physvol?
537  {
538  G4RotationMatrix rot;
539  if (physvol->GetFrameRotation() != 0)
540  {
541  rot = *(physvol->GetFrameRotation());
542  }
543  G4Transform3D P(rot,physvol->GetObjectTranslation());
544 
545  PhysvolWrite(volumeElement,physvol,invR*P*daughterR,ModuleName);
546  }
547  BorderSurfaceCache(GetBorderSurface(physvol));
548  }
549 
550  if (cexport) { ExportEnergyCuts(volumePtr); }
551  // Add optional energy cuts
552 
553  // Here write the auxiliary info
554  //
555  auxiter = auxmap.find(volumePtr);
556  if (auxiter != auxmap.end())
557  {
558  AddAuxInfo(&(auxiter->second), volumeElement);
559  }
560 
561  structureElement->appendChild(volumeElement);
562  // Append the volume AFTER traversing the children so that
563  // the order of volumes will be correct!
564 
565  VolumeMap()[tmplv] = R;
566 
567  AddExtension(volumeElement, volumePtr);
568  // Add any possible user defined extension attached to a volume
569 
570  AddMaterial(volumePtr->GetMaterial());
571  // Add the involved materials and solids!
572 
573  AddSolid(solidPtr);
574 
575  SkinSurfaceCache(GetSkinSurface(volumePtr));
576 
577  return R;
578 }
579 
580 void
581 PHG4GDMLWriteStructure::AddVolumeAuxiliary(PHG4GDMLAuxStructType myaux,
582  const G4LogicalVolume* const lvol)
583 {
584  std::map<const G4LogicalVolume*,
585  PHG4GDMLAuxListType>::iterator pos = auxmap.find(lvol);
586 
587  if (pos == auxmap.end()) { auxmap[lvol] = PHG4GDMLAuxListType(); }
588 
589  auxmap[lvol].push_back(myaux);
590 }
591 
592 void
593 PHG4GDMLWriteStructure::SetEnergyCutsExport(G4bool fcuts)
594 {
595  cexport = fcuts;
596 }
597 
598 void
599 PHG4GDMLWriteStructure::ExportEnergyCuts(const G4LogicalVolume* const lvol)
600 {
601 // PHG4GDMLEvaluator eval;
602  G4ProductionCuts* pcuts = lvol->GetRegion()->GetProductionCuts();
603  G4ProductionCutsTable* ctab = G4ProductionCutsTable::GetProductionCutsTable();
604  G4Gamma* gamma = G4Gamma::Gamma();
605  G4Electron* eminus = G4Electron::Electron();
606  G4Positron* eplus = G4Positron::Positron();
607  G4Proton* proton = G4Proton::Proton();
608 
609  G4double gamma_cut = ctab->ConvertRangeToEnergy(gamma, lvol->GetMaterial(),
610  pcuts->GetProductionCut("gamma"));
611  G4double eminus_cut = ctab->ConvertRangeToEnergy(eminus, lvol->GetMaterial(),
612  pcuts->GetProductionCut("e-"));
613  G4double eplus_cut = ctab->ConvertRangeToEnergy(eplus, lvol->GetMaterial(),
614  pcuts->GetProductionCut("e+"));
615  G4double proton_cut = ctab->ConvertRangeToEnergy(proton, lvol->GetMaterial(),
616  pcuts->GetProductionCut("proton"));
617 
618  PHG4GDMLAuxStructType gammainfo = {"gammaECut",
619  ConvertToString(gamma_cut), "MeV", 0};
620  PHG4GDMLAuxStructType eminusinfo = {"electronECut",
621  ConvertToString(eminus_cut), "MeV", 0};
622  PHG4GDMLAuxStructType eplusinfo = {"positronECut",
623  ConvertToString(eplus_cut), "MeV", 0};
624  PHG4GDMLAuxStructType protinfo = {"protonECut",
625  ConvertToString(proton_cut), "MeV", 0};
626 
627  AddVolumeAuxiliary(gammainfo, lvol);
628  AddVolumeAuxiliary(eminusinfo, lvol);
629  AddVolumeAuxiliary(eplusinfo, lvol);
630  AddVolumeAuxiliary(protinfo, lvol);
631 }
632 
633 G4String PHG4GDMLWriteStructure::ConvertToString(G4double dval)
634 {
635  std::ostringstream os;
636  os << dval;
637  G4String vl = os.str();
638  return vl;
639 }
640