Bug Summary

File:scratch/gluex/scan-build-work/hdgeant4/src/HddsG4Builder.cc
Location:line 733, column 7
Description:Value stored to 'playerN' is never read

Annotated Source Code

1//
2// HddsG4Builder - class implementation
3//
4// author: richard.t.jones at uconn.edu
5// version: may 12, 2012
6//
7
8#define APP_NAME"HddsG4Builder" "HddsG4Builder"
9
10#include <HddsG4Builder.hh>
11#include <G4Box.hh>
12#include <G4Trd.hh>
13#include <G4Trap.hh>
14#include <G4Para.hh>
15#include <G4Tubs.hh>
16#include <G4Cons.hh>
17#include <G4Hype.hh>
18#include <G4Torus.hh>
19#include <G4Orb.hh>
20#include <G4Sphere.hh>
21#include <G4Polycone.hh>
22#include <G4Polyhedra.hh>
23#include <G4EllipticalTube.hh>
24#include <G4SystemOfUnits.hh>
25#include <G4Mag_UsualEqRhs.hh>
26#include <G4MagIntegratorStepper.hh>
27#include <G4ExactHelixStepper.hh>
28#include <G4HelixMixedStepper.hh>
29#include <G4ClassicalRK4.hh>
30#include <G4ChordFinder.hh>
31#include <G4PVPlacement.hh>
32#include <G4PVDivision.hh>
33#include <G4OpticalSurface.hh>
34#include <G4LogicalSkinSurface.hh>
35#include <G4VisAttributes.hh>
36
37#include <xercesc/util/PlatformUtils.hpp>
38#include <xercesc/util/XMLString.hpp>
39#include <xercesc/util/XMLStringTokenizer.hpp>
40#include <xercesc/sax/SAXParseException.hpp>
41#include <xercesc/parsers/XercesDOMParser.hpp>
42#include <xercesc/framework/LocalFileFormatTarget.hpp>
43#include <xercesc/dom/DOM.hpp>
44#include <xercesc/util/XercesDefs.hpp>
45#include <xercesc/sax/ErrorHandler.hpp>
46
47using namespace xercesc;
48
49#include <assert.h>
50#include <stdlib.h>
51#include <ctype.h>
52#include <stdio.h>
53#include <math.h>
54
55#include <iostream>
56#include <fstream>
57#include <sstream>
58#include <string>
59#include <iomanip>
60#include <vector>
61#include <list>
62#include <map>
63
64#define X(str)XString(str).unicode_str() XString(str).unicode_str()
65#define S(str)str.c_str() str.c_str()
66
67#ifdef LINUX_CPUTIME_PROFILING
68extern CPUtimer timer;
69#endif
70
71HddsG4Builder::HddsG4Builder() : fWorldVolume(0) { }
72
73HddsG4Builder::HddsG4Builder(const HddsG4Builder &src)
74{
75 fWorldVolume = src.fWorldVolume;
76 fElements = src.fElements;
77 fMaterials = src.fMaterials;
78 fMagneticRegions = src.fMagneticRegions;
79 fLogicalVolumes = src.fLogicalVolumes;
80 fPhysicalVolumes = src.fPhysicalVolumes;
81 fSensitiveVolumes = src.fSensitiveVolumes;
82 fRotations = src.fRotations;
83 fCurrentMother = src.fCurrentMother;
84 fCurrentPlacement = src.fCurrentPlacement;
85 fCurrentDivision = src.fCurrentDivision;
86 fCurrentPhiCenter = src.fCurrentPhiCenter;
87}
88
89HddsG4Builder::~HddsG4Builder() { }
90
91int HddsG4Builder::createMaterial(DOMElement* el)
92{
93#ifdef LINUX_CPUTIME_PROFILING
94 std::stringstream timestr;
95 timestr << "createMaterial: " << timer.getUserTime() << ", "
96 << timer.getSystemTime() << ", " << timer.getRealTime();
97 timer.resetClocks();
98#endif
99 int imate = CodeWriter::createMaterial(el);
100
101 if (fSubst.fBrewList.size() == 0)
102 {
103 XString matS = fSubst.getName();
104 XString symS = fSubst.getSymbol();
105 double A = fSubst.getAtomicWeight();
106 double Z = fSubst.getAtomicNumber();
107 double dens = fSubst.getDensity();
108 fElements[imate] = new G4Element(matS,symS,Z,A*g/mole);
109 fMaterials[imate] = new G4Material(matS,Z,A*g/mole,dens*g/cm3);
110 }
111 else
112 {
113 XString matS = fSubst.getName();
114 double dens = fSubst.getDensity();
115 int ncomp = fSubst.fBrewList.size();
116 fMaterials[imate] = new G4Material(matS,dens*g/cm3,ncomp);
117 std::list<Substance::Brew>::iterator iter;
118 for (iter = fSubst.fBrewList.begin();
119 iter != fSubst.fBrewList.end(); iter++)
120 {
121 int subimate = iter->sub->fUniqueID;
122 std::map<int,G4Element*>::iterator elemit = fElements.find(subimate);
123 if (iter->natoms)
124 {
125 fMaterials[imate]->AddElement(elemit->second,iter->natoms);
126 }
127 else if (elemit != fElements.end())
128 {
129 fMaterials[imate]->AddElement(elemit->second,iter->wfact);
130 }
131 else
132 {
133 fMaterials[imate]->AddMaterial(fMaterials[subimate],iter->wfact);
134 }
135 }
136 }
137
138 DOMNodeList* propList = el->getElementsByTagName(X("optical_properties")XString("optical_properties").unicode_str());
139 if (propList->getLength() > 0)
140 {
141 DOMElement* propEl = (DOMElement*)propList->item(0);
142 DOMNodeList* specL = propEl->getElementsByTagName(X("specify")XString("specify").unicode_str());
143 int len = specL->getLength();
144 std::vector<double> Ephot;
145 std::vector<double> abslen;
146 std::vector<double> effic;
147 std::vector<double> rindex;
148 std::vector<double> smooth;
149 std::vector<double> reflect;
150 for (int i=0; i < len; ++i)
151 {
152 DOMElement* specEl = (DOMElement*)specL->item(i);
153 XString valS;
154 valS = specEl->getAttribute(X("E")XString("E").unicode_str());
155 Ephot.push_back(atof(S(valS)valS.c_str())*eV);
156 valS = specEl->getAttribute(X("refindex")XString("refindex").unicode_str());
157 rindex.push_back(atof(S(valS)valS.c_str()));
158 valS = specEl->getAttribute(X("abslen")XString("abslen").unicode_str());
159 abslen.push_back(atof(S(valS)valS.c_str())*cm);
160 valS = specEl->getAttribute(X("smooth")XString("smooth").unicode_str());
161 smooth.push_back(atof(S(valS)valS.c_str()));
162 valS = specEl->getAttribute(X("reflect")XString("reflect").unicode_str());
163 reflect.push_back(atof(S(valS)valS.c_str()));
164 valS = specEl->getAttribute(X("effic")XString("effic").unicode_str());
165 effic.push_back(atof(S(valS)valS.c_str()));
166 }
167 G4MaterialPropertiesTable *mpt = new G4MaterialPropertiesTable();
168 if (rindex[0] > 0)
169 {
170 mpt->AddProperty("RINDEX",&Ephot[0],&rindex[0],len);
171 }
172 if (abslen[0] > 0)
173 {
174 mpt->AddProperty("ABSLENGTH",&Ephot[0],&abslen[0],len);
175 }
176 if (smooth[0] > 0)
177 {
178 mpt->AddProperty("POLISH",&Ephot[0],&smooth[0],len);
179 }
180 if (reflect[0] > 0)
181 {
182 mpt->AddProperty("REFLECTIVITY",&Ephot[0],&reflect[0],len);
183 }
184 if (effic[0] > 0)
185 {
186 mpt->AddProperty("EFFICIENCY",&Ephot[0],&effic[0],len);
187 }
188 fMaterials[imate]->SetMaterialPropertiesTable(mpt);
189 }
190#ifdef LINUX_CPUTIME_PROFILING
191 timestr << " ( " << timer.getUserDelta() << " ) ";
192 G4cerr(*G4cerr_p) << timestr.str() << G4endlstd::endl;
193#endif
194 return imate;
195}
196
197int HddsG4Builder::createSolid(DOMElement* el, Refsys& ref)
198{
199#ifdef LINUX_CPUTIME_PROFILING
200 std::stringstream timestr;
201 timestr << "createSolid: " << timer.getUserTime() << ", "
202 << timer.getSystemTime() << ", " << timer.getRealTime();
203 timer.resetClocks();
204#endif
205
206 int ivolu = CodeWriter::createSolid(el,ref);
207 XString nameS(el->getAttribute(X("name")XString("name").unicode_str()));
208 int imate = fSubst.fUniqueID;
209
210 Units unit;
211 unit.getConversions(el);
212
213 G4VSolid *solid;
214 XString shapeS(el->getTagName());
215 if (shapeS == "box")
216 {
217 double xl, yl, zl;
218 XString xyzS(el->getAttribute(X("X_Y_Z")XString("X_Y_Z").unicode_str()));
219 std::stringstream listr(xyzS);
220 listr >> xl >> yl >> zl;
221 solid = new G4Box(S(nameS)nameS.c_str(),
222 xl/2 *cm/unit.cm,
223 yl/2 *cm/unit.cm,
224 zl/2 *cm/unit.cm
225 );
226 }
227 else if (shapeS == "tubs")
228 {
229 double ri, ro, zl, phi0, dphi;
230 XString riozS(el->getAttribute(X("Rio_Z")XString("Rio_Z").unicode_str()));
231 std::stringstream listr(riozS);
232 listr >> ri >> ro >> zl;
233 XString profS(el->getAttribute(X("profile")XString("profile").unicode_str()));
234 listr.clear(), listr.str(profS);
235 listr >> phi0 >> dphi;
236 solid = new G4Tubs(S(nameS)nameS.c_str(),
237 ri * cm/unit.cm,
238 ro * cm/unit.cm,
239 zl/2 * cm/unit.cm,
240 phi0 * deg/unit.deg,
241 dphi * deg/unit.deg
242 );
243 }
244 else if (shapeS == "eltu")
245 {
246
247 double rx, ry, zl;
248 XString rxyzS(el->getAttribute(X("Rxy_Z")XString("Rxy_Z").unicode_str()));
249 std::stringstream listr(rxyzS);
250 listr >> rx >> ry >> zl;
251 solid = new G4EllipticalTube(S(nameS)nameS.c_str(),
252 rx * cm/unit.cm,
253 ry * cm/unit.cm,
254 zl/2 * cm/unit.cm
255 );
256 }
257 else if (shapeS == "trd")
258 {
259 double xm, ym, xp, yp, zl;
260 XString xyzS(el->getAttribute(X("Xmp_Ymp_Z")XString("Xmp_Ymp_Z").unicode_str()));
261 std::stringstream listr(xyzS);
262 listr >> xm >> xp >> ym >> yp >> zl;
263 xm += 1e-100;
264 xp += 1e-100;
265 ym += 1e-100;
266 yp += 1e-100;
267 zl += 1e-100;
268 double alph_xz, alph_yz;
269 XString incS(el->getAttribute(X("inclination")XString("inclination").unicode_str()));
270 listr.clear(), listr.str(incS);
271 listr >> alph_xz >> alph_yz;
272 double x = tan(alph_xz/unit.rad);
273 double y = tan(alph_yz/unit.rad);
274 double r = sqrt(x*x + y*y);
275 solid = new G4Trap(S(nameS)nameS.c_str(),
276 zl/2 * cm/unit.cm,
277 atan2(r,1) * rad,
278 atan2(y,x) * rad,
279 ym/2 * cm/unit.cm,
280 xm/2 * cm/unit.cm,
281 xm/2 * cm/unit.cm,
282 0,
283 yp/2 * cm/unit.cm,
284 xp/2 * cm/unit.cm,
285 xp/2 * cm/unit.cm,
286 0
287 );
288 }
289 else if (shapeS == "pcon")
290 {
291 double phi0, dphi;
292 XString profS(el->getAttribute(X("profile")XString("profile").unicode_str()));
293 std::stringstream listr(profS);
294 listr >> phi0 >> dphi;
295 DOMNodeList* planeList = el->getElementsByTagName(X("polyplane")XString("polyplane").unicode_str());
296 int nplanes = planeList->getLength();
297 std::vector<double> zPlane;
298 std::vector<double> rInner;
299 std::vector<double> rOuter;
300 double zlast = -1e30;
301 double zeps = 0;
302 for (unsigned int p = 0; p < planeList->getLength(); p++)
303 {
304 double ri, ro, zl;
305 DOMNode* node = planeList->item(p);
306 DOMElement* elem = (DOMElement*) node;
307 XString riozS(elem->getAttribute(X("Rio_Z")XString("Rio_Z").unicode_str()));
308 std::stringstream listr1(riozS);
309 listr1 >> ri >> ro >> zl;
310 if (zl < zlast)
311 {
312 G4cerr(*G4cerr_p) << APP_NAME"HddsG4Builder" << " error: Please re-order the polyplanes"
313 << " of volume " << S(nameS)nameS.c_str() << " so that the z-values"
314 << " are non-decreasing"
315 << G4endlstd::endl;
316 exit(1);
317 }
318 zl = (zl < zlast + zeps)? zlast + zeps : zl;
319 zlast = zl;
320 zPlane.push_back(zl * cm/unit.cm);
321 rInner.push_back(ri * cm/unit.cm);
322 rOuter.push_back(ro * cm/unit.cm);
323 }
324 solid = new G4Polycone(S(nameS)nameS.c_str(),
325 phi0 * deg/unit.deg,
326 dphi * deg/unit.deg,
327 nplanes,
328 &zPlane[0],
329 &rInner[0],
330 &rOuter[0]
331 );
332 }
333 else if (shapeS == "pgon")
334 {
335 int segments;
336 XString segS(el->getAttribute(X("segments")XString("segments").unicode_str()));
337 segments = atoi(S(segS)segS.c_str());
338 double phi0, dphi;
339 XString profS(el->getAttribute(X("profile")XString("profile").unicode_str()));
340 std::stringstream listr(profS);
341 listr >> phi0 >> dphi;
342 DOMNodeList* planeList = el->getElementsByTagName(X("polyplane")XString("polyplane").unicode_str());
343 int nplanes = planeList->getLength();
344 std::vector<double> zPlane;
345 std::vector<double> rInner;
346 std::vector<double> rOuter;
347 double zlast = -1e30;
348 double zeps = 0;
349 for (unsigned int p = 0; p < planeList->getLength(); p++)
350 {
351 double ri, ro, zl;
352 DOMNode* node = planeList->item(p);
353 DOMElement* elem = (DOMElement*) node;
354 XString riozS(elem->getAttribute(X("Rio_Z")XString("Rio_Z").unicode_str()));
355 std::stringstream listr1(riozS);
356 listr1 >> ri >> ro >> zl;
357 if (zl < zlast)
358 {
359 G4cerr(*G4cerr_p) << APP_NAME"HddsG4Builder" << " error: Please re-order the polyplanes"
360 << " of volume " << S(nameS)nameS.c_str() << " so that the z-values"
361 << " are non-decreasing."
362 << G4endlstd::endl;
363 exit(1);
364 }
365 zl = (zl < zlast + zeps)? zlast + zeps : zl;
366 zlast = zl;
367 zPlane.push_back(zl * cm/unit.cm);
368 rInner.push_back(ri * cm/unit.cm);
369 rOuter.push_back(ro * cm/unit.cm);
370 }
371 solid = new G4Polyhedra(S(nameS)nameS.c_str(),
372 phi0 * deg/unit.deg,
373 dphi * deg/unit.deg,
374 segments,
375 nplanes,
376 &zPlane[0],
377 &rInner[0],
378 &rOuter[0]
379 );
380 }
381 else if (shapeS == "cons")
382 {
383 double rim, rip, rom, rop, zl;
384 XString riozS(el->getAttribute(X("Rio1_Rio2_Z")XString("Rio1_Rio2_Z").unicode_str()));
385 std::stringstream listr(riozS);
386 listr >> rim >> rom >> rip >> rop >> zl;
387 double phi0, dphi;
388 XString profS(el->getAttribute(X("profile")XString("profile").unicode_str()));
389 listr.clear(), listr.str(profS);
390 listr >> phi0 >> dphi;
391 solid = new G4Cons(S(nameS)nameS.c_str(),
392 zl/2 * cm/unit.cm,
393 rim * cm/unit.cm,
394 rom * cm/unit.cm,
395 rip * cm/unit.cm,
396 rop * cm/unit.cm,
397 phi0 * deg/unit.deg,
398 dphi * deg/unit.deg
399 );
400 }
401 else if (shapeS == "sphere")
402 {
403 double ri, ro;
404 XString rioS(el->getAttribute(X("Rio")XString("Rio").unicode_str()));
405 std::stringstream listr(rioS);
406 listr >> ri >> ro;
407 double theta0, theta1;
408 XString polarS(el->getAttribute(X("polar_bounds")XString("polar_bounds").unicode_str()));
409 listr.clear(), listr.str(polarS);
410 listr >> theta0 >> theta1;
411 double phi0, dphi;
412 XString profS(el->getAttribute(X("profile")XString("profile").unicode_str()));
413 listr.clear(), listr.str(profS);
414 listr >> phi0 >> dphi;
415 solid = new G4Sphere(S(nameS)nameS.c_str(),
416 ri * cm/unit.cm,
417 ro * cm/unit.cm,
418 phi0 * deg/unit.deg,
419 dphi * deg/unit.deg,
420 theta0 * deg/unit.deg,
421 theta1 * deg/unit.deg
422 );
423 }
424 else
425 {
426 G4cerr(*G4cerr_p) << APP_NAME"HddsG4Builder" << " error: volume " << S(nameS)nameS.c_str()
427 << " should be one of the valid shapes, not " << S(shapeS)shapeS.c_str()
428 << G4endlstd::endl;
429 exit(1);
430 }
431
432 vpair_t newvol(ivolu,0);
433 G4FieldManager* fieldmgr = fFieldManagers[ref.fRegionID];
434 fLogicalVolumes[newvol] = new G4LogicalVolume(solid,fMaterials[imate],
435 S(nameS)nameS.c_str(),fieldmgr);
436 G4Material* mate = fMaterials[imate];
437 double dens = mate->GetDensity();
438 double dmod = (int)(dens*97.345/(g/cm3)) % 20;
439 double red = 1-3/(3+dens /(g/cm3));
440 double green = 1-5/(5+dmod);
441 double blue = 1-green;
442 G4VisAttributes* attr = new G4VisAttributes(G4Colour(red,green,blue));
443 fLogicalVolumes[newvol]->SetVisAttributes(attr);
444
445 XString sensiS(el->getAttribute(X("sensitive")XString("sensitive").unicode_str()));
446 if (sensiS == "true")
447 {
448 fSensitiveVolumes[ivolu] = fLogicalVolumes[newvol];
449 }
450
451 G4MaterialPropertiesTable *mpt = fMaterials[imate]->GetMaterialPropertiesTable();
452 if (mpt)
453 {
454 G4MaterialPropertyVector *poli_vector = mpt->GetProperty("POLISH");
455 G4MaterialPropertyVector *refl_vector = mpt->GetProperty("REFLECTIVITY");
456 G4MaterialPropertyVector *refi_vector = mpt->GetProperty("RINDEX");
457 G4MaterialPropertyVector *absl_vector = mpt->GetProperty("ABSLENGTH");
458 G4MaterialPropertyVector *effi_vector = mpt->GetProperty("EFFICIENCY");
459 if (refi_vector == 0) {
460 if (refl_vector == 0) {
461 G4cerr(*G4cerr_p) << APP_NAME"HddsG4Builder" << " error: material " << mate->GetName()
462 << " needs to have either refindex defined (dielectrics)"
463 << " or else reflect (metals) but it has neither!"
464 << G4endlstd::endl;
465 exit(1);
466 }
467 if (absl_vector != 0) {
468 G4cerr(*G4cerr_p) << APP_NAME"HddsG4Builder" << " error: material " << mate->GetName()
469 << " has abslen defined but no refindex; defining abslen"
470 << " requires that refindex also be specified!"
471 << G4endlstd::endl;
472 exit(1);
473 }
474 if (poli_vector != 0 || effi_vector != 0) {
475 G4OpticalSurface *surface = new G4OpticalSurface(S(nameS)nameS.c_str());
476 surface->SetType(dielectric_metal);
477 surface->SetModel(glisur);
478 surface->SetMaterialPropertiesTable(mpt);
479 if (poli_vector != 0) {
480 double polish = poli_vector->GetMaxValue();
481 surface->SetPolish(polish);
482 surface->SetFinish(ground);
483 }
484 else {
485 surface->SetFinish(polished);
486 }
487 new G4LogicalSkinSurface(S(nameS)nameS.c_str(),fLogicalVolumes[newvol],surface);
488 }
489 }
490 else if (refl_vector == 0) {
491 if (poli_vector != 0) {
492 G4OpticalSurface *surface = new G4OpticalSurface(S(nameS)nameS.c_str());
493 surface->SetType(dielectric_dielectric);
494 surface->SetModel(glisur);
495 surface->SetMaterialPropertiesTable(mpt);
496 if (poli_vector != 0) {
497 double polish = poli_vector->GetMaxValue();
498 surface->SetPolish(polish);
499 surface->SetFinish(ground);
500 }
501 else {
502 surface->SetFinish(polished);
503 }
504 new G4LogicalSkinSurface(S(nameS)nameS.c_str(),fLogicalVolumes[newvol],surface);
505 }
506 }
507 else {
508 G4cerr(*G4cerr_p) << APP_NAME"HddsG4Builder" << " error: material " << mate->GetName()
509 << " has optical properties, but neither refindex nor reflect"
510 << " data are provided, you need either one or the other!"
511 << G4endlstd::endl;
512 exit(1);
513 }
514 }
515
516#ifdef LINUX_CPUTIME_PROFILING
517 timestr << " ( " << timer.getUserDelta() << " ) ";
518 G4cerr(*G4cerr_p) << timestr.str() << G4endlstd::endl;
519#endif
520 return ivolu;
521}
522
523int HddsG4Builder::createRotation(Refsys& ref)
524{
525#ifdef LINUX_CPUTIME_PROFILING
526 std::stringstream timestr;
527 timestr << "createRotation: " << timer.getUserTime() << ", "
528 << timer.getSystemTime() << ", " << timer.getRealTime();
529 timer.resetClocks();
530#endif
531 int irot = CodeWriter::createRotation(ref);
532
533 if (irot > 0)
534 {
535 std::vector<double> omega = ref.getRotation();
536 fRotations[irot] = new G4RotationMatrix();
537 fRotations[irot]->rotateX(omega[0]);
538 fRotations[irot]->rotateY(omega[1]);
539 fRotations[irot]->rotateZ(omega[2]);
540 fRotations[irot]->invert();
541 }
542 else if (fRotations.count(0) == 0)
543 {
544 fRotations[0] = new G4RotationMatrix();
545 }
546
547#ifdef LINUX_CPUTIME_PROFILING
548 timestr << " ( " << timer.getUserDelta() << " ) ";
549 G4cerr(*G4cerr_p) << timestr.str() << G4endlstd::endl;
550#endif
551 return irot;
552}
553
554int HddsG4Builder::createVolume(DOMElement* el, Refsys& ref)
555{
556#ifdef LINUX_CPUTIME_PROFILING
557 std::stringstream timestr;
558 timestr << "createVolume: " << timer.getUserTime() << ", "
559 << timer.getSystemTime() << ", " << timer.getRealTime();
560 timer.resetClocks();
561#endif
562 int icopy = CodeWriter::createVolume(el,ref);
563
564 if (fPending)
565 {
566 XString myvoluS(el->getAttribute(X("HDDSvolu")XString("HDDSvolu").unicode_str()));
567 XString motherS(fRef.fMother->getAttribute(X("HDDSvolu")XString("HDDSvolu").unicode_str()));
568 int myvoluI = atoi(S(myvoluS)myvoluS.c_str());
569 int motherI = atoi(S(motherS)motherS.c_str());
570 G4ThreeVector origin(fRef.fOrigin[0]*cm,
571 fRef.fOrigin[1]*cm,
572 fRef.fOrigin[2]*cm);
573 int irot = fRef.fRotation;
574
575 // Apply fix-up in case we are placing the volume inside a phi division
576 // because of a flaw in the way geant4 handles this special case.
577
578 G4RotationMatrix *rot = fRotations[irot];
579 if (fCurrentPhiCenter[motherI] != 0)
580 {
581 rot = new G4RotationMatrix(*rot);
582 rot->rotateZ(-fCurrentPhiCenter[motherI]);
583 origin.rotateZ(fCurrentPhiCenter[motherI]);
584 }
585
586 std::map<vpair_t,G4LogicalVolume*>::iterator mine;
587 for (mine = fLogicalVolumes.find(vpair_t(myvoluI,0));
588 mine != fLogicalVolumes.end() && mine->first.first == myvoluI;
589 ++mine)
590 {
591 int ilayer = mine->first.second;
592 std::map<vpair_t,G4LogicalVolume*>::iterator moms;
593 moms = addNewLayer(motherI,fRef.fRelativeLayer + ilayer);
594 G4PVPlacement *pvol = new G4PVPlacement(rot, origin,
595 mine->second,
596 mine->second->GetName(),
597 moms->second,0,icopy);
598#ifdef CHECK_OVERLAPS_MM
599 pvol->CheckOverlaps(1000,CHECK_OVERLAPS_MM);
600#endif
601#ifdef DEBUG_PLACEMENT
602 XString nameS(el->getAttribute(X("name")XString("name").unicode_str()));
603 G4cout(*G4cout_p) << "volume " << nameS
604 << "->" << mine->second->GetName()
605 << "->" << mine->second->GetSolid()->GetName()
606 << " being placed in mother " << moms->second->GetName()
607 << "->" << moms->second->GetSolid()->GetName() << " at "
608 << fRef.fOrigin[0]*cm << ","
609 << fRef.fOrigin[1]*cm << ","
610 << fRef.fOrigin[2]*cm
611 << G4endlstd::endl;
612#endif
613
614 if (ilayer == 0) {
615 vpair_t mycopy(myvoluI,icopy);
616 fPhysicalVolumes[mycopy] = pvol;
617 fCurrentPlacement[myvoluI] = fPhysicalVolumes.find(mycopy);
618 fCurrentMother[myvoluI] = moms;
619 }
620 }
621
622 fPending = false;
623 }
624
625#ifdef LINUX_CPUTIME_PROFILING
626 timestr << " ( " << timer.getUserDelta() << " ) ";
627 G4cerr(*G4cerr_p) << timestr.str() << G4endlstd::endl;
628#endif
629 return icopy;
630}
631
632std::map<HddsG4Builder::vpair_t,G4LogicalVolume*>::iterator
633HddsG4Builder::addNewLayer(int volume_id, int layer)
634{
635#ifdef LINUX_CPUTIME_PROFILING
636 std::stringstream timestr;
637 timestr << "addNewLayer: " << timer.getUserTime() << ", "
638 << timer.getSystemTime() << ", " << timer.getRealTime();
639 timer.resetClocks();
640#endif
641
642 // Check if the volume already exists on the specified layer,
643 // and if so, return that, otherwise obtain the maximum layer
644 // index that exists already for this volume.
645
646 int max_layer = -1;
647 std::map<vpair_t,G4LogicalVolume*>::iterator volex;
648 for (volex = fLogicalVolumes.find(vpair_t(volume_id,0));
649 volex != fLogicalVolumes.end() &&
650 volex->first.first == volume_id;
651 ++volex)
652 {
653 max_layer = volex->first.second;
654 if (max_layer == layer)
655 {
656 return volex;
657 }
658 }
659 if (max_layer == -1)
660 {
661 G4cerr(*G4cerr_p) << "HddsG4Builder::addNewLayer called for volume "
662 << volume_id << " which does not (yet) exist!"
663 << G4endlstd::endl
664 << "Unable to continue, giving up."
665 << G4endlstd::endl;
666 exit(1);
667 }
668
669 // The volume does not exist yet on the specified layer, so we
670 // want to create a reflection of the existing layer 0 volume
671 // and place it on the specified layer. Before that can happen,
672 // the volume's mother must exist on the specified layer.
673 // Use recursion to ensure that this is always true.
674
675 std::map<vpair_t,G4LogicalVolume*>::iterator moms;
676 if (fCurrentMother.find(volume_id) != fCurrentMother.end())
677 {
678 moms = fCurrentMother[volume_id];
679 moms = addNewLayer(moms->first.first, moms->first.second + layer);
680 }
681 else
682 {
683 moms = fLogicalVolumes.end();
684 }
685
686 // Create a new copy of the existing max_layer volume. If the
687 // new layer index is higher than max_layer then set the interior
688 // filling material to null in the max_layer copy, otherwise
689 // the new copy should have a null material.
690
691 vpair_t newvol(volume_id,layer);
692 vpair_t oldvol(volume_id,max_layer);
693 G4VSolid *solid = fLogicalVolumes[oldvol]->GetSolid()->Clone();
694 G4String name(fLogicalVolumes[vpair_t(volume_id,0)]->GetName());
695 std::stringstream str;
696 str << name << "::" << layer;
697 solid->SetName(str.str());
698 G4Material *material = fLogicalVolumes[oldvol]->GetMaterial();
699 G4FieldManager *fieldmgr = fLogicalVolumes[oldvol]->GetFieldManager();
700 if (layer > max_layer)
701 {
702 fLogicalVolumes[newvol] = new G4LogicalVolume(solid,material,
703 str.str(),fieldmgr);
704 fLogicalVolumes[oldvol]->SetMaterial(0);
705 }
706 else
707 {
708 fLogicalVolumes[newvol] = new G4LogicalVolume(solid,0,str.str(),fieldmgr);
709 }
710
711 // If this is the top-level (world) volume then it cannot be
712 // placed here, so simply return here.
713
714 if (moms == fLogicalVolumes.end())
715 {
716 return fLogicalVolumes.find(newvol);
717 }
718
719 // Now we have a copy of the mother volume on the desired layer,
720 // and we have a new instance of our volume to place inside her.
721 // Two cases are supported in the present code, although more
722 // could be easily added.
723
724 // case 1: G4PVPlacement
725
726 if (fCurrentPlacement.find(volume_id) != fCurrentPlacement.end())
727 {
728 std::map<vpair_t,G4VPhysicalVolume*>::iterator placement;
729 placement = fCurrentPlacement[volume_id];
730 G4PVPlacement *player0 = (G4PVPlacement*)placement->second;
731 vpair_t mycopy(volume_id,placement->first.second);
732 G4PVPlacement *playerN;
733 playerN = new G4PVPlacement(player0->GetRotation(),
Value stored to 'playerN' is never read
734 player0->GetTranslation(),
735 fLogicalVolumes[newvol],
736 str.str(),
737 moms->second,0,
738 placement->first.second);
739#ifdef CHECK_OVERLAPS_MM
740 playerN->CheckOverlaps(1000,CHECK_OVERLAPS_MM);
741#endif
742#ifdef DEBUG_PLACEMENT
743 G4cout(*G4cout_p) << "volume " << str.str()
744 << "->" << fLogicalVolumes[newvol]->GetName()
745 << "->" << fLogicalVolumes[newvol]->GetSolid()->GetName()
746 << " being placed in mother " << moms->second->GetName()
747 << "->" << moms->second->GetSolid()->GetName()
748 << G4endlstd::endl;
749#endif
750 }
751
752 // case 2: G4PVDivision
753
754 if (fCurrentDivision.find(volume_id) != fCurrentDivision.end())
755 {
756 std::map<vpair_t,G4VPhysicalVolume*>::iterator division;
757 division = fCurrentDivision[volume_id];
758 G4PVDivision *division0 = (G4PVDivision*)division->second;
759 EAxis axis;
760 G4int ndiv;
761 G4double width;
762 G4double offset;
763 G4bool consuming;
764 division0->GetReplicationData(axis,ndiv,width,offset,consuming);
765 // axis returned by GetReplicationData is sometimes wrong!
766 axis = division0->GetDivisionAxis();
767 vpair_t mydiv(volume_id,division->first.second);
768 G4PVDivision *divisionN;
769 divisionN = new G4PVDivision(str.str(),
770 fLogicalVolumes[newvol],
771 moms->second,
772 axis, ndiv, width, offset);
773 fLogicalVolumes[newvol]->SetVisAttributes(new G4VisAttributes(false));
774#ifdef DEBUG_PLACEMENT
775 G4cout(*G4cout_p) << ndiv << " copies of division " << str.str()
776 << "->" << fLogicalVolumes[newvol]->GetName()
777 << "->" << fLogicalVolumes[newvol]->GetSolid()->GetName()
778 << " being placed in mother " << moms->second->GetName()
779 << "->" << moms->second->GetSolid()->GetName()
780 << " offset by " << offset << " on axis " << axis
781 << " and repeating every " << width
782 << G4endlstd::endl;
783#endif
784 }
785
786#ifdef LINUX_CPUTIME_PROFILING
787 timestr << " ( " << timer.getUserDelta() << " ) ";
788 G4cerr(*G4cerr_p) << timestr.str() << G4endlstd::endl;
789#endif
790 return fLogicalVolumes.find(newvol);
791}
792
793int HddsG4Builder::createDivision(XString& divStr, Refsys& ref)
794{
795#ifdef LINUX_CPUTIME_PROFILING
796 std::stringstream timestr;
797 timestr << "createDivision: " << timer.getUserTime() << ", "
798 << timer.getSystemTime() << ", " << timer.getRealTime();
799 timer.resetClocks();
800#endif
801 int ndiv = CodeWriter::createDivision(divStr,ref);
802
803 XString motherS(ref.fMother->getAttribute(X("HDDSvolu")XString("HDDSvolu").unicode_str()));
804 XString myvoluS(ref.fPartition.divEl->getAttribute(X("HDDSvolu")XString("HDDSvolu").unicode_str()));
805 int motherI = atoi(S(motherS)motherS.c_str());
806 int myvoluI = atoi(S(myvoluS)myvoluS.c_str());
807 std::map<vpair_t,G4LogicalVolume*>::iterator moms;
808 moms = fLogicalVolumes.find(vpair_t(motherI,0));
809 G4VSolid* solid = moms->second->GetSolid()->Clone();
810 solid->SetName(divStr);
811 EAxis axis;
812 double width,offset;
813
814 if (ref.fPartition.axis == "x")
815 {
816 axis = kXAxis;
817 width = ref.fPartition.step * cm;
818 offset = ref.fPartition.offset * cm;
819 }
820 else if (ref.fPartition.axis == "y")
821 {
822 axis = kYAxis;
823 width = ref.fPartition.step * cm;
824 offset = ref.fPartition.offset * cm;
825 }
826 else if (ref.fPartition.axis == "z")
827 {
828 axis = kZAxis;
829 width = ref.fPartition.step * cm;
830 offset = ref.fPartition.offset * cm;
831 }
832 else if (ref.fPartition.axis == "rho")
833 {
834 axis = kRho;
835 width = ref.fPartition.step * cm;
836 offset = ref.fPartition.offset * cm;
837 }
838 else if (ref.fPartition.axis == "phi")
839 {
840 axis = kPhi;
841 width = ref.fPartition.step * deg;
842 offset = ref.fPartition.offset * deg;
843 }
844 else
845 {
846 XString nameS(ref.fMother->getAttribute(X("name")XString("name").unicode_str()));
847 G4cerr(*G4cerr_p) << APP_NAME"HddsG4Builder" << " error: Volume division is requested "
848 << " for volume " << S(nameS)nameS.c_str() << " but division axis "
849 << "\"" << ref.fPartition.axis << "\""
850 << " is not currently supported."
851 << G4endlstd::endl;
852 exit(1);
853 }
854
855 vpair_t mothvol(moms->first);
856 vpair_t myvol(myvoluI,moms->first.second);
857 G4Material* material = moms->second->GetMaterial();
858 G4FieldManager* fieldmgr = fFieldManagers[ref.fRegionID];
859 fLogicalVolumes[myvol] = new G4LogicalVolume(solid,material,
860 S(divStr)divStr.c_str(),fieldmgr);
861 fLogicalVolumes[myvol]->SetVisAttributes(new G4VisAttributes(false));
862 vpair_t mydiv(myvoluI,0);
863 fPhysicalVolumes[mydiv] = new G4PVDivision(S(divStr)divStr.c_str(),
864 fLogicalVolumes[myvol],
865 fLogicalVolumes[mothvol],
866 axis,ndiv,width,offset);
867 fCurrentDivision[myvoluI] = fPhysicalVolumes.find(mydiv);
868 fCurrentMother[myvoluI] = moms;
869
870 // recompute the limits of the solid representing this division
871
872 G4VPVParameterisation *param = fPhysicalVolumes[mydiv]
873 ->GetParameterisation();
874 if (dynamic_cast<G4Box*>(solid))
875 {
876 param->ComputeDimensions(*(G4Box*)solid,0,fPhysicalVolumes[mydiv]);
877 }
878 else if (dynamic_cast<G4Tubs*>(solid))
879 {
880 param->ComputeDimensions(*(G4Tubs*)solid,0,fPhysicalVolumes[mydiv]);
881 double start = ((G4Tubs*)solid)->GetStartPhiAngle();
882 double delta = ((G4Tubs*)solid)->GetDeltaPhiAngle();
883 fCurrentPhiCenter[myvoluI] = start + delta/2;
884 }
885 else if (dynamic_cast<G4Trd*>(solid))
886 {
887 param->ComputeDimensions(*(G4Trd*)solid,0,fPhysicalVolumes[mydiv]);
888 }
889 else if (dynamic_cast<G4Trap*>(solid))
890 {
891 param->ComputeDimensions(*(G4Trap*)solid,0,fPhysicalVolumes[mydiv]);
892 }
893 else if (dynamic_cast<G4Cons*>(solid))
894 {
895 param->ComputeDimensions(*(G4Cons*)solid,0,fPhysicalVolumes[mydiv]);
896 double start = ((G4Cons*)solid)->GetStartPhiAngle();
897 double delta = ((G4Cons*)solid)->GetDeltaPhiAngle();
898 fCurrentPhiCenter[myvoluI] = start + delta/2;
899 }
900 else if (dynamic_cast<G4Sphere*>(solid))
901 {
902 param->ComputeDimensions(*(G4Sphere*)solid,0,fPhysicalVolumes[mydiv]);
903 double start = ((G4Sphere*)solid)->GetStartPhiAngle();
904 double delta = ((G4Sphere*)solid)->GetDeltaPhiAngle();
905 fCurrentPhiCenter[myvoluI] = start + delta/2;
906 }
907 else if (dynamic_cast<G4Orb*>(solid))
908 {
909 param->ComputeDimensions(*(G4Orb*)solid,0,fPhysicalVolumes[mydiv]);
910 }
911 else if (dynamic_cast<G4Torus*>(solid))
912 {
913 param->ComputeDimensions(*(G4Torus*)solid,0,fPhysicalVolumes[mydiv]);
914 double start = ((G4Torus*)solid)->GetSPhi();
915 double delta = ((G4Torus*)solid)->GetDPhi();
916 fCurrentPhiCenter[myvoluI] = start + delta/2;
917 }
918 else if (dynamic_cast<G4Para*>(solid))
919 {
920 param->ComputeDimensions(*(G4Para*)solid,0,fPhysicalVolumes[mydiv]);
921 }
922 else if (dynamic_cast<G4Polycone*>(solid))
923 {
924 param->ComputeDimensions(*(G4Polycone*)solid,0,fPhysicalVolumes[mydiv]);
925 G4PolyconeHistorical *params;
926 params = ((G4Polycone*)solid)->GetOriginalParameters();
927 double start = params->Start_angle;
928 double delta = params->Opening_angle;
929 fCurrentPhiCenter[myvoluI] = start + delta/2;
930 }
931 else if (dynamic_cast<G4Polyhedra*>(solid))
932 {
933 param->ComputeDimensions(*(G4Polyhedra*)solid,0,fPhysicalVolumes[mydiv]);
934 G4PolyhedraHistorical *params;
935 params = ((G4Polyhedra*)solid)->GetOriginalParameters();
936 double start = params->Start_angle;
937 double delta = params->Opening_angle;
938 fCurrentPhiCenter[myvoluI] = start + delta/2;
939 }
940 else if (dynamic_cast<G4Hype*>(solid))
941 {
942 param->ComputeDimensions(*(G4Hype*)solid,0,fPhysicalVolumes[mydiv]);
943 }
944 else
945 {
946 G4cerr(*G4cerr_p) << "HddsG4Builder::createDivision called for unsupported "
947 << "solid type!" << G4endlstd::endl
948 << "Unable to continue, giving up."
949 << G4endlstd::endl;
950 exit(1);
951 }
952#ifdef DEBUG_PLACEMENT
953 G4cout(*G4cout_p) << ndiv << " copies of division " << divStr
954 << "->" << fLogicalVolumes[myvol]->GetName()
955 << "->" << fLogicalVolumes[myvol]->GetSolid()->GetName()
956 << " being placed in mother " << fLogicalVolumes[mothvol]->GetName()
957 << "->" << fLogicalVolumes[mothvol]->GetSolid()->GetName()
958 << " offset by " << offset << " on axis " << axis
959 << " and repeating every " << width
960 << G4endlstd::endl;
961#endif
962
963#ifdef LINUX_CPUTIME_PROFILING
964 timestr << " ( " << timer.getUserDelta() << " ) ";
965 G4cerr(*G4cerr_p) << timestr.str() << G4endlstd::endl;
966#endif
967 return ndiv;
968}
969
970int HddsG4Builder::createRegion(DOMElement* el, Refsys& ref)
971{
972#ifdef LINUX_CPUTIME_PROFILING
973 std::stringstream timestr;
974 timestr << "createRegion: " << timer.getUserTime() << ", "
975 << timer.getSystemTime() << ", " << timer.getRealTime();
976 timer.resetClocks();
977#endif
978 int iregion = CodeWriter::createRegion(el,ref);
979
980 XString motherS(ref.fMother->getAttribute(X("HDDSvolu")XString("HDDSvolu").unicode_str()));
981 int motherI = atoi(S(motherS)motherS.c_str());
982 std::map<vpair_t,G4LogicalVolume*>::iterator moms;
983 moms = fLogicalVolumes.find(vpair_t(motherI,0));
984
985 G4double *R = (G4double*)ref.fMRmatrix;
986 G4double *O = ref.fMOrigin;
987 G4AffineTransform xform(CLHEP::HepRotation(CLHEP::HepRep3x3(R)),
988 CLHEP::Hep3Vector(O[0],O[1],O[2]));
989
990 if (ref.fRegion)
991 {
992 DOMNodeList* noBfieldL;
993 DOMNodeList* uniBfieldL;
994 DOMNodeList* mapBfieldL;
995 DOMNodeList* compBfieldL;
996 DOMNodeList* swimL;
997 noBfieldL = ref.fRegion->getElementsByTagName(X("noBfield")XString("noBfield").unicode_str());
998 uniBfieldL = ref.fRegion->getElementsByTagName(X("uniformBfield")XString("uniformBfield").unicode_str());
999 mapBfieldL = ref.fRegion->getElementsByTagName(X("mappedBfield")XString("mappedBfield").unicode_str());
1000 compBfieldL = ref.fRegion->getElementsByTagName(X("computedBfield")XString("computedBfield").unicode_str());
1001 swimL = ref.fRegion->getElementsByTagName(X("swim")XString("swim").unicode_str());
1002 double maxArcStep = 0;
1003 XString methodS("helix");
1004 if (swimL->getLength() > 0)
1005 {
1006 DOMElement* swimEl = (DOMElement*)swimL->item(0);
1007 methodS = swimEl->getAttribute(X("method")XString("method").unicode_str());
1008 Units unit;
1009 unit.getConversions(swimEl);
1010 XString stepS(swimEl->getAttribute(X("maxArcStep")XString("maxArcStep").unicode_str()));
1011 maxArcStep = atof(S(stepS)stepS.c_str())/unit.rad;
1012 }
1013 if (noBfieldL->getLength() > 0)
1014 {
1015 fMagneticRegions[iregion] = 0;
1016 G4ThreeVector Bvec(0,0,1e-99);
1017 G4UniformMagField *fld = new G4UniformMagField(Bvec);
1018 G4Mag_EqRhs *eqn = new G4Mag_UsualEqRhs(fld);
1019 G4MagIntegratorStepper *stepper = new G4ExactHelixStepper(eqn);
1020 G4ChordFinder *cfinder = new G4ChordFinder(fld, 1e+99, stepper);
1021 fFieldManagers[iregion] = new G4FieldManager(fld, cfinder);
1022 }
1023 else if (uniBfieldL->getLength() > 0)
1024 {
1025 Units funit;
1026 DOMElement* uniBfieldEl = (DOMElement*)uniBfieldL->item(0);
1027 funit.getConversions(uniBfieldEl);
1028 XString bvecS(uniBfieldEl->getAttribute(X("Bx_By_Bz")XString("Bx_By_Bz").unicode_str()));
1029 std::stringstream str(S(bvecS)bvecS.c_str());
1030 double B[3];
1031 str >> B[0] >> B[1] >> B[2];
1032 double u = kilogauss/funit.kG;
1033 G4ThreeVector Bvec(B[0],B[1],B[2]);
1034 GlueXUniformMagField *fld = new GlueXUniformMagField(Bvec,u,xform);
1035 fMagneticRegions[iregion] = fld;
1036 G4Mag_EqRhs *eqn = new G4Mag_UsualEqRhs(fld);
1037 G4MagIntegratorStepper *stepper = new G4ExactHelixStepper(eqn);
1038 G4ChordFinder *cfinder = new G4ChordFinder(fld, 0.01, stepper);
1039 fFieldManagers[iregion] = new G4FieldManager(fld, cfinder);
1040 }
1041 else if (mapBfieldL->getLength() > 0)
1042 {
1043 Units funit;
1044 DOMElement* mapBfieldEl = (DOMElement*)mapBfieldL->item(0);
1045 funit.getConversions(mapBfieldEl);
1046 XString bvecS(mapBfieldEl->getAttribute(X("maxBfield")XString("maxBfield").unicode_str()));
1047 std::stringstream str(S(bvecS)bvecS.c_str());
1048 double Bmax;
1049 str >> Bmax;
1050 double u = kilogauss/funit.kG;
1051 GlueXMappedMagField *fld = new GlueXMappedMagField(Bmax,u,xform);
1052 fMagneticRegions[iregion] = fld;
1053 G4Mag_EqRhs *eqn = new G4Mag_UsualEqRhs(fld);
1054 G4MagIntegratorStepper *stepper;
1055 if (methodS == "RungeKutta") {
1056 stepper = new G4ClassicalRK4(eqn);
1057 }
1058 else {
1059 stepper = new G4HelixMixedStepper(eqn);
1060 }
1061 G4ChordFinder *cfinder = new G4ChordFinder(fld, 0.01, stepper);
1062 if (maxArcStep > 0) {
1063 double rmin = 0.1 / (0.03 * Bmax * u) * meter;
1064 double max_miss = rmin * (1 - cos(maxArcStep / 2));
1065 cfinder->SetDeltaChord(max_miss);
1066 }
1067 fFieldManagers[iregion] = new G4FieldManager(fld, cfinder);
1068 }
1069 else if (compBfieldL->getLength() > 0)
1070 {
1071 Units funit;
1072 DOMElement* compBfieldEl = (DOMElement*)compBfieldL->item(0);
1073 funit.getConversions(compBfieldEl);
1074 XString bvecS(compBfieldEl->getAttribute(X("maxBfield")XString("maxBfield").unicode_str()));
1075 std::stringstream str(S(bvecS)bvecS.c_str());
1076 double Bmax;
1077 str >> Bmax;
1078 double u = kilogauss/funit.kG;
1079 GlueXComputedMagField *fld = new GlueXComputedMagField(Bmax,u,xform);
1080 fld->SetFunction(XString(compBfieldEl->getAttribute(X("function")XString("function").unicode_str())));
1081 fMagneticRegions[iregion] = fld;
1082 G4Mag_EqRhs *eqn = new G4Mag_UsualEqRhs(fld);
1083 G4MagIntegratorStepper *stepper;
1084 if (methodS == "RungeKutta") {
1085 stepper = new G4ClassicalRK4(eqn);
1086 }
1087 else {
1088 stepper = new G4HelixMixedStepper(eqn);
1089 }
1090 G4ChordFinder *cfinder = new G4ChordFinder(fld, 0.01, stepper);
1091 if (maxArcStep > 0) {
1092 double rmin = 0.1 / (0.03 * Bmax * u) * meter;
1093 double max_miss = rmin * (1 - cos(maxArcStep / 2));
1094 cfinder->SetDeltaChord(max_miss);
1095 }
1096 fFieldManagers[iregion] = new G4FieldManager(fld, cfinder);
1097 }
1098 }
1099
1100#ifdef LINUX_CPUTIME_PROFILING
1101 timestr << " ( " << timer.getUserDelta() << " ) ";
1102 G4cerr(*G4cerr_p) << timestr.str() << G4endlstd::endl;
1103#endif
1104 return iregion;
1105}
1106
1107void HddsG4Builder::createSetFunctions(DOMElement* el, const XString& ident)
1108{
1109#ifdef LINUX_CPUTIME_PROFILING
1110 std::stringstream timestr;
1111 timestr << "createSetFunctions: " << timer.getUserTime() << ", "
1112 << timer.getSystemTime() << ", " << timer.getRealTime();
1113 timer.resetClocks();
1114#endif
1115 CodeWriter::createSetFunctions(el,ident);
1116
1117#ifdef LINUX_CPUTIME_PROFILING
1118 timestr << " ( " << timer.getUserDelta() << " ) ";
1119 G4cerr(*G4cerr_p) << timestr.str() << G4endlstd::endl;
1120#endif
1121}
1122
1123void HddsG4Builder::createGetFunctions(DOMElement* el, const XString& ident)
1124{
1125#ifdef LINUX_CPUTIME_PROFILING
1126 std::stringstream timestr;
1127 timestr << "createGetFunctions: " << timer.getUserTime() << ", "
1128 << timer.getSystemTime() << ", " << timer.getRealTime();
1129 timer.resetClocks();
1130#endif
1131 CodeWriter::createGetFunctions(el,ident);
1132
1133 std::vector<int> table;
1134 std::vector<int> start;
1135 start.push_back(0);
1136
1137 XString funcNameStr;
1138 XString identCaps(ident);
1139 identCaps[0] = toupper(identCaps[0]);
1140 funcNameStr = "get" + identCaps;
1141 for (int ivolu = 1; ivolu <= Refsys::fVolumes; ivolu++)
1142 {
1143 start.push_back(0);
1144 int ncopy = Refsys::fIdentifierTable[ivolu]["copy counter"].back();
1145 std::map<std::string,std::vector<int> >::iterator idlist =
1146 Refsys::fIdentifierTable[ivolu].find(ident);
1147 if (idlist != Refsys::fIdentifierTable[ivolu].end())
1148 {
1149 if (ncopy != (int)idlist->second.size())
1150 {
1151 G4cerr(*G4cerr_p) << APP_NAME"HddsG4Builder" << " warning: volume " << ivolu
1152 << " has " << ncopy << " copies, but "
1153 << idlist->second.size() << " "
1154 << ident << " identifiers!" << G4endlstd::endl;
1155 for (int idx = 0; idx < (int)idlist->second.size(); idx++)
1156 {
1157 G4cerr(*G4cerr_p) << idlist->second[idx] << " ";
1158 if (idx/20*20 == idx)
1159 G4cerr(*G4cerr_p) << G4endlstd::endl;
1160 }
1161 G4cerr(*G4cerr_p) << G4endlstd::endl;
1162 }
1163 start[ivolu] = table.size() + 1;
1164 for (int icopy = 0; icopy < ncopy; icopy++)
1165 {
1166 table.push_back(idlist->second[icopy]);
1167 }
1168 }
1169 else
1170 {
1171 start[ivolu] = 0;
1172 }
1173 }
1174
1175#ifdef LINUX_CPUTIME_PROFILING
1176 timestr << " ( " << timer.getUserDelta() << " ) ";
1177 G4cerr(*G4cerr_p) << timestr.str() << G4endlstd::endl;
1178#endif
1179}
1180
1181void HddsG4Builder::createMapFunctions(DOMElement* el, const XString& ident)
1182{
1183#ifdef LINUX_CPUTIME_PROFILING
1184 std::stringstream timestr;
1185 timestr << "createMapFunctions: " << timer.getUserTime() << ", "
1186 << timer.getSystemTime() << ", " << timer.getRealTime();
1187 timer.resetClocks();
1188#endif
1189 CodeWriter::createMapFunctions(el,ident);
1190
1191 if (el == 0)
1192 {
1193 return;
1194 }
1195 DOMNodeList* regionsL = el->getOwnerDocument()->getDocumentElement()
1196 ->getElementsByTagName(X("regions")XString("regions").unicode_str());
1197 if (regionsL->getLength() == 0)
1198 {
1199 return;
1200 }
1201 DOMElement* regionsEl = (DOMElement*)regionsL->item(0);
1202 DOMNodeList* regionL = regionsEl->getElementsByTagName(X("region")XString("region").unicode_str());
1203 for (int ireg=0; ireg < (int)regionL->getLength(); ++ireg)
1204 {
1205 DOMElement* regionEl = (DOMElement*)regionL->item(ireg);
1206 DOMNodeList* mapfTagL = regionEl->getElementsByTagName(X("mappedBfield")XString("mappedBfield").unicode_str());
1207 if (mapfTagL->getLength() == 0)
1208 continue;
1209 DOMElement* mapfEl = (DOMElement*)mapfTagL->item(0);
1210 XString nameS(regionEl->getAttribute(X("name")XString("name").unicode_str()));
1211 XString iregionS(regionEl->getAttribute(X("HDDSregion")XString("HDDSregion").unicode_str()));
1212 if (iregionS.size() == 0)
1213 continue;
1214 int iregion = atoi(S(iregionS)iregionS.c_str());
1215 GlueXMappedMagField *magfield = (GlueXMappedMagField*)
1216 fMagneticRegions[iregion];
1217 int axorder[] = {0,0,0,0};
1218 int axsamples[] = {0,0,0,0};
1219
1220 XString gridtype;
1221 DOMNodeList* gridL = mapfEl->getElementsByTagName(X("grid")XString("grid").unicode_str());
1222 int ngrid;
1223 for (ngrid = 0; ngrid < (int)gridL->getLength(); ++ngrid)
1224 {
1225 int axsense[] = {1,1,1,1};
1226 double axlower[4], axupper[4];
1227 DOMElement* gridEl = (DOMElement*)gridL->item(ngrid);
1228 XString typeS(gridEl->getAttribute(X("type")XString("type").unicode_str()));
1229 if (gridtype.size() > 0 && typeS != gridtype)
1230 {
1231 G4cerr(*G4cerr_p) << APP_NAME"HddsG4Builder" << " error: mappedBfield in region " << S(nameS)nameS.c_str()
1232 << " superimposes incompatible grid types." << G4endlstd::endl;
1233 exit(1);
1234 }
1235 gridtype = typeS;
1236
1237 DOMNodeList* samplesL = gridEl->getElementsByTagName(X("samples")XString("samples").unicode_str());
1238 if (samplesL->getLength() != 3)
1239 {
1240 G4cerr(*G4cerr_p) << APP_NAME"HddsG4Builder" << " error: mappedBfield in region " << S(nameS)nameS.c_str()
1241 << " does not have samples for three axes." << G4endlstd::endl;
1242 exit(1);
1243 }
1244
1245 for (int iax = 1; iax <= 3; ++iax)
1246 {
1247 DOMElement* sampleEl = (DOMElement*)samplesL->item(iax-1);
1248 XString nS(sampleEl->getAttribute(X("n")XString("n").unicode_str()));
1249 XString axisS(sampleEl->getAttribute(X("axis")XString("axis").unicode_str()));
1250 XString boundsS(sampleEl->getAttribute(X("bounds")XString("bounds").unicode_str()));
1251 XString senseS(sampleEl->getAttribute(X("sense")XString("sense").unicode_str()));
1252 Units sunit;
1253 double bound[2];
1254 sunit.getConversions(sampleEl);
1255 std::stringstream listr(boundsS);
1256 listr >> bound[0] >> bound[1];
1257 int iaxis=0;
1258 if (gridtype == "cartesian")
1259 {
1260 if (axisS == "x" &&
1261 (axorder[0] == 0 || axorder[0] == iax))
1262 {
1263 iaxis = 1;
1264 axorder[0] = iax;
1265 bound[0] /= sunit.cm;
1266 bound[1] /= sunit.cm;
1267 }
1268 else if (axisS == "y" &&
1269 (axorder[1] == 0 || axorder[1] == iax))
1270 {
1271 iaxis = 2;
1272 axorder[1] = iax;
1273 bound[0] /= sunit.cm;
1274 bound[1] /= sunit.cm;
1275 }
1276 else if (axisS == "z" &&
1277 (axorder[2] == 0 || axorder[2] == iax))
1278 {
1279 iaxis = 3;
1280 axorder[2] = iax;
1281 bound[0] /= sunit.cm;
1282 bound[1] /= sunit.cm;
1283 }
1284 else
1285 {
1286 G4cerr(*G4cerr_p) << APP_NAME"HddsG4Builder" << " error: grid in region " << S(nameS)nameS.c_str()
1287 << " contains an incompatible set of samples."
1288 << G4endlstd::endl;
1289 exit(1);
1290 }
1291 }
1292 else if (gridtype == "cylindrical")
1293 {
1294 if (axisS == "r" && axorder[0] == 0)
1295 {
1296 iaxis = 1;
1297 axorder[0] = iax;
1298 bound[0] /= sunit.cm;
1299 bound[1] /= sunit.cm;
1300 }
1301 else if (axisS == "phi" && axorder[1] == 0)
1302 {
1303 iaxis = 2;
1304 axorder[1] = iax;
1305 bound[0] /= sunit.rad;
1306 bound[1] /= sunit.rad;
1307 }
1308 else if (axisS == "z" && axorder[2] == 0)
1309 {
1310 iaxis = 3;
1311 axorder[2] = iax;
1312 bound[0] /= sunit.cm;
1313 bound[1] /= sunit.cm;
1314 }
1315 else
1316 {
1317 G4cerr(*G4cerr_p) << APP_NAME"HddsG4Builder" << " error: grid in region " << S(nameS)nameS.c_str()
1318 << " contains an incompatible set of samples."
1319 << G4endlstd::endl;
1320 exit(1);
1321 }
1322 }
1323
1324 int n = atoi(S(nS)nS.c_str());
1325 if (axsamples[iaxis] == 0 || axsamples[iaxis] == n)
1326 {
1327 axsamples[iaxis] = n;
1328 }
1329 else
1330 {
1331 G4cerr(*G4cerr_p) << APP_NAME"HddsG4Builder" << " error: mappedBfield in region " << S(nameS)nameS.c_str()
1332 << " combines incompatible grid elements."
1333 << G4endlstd::endl;
1334 exit(1);
1335 }
1336
1337 axlower[iaxis] = bound[0];
1338 axupper[iaxis] = bound[1];
1339 if (senseS == "reverse")
1340 {
1341 axsense[iaxis] = -1;
1342 }
1343 }
1344
1345 if (gridtype == "cartesian")
1346 {
1347 magfield->AddCartesianGrid(axsamples, axorder, axsense,
1348 axlower, axupper);
1349 }
1350 else if (gridtype == "cylindrical")
1351 {
1352 magfield->AddCylindricalGrid(axsamples, axorder, axsense,
1353 axlower, axupper);
1354 }
1355 else
1356 {
1357 G4cerr(*G4cerr_p) << APP_NAME"HddsG4Builder" << " error: unrecognized grid type "
1358 << S(gridtype)gridtype.c_str() << G4endlstd::endl;
1359 exit(1);
1360 }
1361 }
1362
1363 XString mapS(mapfEl->getAttribute(X("map")XString("map").unicode_str()));
1364 XString encS(mapfEl->getAttribute(X("encoding")XString("encoding").unicode_str()));
1365 if (encS != "utf-8")
1366 {
1367 G4cerr(*G4cerr_p) << APP_NAME"HddsG4Builder" << " error: mappedBfield in region " << S(nameS)nameS.c_str()
1368 << " uses unsupported encoding " << encS << G4endlstd::endl;
1369 exit(1);
1370 }
1371 else if (mapS.substr(0,7) != "file://")
1372 {
1373 G4cerr(*G4cerr_p) << APP_NAME"HddsG4Builder" << " error: mappedBfield in region " << S(nameS)nameS.c_str()
1374 << " uses unsupported map URL " << mapS << G4endlstd::endl;
1375 exit(1);
1376 }
1377 mapS.erase(0,7);
1378 magfield->ReadMapFile(S(mapS)mapS.c_str());
1379 }
1380
1381 // Apply a post-build fix to ensure that every point in the geometry
1382 // has a consistent magnetic field on all layers, using recursion.
1383
1384 addReflections(1);
1385
1386#ifdef LINUX_CPUTIME_PROFILING
1387 timestr << " ( " << timer.getUserDelta() << " ) ";
1388 G4cerr(*G4cerr_p) << timestr.str() << G4endlstd::endl;
1389#endif
1390}
1391
1392void HddsG4Builder::translate(DOMElement* topel)
1393{
1394#ifdef LINUX_CPUTIME_PROFILING
1395 std::stringstream timestr;
1396 timestr << "translate: " << timer.getUserTime() << ", "
1397 << timer.getSystemTime() << ", " << timer.getRealTime();
1398 timer.resetClocks();
1399#endif
1400
1401 CodeWriter::translate(topel);
1402 DOMDocument* document = topel->getOwnerDocument();
1403 DOMElement* worldEl = document->getElementById(X("everything")XString("everything").unicode_str());
1404 XString topvolS(worldEl->getAttribute(X("envelope")XString("envelope").unicode_str()));
1405 DOMElement* topvolEl = document->getElementById(X(topvolS)XString(topvolS).unicode_str());
1406 XString ivoluS(topvolEl->getAttribute(X("HDDSvolu")XString("HDDSvolu").unicode_str()));
1407 if (ivoluS.size() == 0)
1408 {
1409 G4cerr(*G4cerr_p) << APP_NAME"HddsG4Builder" << " error: top-level composition named "
1410 << "\"everything\" was not found in input document!"
1411 << G4endlstd::endl;
1412 exit(1);
1413 }
1414 fWorldVolume = atoi(S(ivoluS)ivoluS.c_str());
1415
1416#ifdef LINUX_CPUTIME_PROFILING
1417 timestr << " ( " << timer.getUserDelta() << " ) ";
1418 G4cerr(*G4cerr_p) << timestr.str() << G4endlstd::endl;
1419#endif
1420}
1421
1422G4LogicalVolume* HddsG4Builder::getWorldVolume(int parallel) const
1423{
1424 int worlds = 0;
1425 std::map<vpair_t,G4LogicalVolume*>::const_iterator paraworld;
1426 for (paraworld = fLogicalVolumes.find(vpair_t(fWorldVolume,0));
1427 paraworld != fLogicalVolumes.end() &&
1428 paraworld->first.first == fWorldVolume;
1429 ++paraworld, ++worlds)
1430 {}
1431 if (parallel < worlds)
1432 {
1433 for (int p=0; p <= parallel; ++p, --paraworld) {}
1434 return paraworld->second;
1435 }
1436 return 0;
1437}
1438
1439void HddsG4Builder::addReflections(int volume_id)
1440{
1441#ifdef LINUX_CPUTIME_PROFILING
1442 std::stringstream timestr;
1443 timestr << "addReflections: " << timer.getUserTime() << ", "
1444 << timer.getSystemTime() << ", " << timer.getRealTime();
1445 timer.resetClocks();
1446#endif
1447 std::vector<int> layer;
1448 std::map<vpair_t,G4LogicalVolume*>::iterator witer;
1449 for (witer = fLogicalVolumes.find(vpair_t(1,0));
1450 witer != fLogicalVolumes.end() && witer->first.first == 1;
1451 ++witer)
1452 {
1453 layer.push_back(witer->first.second);
1454 }
1455
1456 std::map<int,G4FieldManager*> lastManager;
1457 std::map<int,G4FieldManager*>::iterator miter;
1458 std::map<vpair_t,G4LogicalVolume*>::iterator mine;
1459 for (mine = fLogicalVolumes.find(vpair_t(volume_id,0));
1460 mine != fLogicalVolumes.end() && mine->first.first == volume_id;
1461 ++mine)
1462 {
1463 G4FieldManager* fieldmgr = mine->second->GetFieldManager();
1464 for (int child = 0; child < (int)mine->second->GetNoDaughters(); ++child) {
1465 G4VPhysicalVolume* pvol = mine->second->GetDaughter(child);
1466 G4LogicalVolume* lvol = pvol->GetLogicalVolume();
1467 G4FieldManager* mgr = lvol->GetFieldManager();
1468 int childId = getVolumeId(lvol);
1469 miter = lastManager.find(childId);
1470 if ((miter != lastManager.end() && mgr != miter->second) ||
1471 (miter == lastManager.end() && mine->first.second != 0 &&
1472 mgr != 0 && mgr != fieldmgr))
1473 {
1474 G4cerr(*G4cerr_p) << "HddsG4Builder::addReflections warning - "
1475 << "local magnetic field assigned to volume "
1476 << lvol->GetName()
1477 << " placed on a buried layer in the geometry."
1478 << G4endlstd::endl
1479 << "This is a violation of the assumption made "
1480 << "in the Geant4 tracking algorithm that fields "
1481 << "at a given point are the same on all layers."
1482 << G4endlstd::endl
1483 << "Expect inconsistent results."
1484 << G4endlstd::endl;
1485 }
1486 else {
1487 lastManager[childId] = mgr;
1488 }
1489
1490 // Use recursion to visit every child logical volume once
1491
1492 if (miter == lastManager.end()) {
1493 addReflections(childId);
1494 }
1495
1496 // Any volume in the hierarchy that has a different magnetic field
1497 // manager than its mother volume needs to be reflected onto all
1498 // layers to maintain consistency of the field value at all points.
1499
1500 if (mgr != 0 && mgr != fieldmgr) {
1501 vpair_t cpy(childId, pvol->GetCopyNo());
1502 std::map<vpair_t,G4VPhysicalVolume*>::iterator piter;
1503 piter = fPhysicalVolumes.find(cpy);
1504 if (piter == fPhysicalVolumes.end()) {
1505 G4cerr(*G4cerr_p) << "HddsG4Builder::addReflections error - "
1506 << "physical volume " << pvol->GetName()
1507 << " copy " << cpy.second
1508 << " not found in physical volume lookup table, "
1509 << "cannot continue." << G4endlstd::endl;
1510 exit(1);
1511 }
1512 else if (piter->second->IsReplicated()) {
1513 fCurrentDivision[childId] = piter;
1514 }
1515 else {
1516 fCurrentPlacement[childId] = piter;
1517 }
1518 fCurrentMother[childId] = mine;
1519 std::vector<int>::iterator liter;
1520 for (liter = layer.begin(); liter != layer.end(); ++liter) {
1521 if (*liter > mine->first.second) {
1522#if DEBUG_REFLECTION
1523 G4cout(*G4cout_p) << "HddsG4Builder::addReflections info - "
1524 << "reflected volume " << lvol->GetName()
1525 << " copy " << cpy.second
1526 << " in mother " << mine->second->GetName()
1527 << " onto layer " << *liter
1528 << " because child field manager " << mgr
1529 << " is different from mother field manager " << fieldmgr
1530 << G4endlstd::endl;
1531#endif
1532 addNewLayer(childId, *liter);
1533 }
1534 }
1535 }
1536 }
1537 }
1538#ifdef LINUX_CPUTIME_PROFILING
1539 timestr << " ( " << timer.getUserDelta() << " ) ";
1540 G4cerr(*G4cerr_p) << timestr.str() << G4endlstd::endl;
1541#endif
1542}
1543
1544int HddsG4Builder::getVolumeId(G4LogicalVolume* vol) const
1545{
1546 std::map<vpair_t,G4LogicalVolume*>::const_iterator iter;
1547 for (iter = fLogicalVolumes.begin(); iter != fLogicalVolumes.end(); ++iter)
1548 {
1549 if (iter->second == vol) {
1550 return iter->first.first;
1551 }
1552 }
1553 G4cerr(*G4cerr_p) << "HddsG4Builder::getVolumeId error - "
1554 << "logical volume " << vol->GetName()
1555 << " not found in volume lookup table, "
1556 << "cannot continue." << G4endlstd::endl;
1557 exit(1);
1558}
1559
1560G4LogicalVolume* HddsG4Builder::getVolume(const G4String volname) const
1561{
1562 std::map<vpair_t,G4LogicalVolume*>::const_iterator iter;
1563 for (iter = fLogicalVolumes.begin(); iter != fLogicalVolumes.end(); ++iter)
1564 {
1565 if (iter->second->GetName() == volname) {
1566 return iter->second;
1567 }
1568 }
1569 return 0;
1570}
1571
1572const std::map<int, G4LogicalVolume*> HddsG4Builder::getSensitiveVolumes() const
1573{
1574 return fSensitiveVolumes;
1575}