Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/persistency/gdml/src/G4STRead.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 // G4STRead implementation
 27 //
 28 // Author: Zoltan Torzsok, November 2007
 29 // --------------------------------------------------------------------
 30 
 31 #include <fstream>
 32 
 33 #include "G4STRead.hh"
 34 
 35 #include "G4Material.hh"
 36 #include "G4Box.hh"
 37 #include "G4QuadrangularFacet.hh"
 38 #include "G4TriangularFacet.hh"
 39 #include "G4TessellatedSolid.hh"
 40 #include "G4LogicalVolume.hh"
 41 #include "G4PVPlacement.hh"
 42 #include "G4AffineTransform.hh"
 43 #include "G4VoxelLimits.hh"
 44 
 45 // --------------------------------------------------------------------
 46 void G4STRead::TessellatedRead(const std::string& line)
 47 {
 48   if(tessellatedList.size() > 0)
 49   {
 50     tessellatedList.back()->SetSolidClosed(true);
 51     // Finish the previous solid at first!
 52   }
 53 
 54   std::istringstream stream(line.substr(2));
 55 
 56   G4String name;
 57   stream >> name;
 58 
 59   G4TessellatedSolid* tessellated = new G4TessellatedSolid(name);
 60   volumeMap[tessellated] =
 61     new G4LogicalVolume(tessellated, solid_material, name + "_LV", 0, 0, 0);
 62   tessellatedList.push_back(tessellated);
 63 
 64 #ifdef G4VERBOSE
 65   G4cout << "G4STRead: Reading solid: " << name << G4endl;
 66 #endif
 67 }
 68 
 69 // --------------------------------------------------------------------
 70 void G4STRead::FacetRead(const std::string& line)
 71 {
 72   if(tessellatedList.size() == 0)
 73   {
 74     G4Exception("G4STRead::FacetRead()", "ReadError", FatalException,
 75                 "A solid must be defined before defining a facet!");
 76   }
 77 
 78   if(line[2] == '3')  // Triangular facet
 79   {
 80     G4double x1, y1, z1;
 81     G4double x2, y2, z2;
 82     G4double x3, y3, z3;
 83 
 84     std::istringstream stream(line.substr(4));
 85     stream >> x1 >> y1 >> z1 >> x2 >> y2 >> z2 >> x3 >> y3 >> z3;
 86     tessellatedList.back()->AddFacet(new G4TriangularFacet(
 87       G4ThreeVector(x1, y1, z1), G4ThreeVector(x2, y2, z2),
 88       G4ThreeVector(x3, y3, z3), ABSOLUTE));
 89   }
 90   else if(line[2] == '4')  // Quadrangular facet
 91   {
 92     G4double x1, y1, z1;
 93     G4double x2, y2, z2;
 94     G4double x3, y3, z3;
 95     G4double x4, y4, z4;
 96 
 97     std::istringstream stream(line.substr(4));
 98     stream >> x1 >> y1 >> z1 >> x2 >> y2 >> z2 >> x3 >> y3 >> z3 >> x4 >> y4 >>
 99       z4;
100     tessellatedList.back()->AddFacet(new G4QuadrangularFacet(
101       G4ThreeVector(x1, y1, z1), G4ThreeVector(x2, y2, z2),
102       G4ThreeVector(x3, y3, z3), G4ThreeVector(x4, y4, z4), ABSOLUTE));
103   }
104   else
105   {
106     G4Exception("G4STRead::FacetRead()", "ReadError", FatalException,
107                 "Number of vertices per facet should be either 3 or 4!");
108   }
109 }
110 
111 // --------------------------------------------------------------------
112 void G4STRead::PhysvolRead(const std::string& line)
113 {
114   G4int level;
115   G4String name;
116   G4double r1, r2, r3;
117   G4double r4, r5, r6;
118   G4double r7, r8, r9;
119   G4double pX, pY, pZ;
120   G4double n1, n2, n3, n4, n5;
121 
122   std::istringstream stream(line.substr(2));
123   stream >> level >> name >> r1 >> r2 >> r3 >> n1 >> r4 >> r5 >> r6 >> n2 >>
124     r7 >> r8 >> r9 >> n3 >> pX >> pY >> pZ >> n4 >> n5;
125   std::string::size_type idx = name.rfind("_");
126   if(idx != std::string::npos)
127   {
128     name.resize(idx);
129   }
130   else
131   {
132     G4Exception("G4STRead::PhysvolRead()", "ReadError", FatalException,
133                 "Invalid input stream!");
134     return;
135   }
136 
137   G4cout << "G4STRead: Placing tessellated solid: " << name << G4endl;
138 
139   G4TessellatedSolid* tessellated = nullptr;
140 
141   for(std::size_t i = 0; i < tessellatedList.size(); ++i)
142   {  // Find the volume for this physvol!
143     if(tessellatedList[i]->GetName() == G4String(name))
144     {
145       tessellated = tessellatedList[i];
146       break;
147     }
148   }
149 
150   if(tessellated == nullptr)
151   {
152     G4String error_msg = "Referenced solid '" + name + "' not found!";
153     G4Exception("G4STRead::PhysvolRead()", "ReadError", FatalException,
154                 error_msg);
155   }
156   if(volumeMap.find(tessellated) == volumeMap.cend())
157   {
158     G4String error_msg = "Referenced solid '" + name +
159                          "' is not associated with a logical volume!";
160     G4Exception("G4STRead::PhysvolRead()", "InvalidSetup", FatalException,
161                 error_msg);
162   }
163   const G4RotationMatrix rot(G4ThreeVector(r1, r2, r3),
164                              G4ThreeVector(r4, r5, r6),
165                              G4ThreeVector(r7, r8, r9));
166   const G4ThreeVector pos(pX, pY, pZ);
167 
168   new G4PVPlacement(G4Transform3D(rot.inverse(), pos), volumeMap[tessellated],
169                     name + "_PV", world_volume, 0, 0);
170   // Note: INVERSE of rotation is needed!!!
171 
172   G4double minx, miny, minz;
173   G4double maxx, maxy, maxz;
174   const G4VoxelLimits limits;
175 
176   tessellated->CalculateExtent(kXAxis, limits, G4AffineTransform(rot, pos),
177                                minx, maxx);
178   tessellated->CalculateExtent(kYAxis, limits, G4AffineTransform(rot, pos),
179                                miny, maxy);
180   tessellated->CalculateExtent(kZAxis, limits, G4AffineTransform(rot, pos),
181                                minz, maxz);
182 
183   if(world_extent.x() < std::fabs(minx))
184   {
185     world_extent.setX(std::fabs(minx));
186   }
187   if(world_extent.y() < std::fabs(miny))
188   {
189     world_extent.setY(std::fabs(miny));
190   }
191   if(world_extent.z() < std::fabs(minz))
192   {
193     world_extent.setZ(std::fabs(minz));
194   }
195   if(world_extent.x() < std::fabs(maxx))
196   {
197     world_extent.setX(std::fabs(maxx));
198   }
199   if(world_extent.y() < std::fabs(maxy))
200   {
201     world_extent.setY(std::fabs(maxy));
202   }
203   if(world_extent.z() < std::fabs(maxz))
204   {
205     world_extent.setZ(std::fabs(maxz));
206   }
207 }
208 
209 // --------------------------------------------------------------------
210 void G4STRead::ReadGeom(const G4String& name)
211 {
212 #ifdef G4VERBOSE
213   G4cout << "G4STRead: Reading '" << name << "'..." << G4endl;
214 #endif
215   std::ifstream GeomFile(name);
216 
217   if(!GeomFile)
218   {
219     G4String error_msg = "Cannot open file: " + name;
220     G4Exception("G4STRead::ReadGeom()", "ReadError", FatalException, error_msg);
221   }
222 
223   tessellatedList.clear();
224   volumeMap.clear();
225   std::string line;
226 
227   while(getline(GeomFile, line))
228   {
229     if(line[0] == 'f')
230     {
231       TessellatedRead(line);
232     }
233     else if(line[0] == 'p')
234     {
235       FacetRead(line);
236     }
237   }
238 
239   if(tessellatedList.size() > 0)  // Finish the last solid!
240   {
241     tessellatedList.back()->SetSolidClosed(true);
242   }
243 
244   G4cout << "G4STRead: Reading '" << name << "' done." << G4endl;
245 }
246 
247 // --------------------------------------------------------------------
248 void G4STRead::ReadTree(const G4String& name)
249 {
250 #ifdef G4VERBOSE
251   G4cout << "G4STRead: Reading '" << name << "'..." << G4endl;
252 #endif
253   std::ifstream TreeFile(name);
254 
255   if(!TreeFile)
256   {
257     G4String error_msg = "Cannot open file: " + name;
258     G4Exception("G4STRead::ReadTree()", "ReadError", FatalException, error_msg);
259   }
260 
261   std::string line;
262 
263   while(getline(TreeFile, line))
264   {
265     if(line[0] == 'g')
266     {
267       PhysvolRead(line);
268     }
269   }
270 
271   G4cout << "G4STRead: Reading '" << name << "' done." << G4endl;
272 }
273 
274 // --------------------------------------------------------------------
275 G4LogicalVolume* G4STRead::Read(const G4String& name,
276                                 G4Material* mediumMaterial,
277                                 G4Material* solidMaterial)
278 {
279   if(mediumMaterial == nullptr)
280   {
281     G4Exception("G4STRead::Read()", "InvalidSetup", FatalException,
282                 "Pointer to medium material is not valid!");
283   }
284   if(solidMaterial == nullptr)
285   {
286     G4Exception("G4STRead::Read()", "InvalidSetup", FatalException,
287                 "Pointer to solid material is not valid!");
288   }
289 
290   solid_material = solidMaterial;
291 
292   world_box = new G4Box("TessellatedWorldBox", kInfinity, kInfinity, kInfinity);
293   // We don't know the extent of the world yet!
294   world_volume = new G4LogicalVolume(world_box, mediumMaterial,
295                                      "TessellatedWorldLV", 0, 0, 0);
296   world_extent = G4ThreeVector(0, 0, 0);
297 
298   ReadGeom(name + ".geom");
299   ReadTree(name + ".tree");
300 
301   // Now setting the world extent ...
302   //
303   if(world_box->GetXHalfLength() > world_extent.x())
304   {
305     world_box->SetXHalfLength(world_extent.x());
306   }
307   if(world_box->GetYHalfLength() > world_extent.y())
308   {
309     world_box->SetYHalfLength(world_extent.y());
310   }
311   if(world_box->GetZHalfLength() > world_extent.z())
312   {
313     world_box->SetZHalfLength(world_extent.z());
314   }
315 
316   return world_volume;
317 }
318