GEMC  2.3
Geant4 Monte-Carlo Framework
detector.cc
Go to the documentation of this file.
1 // G4 headers
2 #include "G4Box.hh"
3 #include "G4Cons.hh"
4 #include "G4Ellipsoid.hh"
5 #include "G4IntersectionSolid.hh"
6 #include "G4NistManager.hh"
7 #include "G4Para.hh"
8 #include "G4Polycone.hh"
9 #include "G4Polyhedra.hh"
10 #include "G4Torus.hh"
11 #include "G4Trd.hh"
12 #include "G4Trap.hh"
13 #include "G4Tubs.hh"
14 #include "G4CutTubs.hh"
15 #include "G4EllipticalTube.hh"
16 #include "G4Paraboloid.hh"
17 #include "G4Hype.hh"
18 #include "G4Sphere.hh"
19 #include "G4GenericTrap.hh"
20 #include "G4SubtractionSolid.hh"
21 #include "G4UnionSolid.hh"
22 #include "G4IntersectionSolid.hh"
23 #include "G4PVReplica.hh"
24 #include "G4UnitsTable.hh"
25 #include "G4RotationMatrix.hh"
26 #include "G4TwoVector.hh"
27 
28 // gemc headers
29 #include "detector.h"
30 
31 // C++ headers
32 #include <string>
33 #include <vector>
34 using namespace std;
35 
36 // CLHEP units
37 #include "CLHEP/Units/PhysicalConstants.h"
38 using namespace CLHEP;
39 
40 detector::detector()
41 {
42  SolidV = NULL;
43  LogicV = NULL;
44  PhysicalV = NULL;
45 }
46 
47 int detector::create_solid(goptions gemcOpt, map<string, detector> *Map)
48 {
49  int built = 0;
50  if(SolidV) delete SolidV;
51 
52  string hd_msg = gemcOpt.optMap["LOG_MSG"].args + " Solid: >> ";
53  double VERB = gemcOpt.optMap["G4P_VERBOSITY"].arg ;
54  string catch_v = gemcOpt.optMap["CATCH"].args;
55 
56  if(type.find("ReplicaOf") != string::npos)
57  {
58  if(VERB>4 || name.find(catch_v) != string::npos)
59  cout << hd_msg << " " << name << " is a Replica. Solid Volume will not be built." << endl;
60  return 0;
61  }
62 
63  // ####
64  // Box
65  // ####
66  if(type == "Box")
67  {
68  if(dimensions.size() != 3)
69  {
70  cout << hd_msg << " Fatal Error: the number of dimensions for " << name
71  << " is " << dimensions.size() << ":" << endl;
72  for(unsigned int i=0; i<dimensions.size(); i++) cout << " dimension " << i + 1 << ": " << dimensions[i] << endl;
73  cout << " This does not match a G4Box. Exiting" << endl << endl;
74  exit(0);
75  }
76 
77  // checking BOX dimensions
78  for(unsigned int i=0; i<dimensions.size(); i++)
79  if(dimensions[i] == 0)
80  cout << " !!! Warning: BOX has one side null!" << endl;
81 
82  SolidV = new G4Box(name,
83  dimensions[0],
84  dimensions[1],
85  dimensions[2]);
86 
87  built = 1;
88  }
89 
90  // ##############
91  // Parallelepiped
92  // ##############
93  if(type == "Parallelepiped")
94  {
95  if(dimensions.size() != 6)
96  {
97  cout << hd_msg << " Fatal Error: the number of dimensions for " << name
98  << " is " << dimensions.size() << ":" << endl;
99  for(unsigned int i=0; i<dimensions.size(); i++) cout << " dimension " << i + 1 << ": " << dimensions[i] << endl;
100  cout << " This does not match a G4Para. Exiting" << endl << endl;
101  exit(0);
102  }
103 
104  SolidV = new G4Para(name,
105  dimensions[0],
106  dimensions[1],
107  dimensions[2],
108  dimensions[3],
109  dimensions[4],
110  dimensions[5]);
111 
112  built = 1;
113  }
114 
115 
116  // ######
117  // Sphere
118  // ######
119  if(type == "Sphere")
120  {
121  if(dimensions.size() != 6)
122  {
123  cout << hd_msg << " Fatal Error: the number of dimensions for " << name
124  << " is " << dimensions.size() << ":" << endl;
125  for(unsigned int i=0; i<dimensions.size(); i++) cout << " dimension " << i + 1 << ": " << dimensions[i] << endl;
126  cout << " This does not match a G4Sphere. Exiting." << endl << endl;
127  exit(0);
128  }
129 
130  SolidV = new G4Sphere(name,
131  dimensions[0],
132  dimensions[1],
133  dimensions[2],
134  dimensions[3],
135  dimensions[4],
136  dimensions[5]);
137 
138  built = 1;
139  }
140 
141 
142  // #########
143  // Ellipsoid
144  // #########
145  if(type == "Ellipsoid")
146  {
147  if(dimensions.size() != 5)
148  {
149  cout << hd_msg << " Fatal Error: the number of dimensions for " << name
150  << " is " << dimensions.size() << ":" << endl;
151  for(unsigned int i=0; i<dimensions.size(); i++) cout << " dimension " << i + 1 << ": " << dimensions[i] << endl;
152  cout << " This does not match a G4Ellipsoid. Exiting." << endl << endl;
153  exit(0);
154  }
155 
156  SolidV = new G4Ellipsoid(name,
157  dimensions[0],
158  dimensions[1],
159  dimensions[2],
160  dimensions[3],
161  dimensions[4]);
162 
163  built = 1;
164  }
165 
166  // ##########
167  // Paraboloid
168  // ##########
169  if(type == "Paraboloid")
170  {
171  if(dimensions.size() != 3)
172  {
173  cout << hd_msg << " Fatal Error: the number of dimensions for " << name
174  << " is " << dimensions.size() << ":" << endl;
175  for(unsigned int i=0; i<dimensions.size(); i++) cout << " dimension " << i + 1 << ": " << dimensions[i] << endl;
176  cout << " This does not match a G4Paraboloid. Exiting." << endl << endl;
177  exit(0);
178  }
179 
180  SolidV = new G4Paraboloid(name,
181  dimensions[0],
182  dimensions[1],
183  dimensions[2]);
184 
185  built = 1;
186  }
187 
188 
189  // ##################
190  // Hyperbolic Profile
191  // ##################
192  if(type == "Hype")
193  {
194  if(dimensions.size() != 5)
195  {
196  cout << hd_msg << " Fatal Error: the number of dimensions for " << name
197  << " is " << dimensions.size() << ":" << endl;
198  for(unsigned int i=0; i<dimensions.size(); i++) cout << " dimension " << i + 1 << ": " << dimensions[i] << endl;
199  cout << " This does not match a G4Hype. Exiting." << endl << endl;
200  exit(0);
201  }
202 
203  SolidV = new G4Hype(name,
204  dimensions[0],
205  dimensions[1],
206  dimensions[2],
207  dimensions[3],
208  dimensions[4]);
209 
210  built = 1;
211  }
212 
213 
214 
215  // ####
216  // Tube
217  // ####
218  if(type == "Tube")
219  {
220  if(dimensions.size() != 5)
221  {
222  cout << hd_msg << " Fatal Error: the number of dimensions for " << name
223  << " is " << dimensions.size() << ":" << endl;
224  for(unsigned int i=0; i<dimensions.size(); i++) cout << " dimension " << i + 1 << ": " << dimensions[i] << endl;
225  cout << " This does not match a G4Tubs. Exiting" << endl << endl;
226  exit(0);
227  }
228 
229  SolidV = new G4Tubs(name,
230  dimensions[0],
231  dimensions[1],
232  dimensions[2],
233  dimensions[3],
234  dimensions[4]);
235 
236  built = 1;
237  }
238 
239  // ####
240  // Cut Tube
241  // ####
242  if(type == "CTube")
243  {
244  if(dimensions.size() != 11)
245  {
246  cout << hd_msg << " Fatal Error: the number of dimensions for " << name
247  << " is " << dimensions.size() << ":" << endl;
248  for(unsigned int i=0; i<dimensions.size(); i++) cout << " dimension " << i + 1 << ": " << dimensions[i] << endl;
249  cout << " This does not match a G4CutTubs. Exiting" << endl << endl;
250  exit(0);
251  }
252 
253  SolidV = new G4CutTubs(name,
254  dimensions[0],
255  dimensions[1],
256  dimensions[2],
257  dimensions[3],
258  dimensions[4],
259  G4ThreeVector(dimensions[5], dimensions[6], dimensions[7]),
260  G4ThreeVector(dimensions[8], dimensions[9], dimensions[10]));
261 
262  built = 1;
263  }
264 
265  // ###############
266  // G4ElipticalTube
267  // ###############
268  if(type == "EllipticalTube" || type == "Eltu")
269  {
270  if(dimensions.size() != 3)
271  {
272  cout << hd_msg << " Fatal Error: the number of dimensions for " << name
273  << " is " << dimensions.size() << ":" << endl;
274  for(unsigned int i=0; i<dimensions.size(); i++) cout << " dimension " << i + 1 << ": " << dimensions[i] << endl;
275  cout << " This does not match a G4ElipticalTube. Exiting" << endl << endl;
276  exit(0);
277  }
278 
279  SolidV = new G4EllipticalTube(name,
280  dimensions[0], // dx: The equation of the surface in x/y is 1.0 = (x/dx)**2 +(y/dy)**2
281  dimensions[1], // dy
282  dimensions[2]); // dz = half length in z
283  built = 1;
284 
285  }
286 
287 
288 
289  // ####
290  // Cone
291  // ####
292  if(type == "Cons")
293  {
294  if(dimensions.size() != 7)
295  {
296  cout << hd_msg << " Fatal Error: the number of dimensions for " << name
297  << " is " << dimensions.size() << ":" << endl;
298  for(unsigned int i=0; i<dimensions.size(); i++) cout << " dimension " << i + 1 << ": " << dimensions[i] << endl;
299  cout << " This does not match a G4Cons. Exiting" << endl << endl;
300  exit(0);
301  }
302 
303  SolidV = new G4Cons(name,
304  dimensions[0],
305  dimensions[1],
306  dimensions[2],
307  dimensions[3],
308  dimensions[4],
309  dimensions[5],
310  dimensions[6]);
311 
312  built = 1;
313  }
314 
315  // #####
316  // Torus
317  // #####
318  if(type == "Torus")
319  {
320  if(dimensions.size() != 5)
321  {
322  cout << hd_msg << " Fatal Error: the number of dimensions for " << name
323  << " is " << dimensions.size() << ":" << endl;
324  for(unsigned int i=0; i<dimensions.size(); i++) cout << " dimension " << i + 1 << ": " << dimensions[i] << endl;
325  cout << " This does not match a G4Torus. Exiting" << endl << endl;
326  exit(0);
327  }
328  SolidV = new G4Torus(name,
329  dimensions[0],
330  dimensions[1],
331  dimensions[2],
332  dimensions[3],
333  dimensions[4]);
334 
335  built = 1;
336  }
337 
338 
339 
340  // #########
341  // Trapezoid
342  // #########
343  if(type == "Trd")
344  {
345  if(dimensions.size() != 5)
346  {
347  cout << hd_msg << " Fatal Error: the number of dimensions for " << name
348  << " is " << dimensions.size() << ":" << endl;
349  for(unsigned int i=0; i<dimensions.size(); i++) cout << " dimension " << i + 1 << ": " << dimensions[i] << endl;
350  cout << " This does not match a G4Trd. Exiting" << endl << endl;
351  exit(0);
352  }
353  SolidV = new G4Trd(name,
354  dimensions[0],
355  dimensions[1],
356  dimensions[2],
357  dimensions[3],
358  dimensions[4]);
359  built = 1;
360  }
361 
362 
363  // ##################
364  // Inclined Trapezoid
365  // ##################
366  if(type == "ITrd")
367  {
368  if(dimensions.size() != 7)
369  {
370  cout << hd_msg << " Fatal Error: the number of dimensions for " << name
371  << " is " << dimensions.size() << ":" << endl;
372  for(unsigned int i=0; i<dimensions.size(); i++) cout << " dimension " << i + 1 << ": " << dimensions[i] << endl;
373  cout << " This does not match a ITrd. Exiting" << endl << endl;
374  exit(0);
375  }
376 
377  double alph_xz = dimensions[5];
378  double alph_yz = dimensions[6];
379  double x = tan(alph_xz);
380  double y = tan(alph_yz);
381  double r = sqrt(x*x + y*y);
382 
383  // Calculating the 11 G4Trap parameters
384  double pDz = dimensions[4];
385  double pTheta = atan2(r, 1);
386  double pPhi = atan2(y, x);
387  double pDy1 = dimensions[2];
388  double pDx1 = dimensions[0];
389  double pDx2 = pDx1;
390  double pAlp1 = 0;
391  double pDy2 = dimensions[3];
392  double pDx3 = dimensions[1];
393  double pDx4 = pDx3;
394  double pAlp2 = 0;
395 
396  SolidV = new G4Trap(name,
397  pDz,
398  pTheta,
399  pPhi,
400  pDy1,
401  pDx1,
402  pDx2,
403  pAlp1,
404  pDy2,
405  pDx3,
406  pDx4,
407  pAlp2);
408 
409  built = 1;
410  }
411 
412 
413  // ########################
414  // Geant4 Generic Trapezoid
415  // ########################
416  if(type == "G4Trap")
417  {
418  if(dimensions.size() != 11)
419  {
420  cout << hd_msg << " Fatal Error: the number of dimensions for " << name
421  << " is " << dimensions.size() << ":" << endl;
422  for(unsigned int i=0; i<dimensions.size(); i++) cout << " dimension " << i + 1 << ": " << dimensions[i] << endl;
423  cout << " This does not match a G4Trd. Exiting" << endl << endl;
424  exit(0);
425  }
426  SolidV = new G4Trap(name,
427  dimensions[0],
428  dimensions[1],
429  dimensions[2],
430  dimensions[3],
431  dimensions[4],
432  dimensions[5],
433  dimensions[6],
434  dimensions[7],
435  dimensions[8],
436  dimensions[9],
437  dimensions[10]);
438  built = 1;
439  }
440 
441  // ##########################
442  // Geant4 Arbitrary Trapezoid
443  // ##########################
444  if(type == "G4GenericTrap")
445  {
446 
447  int nz = (dimensions.size() -1 ) / 2;
448  double dz = dimensions[0];
449  vector<G4TwoVector> vertices;
450 
451  for(int v=0; v<nz; v++)
452  vertices.push_back(G4TwoVector(dimensions[2*v+1], dimensions[2*v+2]));
453 
454  SolidV = new G4GenericTrap(name,
455  dz,
456  vertices);
457  built = 1;
458  }
459 
460  // #######################################
461  // G4Trap Constructor with 4+4 coordinates
462  // #######################################
463  if(type == "G4TrapC")
464  {
465  if(dimensions.size() != 24)
466  {
467  cout << hd_msg << " Fatal Error: the number of dimensions for " << name
468  << " is " << dimensions.size() << ":" << endl;
469  for(unsigned int i=0; i<dimensions.size(); i++) cout << " dimension " << i + 1 << ": " << dimensions[i] << endl;
470  cout << " This does not match a G4Trd. Exiting" << endl << endl;
471  exit(0);
472  }
473 
474  G4ThreeVector points[8];
475  for(int v=0; v<8; v++)
476  {
477  points[v].setX(dimensions[v*3+0]/cm * cm);
478  points[v].setY(dimensions[v*3+1]/cm * cm);
479  points[v].setZ(dimensions[v*3+2]/cm * cm);
480  }
481 
482  SolidV = new G4Trap(name, // name
483  points); // coordinates
484 
485  built = 1;
486  }
487 
488 
489 
490  // ####################
491  // Polyhedra (PGON)
492  // First G4 constructor
493  // WARNING: The order is NOT the same
494  // as the geant4 constructor.
495  // Let's revert to something that is compatible?
496  // ####################
497  if(type == "Pgon")
498  {
499  if(dimensions.size() < 8)
500  {
501  cout << hd_msg << " Fatal Error: the number of dimensions for " << name
502  << " is " << dimensions.size() << ":" << endl;
503  for(unsigned int i=0; i<dimensions.size(); i++) cout << " dimension " << i + 1 << ": " << dimensions[i] << endl;
504  cout << " This does not match a G4Polyhedra. Exiting" << endl << endl;
505  exit(0);
506  }
507  int nsides = (int) dimensions[2];
508  int zplanes = (int) dimensions[3];
509  if(nsides < 1)
510  {
511  cout << hd_msg << " Fatal Error: no sides for " << name
512  << "... should be a G4Polyhedra. Exiting" << endl << endl;
513  exit(0);
514  }
515  double* zPlane = new double[zplanes];
516  double* rInner = new double[zplanes];
517  double* rOuter = new double[zplanes];
518 
519  for(int zpl=0; zpl<zplanes; zpl++)
520  {
521  rInner[zpl] = dimensions[4 + 0*zplanes + zpl] ;
522  rOuter[zpl] = dimensions[4 + 1*zplanes + zpl] ;
523  zPlane[zpl] = dimensions[4 + 2*zplanes + zpl] ;
524  }
525  SolidV = new G4Polyhedra(name,
526  dimensions[0],
527  dimensions[1],
528  nsides,
529  zplanes,
530  zPlane,
531  rInner,
532  rOuter);
533  built = 1;
534  }
535 
536 
537  // ####################
538  // Polyhedra (PCON)
539  // First G4 constructor
540  // ####################
541  if(type == "Polycone")
542  {
543  if(dimensions.size() < 7)
544  {
545  cout << hd_msg << " Fatal Error: the number of dimensions for " << name
546  << " is " << dimensions.size() << ":" << endl;
547  for(unsigned int i=0; i<dimensions.size(); i++) cout << " dimension " << i + 1 << ": " << dimensions[i] << endl;
548  cout << " This does not match a G4Polycone. Exiting" << endl << endl;
549  exit(0);
550  }
551  int zplanes = (int) dimensions[2];
552 
553  double* zPlane = new double[zplanes];
554  double* rInner = new double[zplanes];
555  double* rOuter = new double[zplanes];
556 
557  for(int zpl=0; zpl<zplanes; zpl++)
558  {
559  // notice that the order should be different
560  // why didn't I put z first in the constructor as in geant4?
561  // might have to change it later
562  rInner[zpl] = dimensions[3 + 0*zplanes + zpl] ;
563  rOuter[zpl] = dimensions[3 + 1*zplanes + zpl] ;
564  zPlane[zpl] = dimensions[3 + 2*zplanes + zpl] ;
565  }
566  SolidV = new G4Polycone(name,
567  dimensions[0],
568  dimensions[1],
569  zplanes,
570  zPlane,
571  rInner,
572  rOuter);
573  built = 1;
574  }
575 
576 
577  // ############################
578  // CopyPlacement:
579  // Point LogicV to the original
580  // ############################
581  if(type.find("CopyOf") != string::npos && type.find("CopyOf") == 0)
582  {
583  hd_msg = gemcOpt.optMap["LOG_MSG"].args + " Copy: >> ";
584  string original(type, 6, 190);
585 
586  // Look for original
587  map<string, detector>::iterator it = (*Map).find(TrimSpaces(original));
588  if(it == (*Map).end())
589  {
590  cout << hd_msg << " <" << original << "> not found. Exiting." << endl << endl;
591  exit(0);
592  }
593  else
594  {
595  if(VERB>4 || name.find(catch_v) != string::npos)
596  {
597  cout << hd_msg << " " << name << " is a copy of <" << TrimSpaces(original) << ">. Pointing to its logical volume." << endl;
598  }
599  SetLogical(it->second.GetLogical());
600  }
601  built = 1;
602  }
603 
604 
605  // ################
606  // Solid Operations
607  // ################
608  if(type.find("Operation:") != string::npos && type.find("Operation:") == 0)
609  {
610  hd_msg = gemcOpt.optMap["LOG_MSG"].args + " Operation: >> ";
611 
612  bool translationFirst = false;
613 
614  // If Operation:~ it will perform the translation first
615  size_t posTld = type.find("~");
616  if( posTld != string::npos )
617  {
618  translationFirst = true;
619  type.replace( posTld, 1, " " );
620  }
621 
623  bool absolutecoordinates = false;
624 
625  // If Operation:@ it will assume that position of second object is given in the common mother volume of both objects
626 
627  size_t pos_at = type.find("@");
628  if( pos_at != string::npos )
629  {
630  absolutecoordinates = true;
631  type.replace( pos_at, 1, " " );
632  }
634 
635 
636  size_t posp = 0;
637  size_t posm = 0;
638  size_t post = 0;
639  size_t pos = 0;
640  // harcoded max size of 200 here
641  string operation(type, 10, 190);
642  posp = operation.find("+");
643  posm = operation.find("-");
644  post = operation.find("*");
645  if (posp != string::npos) pos = posp;
646  else if(posm != string::npos) pos = posm;
647  else if(post != string::npos) pos = post;
648  if(!posp && !posm && !post)
649  {
650  cout << hd_msg << " Operation " << operation << " for " << name << " not recognized. Exiting." << endl;
651  exit(0);
652  }
653 
654  // Locating solids
655  string solid1, solid2;
656  string tsolid1, tsolid2;
657  solid1.assign(operation, 0, pos);
658  solid2.assign(operation, pos+1, operation.size());
659  tsolid1 = TrimSpaces(solid1);
660  tsolid2 = TrimSpaces(solid2);
661 
662  // Locating second solid transformation
663  map<string, detector>::iterator it1 = (*Map).find(tsolid1);
664  map<string, detector>::iterator it2 = (*Map).find(tsolid2);
665  if(it1 == (*Map).end())
666  {
667  cout << hd_msg << " " << tsolid1 << " Not found. Exiting." << endl << endl;
668  exit(0);
669  }
670  if(it2 == (*Map).end())
671  {
672  cout << hd_msg << " " << tsolid2 << " Not found. Exiting." << endl << endl;
673  exit(0);
674  }
675 
676  // Define rotational and translational transformations then combine them
677  G4RotationMatrix rotate = it2->second.rot ;
678  G4ThreeVector translate = it2->second.pos;
679  G4RotationMatrix invRot = rotate.invert() ;
680  G4Transform3D transf1( invRot, G4ThreeVector( 0, 0, 0 ) );
681  G4Transform3D transf2( G4RotationMatrix(), translate );
682  G4Transform3D transform = transf2 * transf1 ;
683 
684  if( absolutecoordinates && TrimSpaces(it1->second.mother) == TrimSpaces(it2->second.mother) )
685  {
686  //assume that second object position and rotation are given in absolute (mother) coordinates:
687 
688  G4RotationMatrix invrot1 = (it1->second.rot).inverse();
689  G4RotationMatrix rotate2 = it2->second.rot;
690 
691  G4ThreeVector net_translation = it2->second.pos - it1->second.pos;
692 
693  // The net rotation should be the INVERSE of the rotation of object 2 relative to object 1.
694  // rotate1/2 is the rotation of object 1/2 relative to their common mother.
695  // If R is the rotation of 2 relative to 1, then R2 = R * R1 --> R = R2 R1^-1
696  G4RotationMatrix invnet_rotation = (rotate2 * invrot1).invert();
697 
698  // In order to express the relative position of object 2 in the coordinate system of object one, we must rotate it
699  // applying the same rotation as that used to position object 1, according to the GEANT4 framework:
700  // I do not quite understand WHY this works, but through trial and error, I have
701  // discovered that the combination of operations below is what works:
702  net_translation *= it1->second.rot;
703  transform = G4Transform3D( invnet_rotation, net_translation );
704  //We don't want there to be any possibility to overwrite "transform" in this special case, so we force translationFirst to false here:
705  translationFirst = false;
706  }
707 
708 
709  // If there was tilda in the operation string then the rotation and translation are switched
710  // with respect to the default behaviour of G4UnionSolid with separate rotational matrix
711  // and translatin vector
712  if( translationFirst )
713  {
714  transform = transf1 * transf2 ;
715  }
716 
717  if(posp != string::npos)
718  {
719  SolidV = new G4UnionSolid( name, it1->second.GetSolid(), it2->second.GetSolid(), transform );
720  }
721  if(posm != string::npos)
722  {
723  SolidV = new G4SubtractionSolid( name, it1->second.GetSolid(), it2->second.GetSolid(), transform );
724  }
725  if(post != string::npos)
726  {
727  SolidV = new G4IntersectionSolid(name, it1->second.GetSolid(), it2->second.GetSolid(), transform );
728  }
729 
730  if(VERB>4 || name.find(catch_v) != string::npos)
731  {
732  cout << hd_msg << " " << name << " is the " << (pos==posp ? " sum " : " difference ") << " of " << tsolid1 << " and " << tsolid2 << endl;;
733  }
734  built = 1;
735  }
736 
737 
738 
739 
740  if(VERB>4 || name.find(catch_v) != string::npos)
741  {
742  cout << hd_msg << " " << name << " solid " << type << " built." << endl;
743  }
744 
745 
746  if(built==0)
747  {
748  cout << hd_msg << " " << name << " solid >" << type << "< not recognized. Exiting." << endl;
749  exit(0);
750  }
751  return 1;
752 }
753 
754 
755 int detector::create_logical_volume(map<string, G4Material*> *MMats, goptions gemcOpt)
756 {
757  string hd_msg = gemcOpt.optMap["LOG_MSG"].args + " Logical: >> ";
758  double VERB = gemcOpt.optMap["G4P_VERBOSITY"].arg ;
759  string catch_v = gemcOpt.optMap["CATCH"].args;
760  string defmat = gemcOpt.optMap["DEFAULT_MATERIAL"].args;
761 
762  // don't build the logical volumes for components or replicas
763  if(material == "Component" || material == "OfReplica")
764  {
765  if(VERB>4 || name.find(catch_v) != string::npos)
766  cout << hd_msg << " " << name << " is a Solid Component or a Replicant. Logical Volume will not be built." << endl;
767  return 0;
768  }
769 
770  // Check if Material Exists
771  map<string, G4Material*>::iterator i = MMats->find(material);
772 
773  // if material is not defined, look in G4 table.
774  if(i == MMats->end() && LogicV == 0)
775  {
776  G4NistManager* matman = G4NistManager::Instance();
777  if(matman->FindOrBuildMaterial(material)) (*MMats)[material] = matman->FindOrBuildMaterial(material);
778  }
779  i = MMats->find(material);
780 
781  // if material is still not defined, use air
782  if(i == MMats->end() && LogicV == 0)
783  {
784  if(defmat == "none")
785  {
786  cout << hd_msg << " Warning: material >" << material << "< is not defined. Exiting" << endl;
787  cout << hd_msg << " You can set the DEFAULT_MATERIAL flag to replace an undefined material. " << endl;
788  exit(0);
789  }
790  else
791  {
792  material = defmat;
793  if(MMats->find(material)== MMats->end())
794  {
795  cout << hd_msg << " Warning: " << defmat << " set with DEFAULT_MATERIAL is not found. Exiting" << endl;
796  exit(0);
797 
798  }
799  }
800  }
801 
802  // Logical Volume Basic Constructor
803  // If LogicV exists already, this is a copy
804  if(LogicV == 0)
805  LogicV = new G4LogicalVolume(SolidV, (*MMats)[material], name, 0, 0, 0, true);
806 
807  if(name == "root") LogicV->SetVisAttributes(G4VisAttributes::GetInvisible());
808  else LogicV->SetVisAttributes(VAtts);
809 
810  if(VERB>4 || name.find(catch_v) != string::npos)
811  {
812  cout << hd_msg << " " << name << " Logical Volume built." << endl;
813  }
814 
815  return 1;
816 }
817 
818 
819 int detector::create_physical_volumes(goptions gemcOpt, G4LogicalVolume *mamma)
820 {
821  string hd_msg = gemcOpt.optMap["LOG_MSG"].args + " Physical: >> ";
822  double VERB = gemcOpt.optMap["G4P_VERBOSITY"].arg ;
823  bool OVERL = gemcOpt.optMap["CHECK_OVERLAPS"].arg > 0 ;
824  string catch_v = gemcOpt.optMap["CATCH"].args;
825  if(PhysicalV) delete PhysicalV;
826 
827  // don't build physical volumes for components or replicas.
828  // Replicas are built in the dedicated routine
829  if(material == "Component" || material == "OfReplica")
830  {
831  if(VERB>4 || name.find(catch_v) != string::npos)
832  cout << hd_msg << " " << name << " is a Solid Component or a Replicant. Physical Volume will not be built, or replicas will be built instead." << endl;
833  return 1;
834  }
835 
836 
837  if(name == "root")
838  PhysicalV = new G4PVPlacement(0,
839  G4ThreeVector(),
840  LogicV,
841  name.c_str(),
842  0,
843  false,
844  0);
845 
846  else
847  PhysicalV = new G4PVPlacement(&rot,
848  pos,
849  LogicV,
850  name.c_str(),
851  mamma,
852  false,
853  ncopy,
854  OVERL);
855 
856  if(VERB>4 || name.find(catch_v) != string::npos)
857  {
858  if(mamma)
859  cout << hd_msg << " " << name << " Physical Volume(s) built inside " << mamma->GetName() << "." << endl;
860  }
861 
862  return 1;
863 }
864 
865 int detector::create_replicas(goptions gemcOpt, G4LogicalVolume *mamma, detector replicant)
866 {
867  string hd_msg = gemcOpt.optMap["LOG_MSG"].args + " Physical: >> ";
868  double VERB = gemcOpt.optMap["G4P_VERBOSITY"].arg ;
869  string catch_v = gemcOpt.optMap["CATCH"].args;
870 
871  // deleting the replicant detector physicsal volume as well
872  if(replicant.PhysicalV) delete replicant.PhysicalV;
873 
874  // reset this physical volume
875  if(PhysicalV) delete PhysicalV;
876 
877 
878  if(type.find("ReplicaOf:") == 0)
879  {
880  if(dimensions.size() != 4)
881  {
882  cout << hd_msg << " Fatal Error: the number of parameters for " << name
883  << " is " << dimensions.size() << ":" << endl;
884  for(unsigned int i=0; i<dimensions.size(); i++)
885  cout << " parameter " << i + 1 << ": " << dimensions[i] << endl;
886  cout << " This does not match a G4Replicas (4). Exiting" << endl << endl;
887  exit(0);
888  }
889 
890  EAxis pAxis;
891  if(dimensions[0] == 1) pAxis = kXAxis;
892  if(dimensions[0] == 2) pAxis = kYAxis;
893  if(dimensions[0] == 3) pAxis = kZAxis;
894  int nreps = (int) get_number(dimensions[1]);
895  double width = get_number(dimensions[2]);
896  double offset = get_number(dimensions[3]);
897 
898  // the logical volume is built
899  // the mother is built as well
900  PhysicalV = new G4PVReplica(name,
901  replicant.LogicV,
902  mamma,
903  pAxis,
904  nreps,
905  width,
906  offset);
907 
908  }
909 
910  if(VERB>4 || name.find(catch_v) != string::npos)
911  {
912  if(mamma)
913  cout << hd_msg << " " << name << " Physical Volume(s) built inside " << mamma->GetName() << "." << endl;
914  }
915 
916  return 1;
917 }
918 
919 
920 
921 // Returns dimensions nomenclature for different solid type
922 vector< vector<string> > dimensionstype(string solidtype)
923 {
924  vector< vector<string> > dtypes;
925  vector<string> dt;
926  dt.resize(2);
927 
928  if(solidtype == "Box")
929  {
930  dt[0] = "half length in X";
931  dt[1] = "Length";
932  dtypes.push_back(dt);
933  dt[0] = "half length in Y";
934  dtypes.push_back(dt);
935  dt[0] = "half length in Z";
936  dtypes.push_back(dt);
937  }
938 
939  if(solidtype == "Sphere")
940  {
941  dt[1] = "Length";
942  dt[0] = "Inner radius";
943  dtypes.push_back(dt);
944  dt[0] = "Outer radius";
945  dtypes.push_back(dt);
946  dt[1] = "Angle";
947  dt[0] = "Starting Phi angle of the segment";
948  dtypes.push_back(dt);
949  dt[0] = "Delta Phi angle of the segment";
950  dtypes.push_back(dt);
951  dt[0] = "Starting Theta angle of the segment";
952  dtypes.push_back(dt);
953  dt[0] = "Delta Theta angle of the segment";
954  dtypes.push_back(dt);
955  }
956  if(solidtype == "Tube")
957  {
958  dt[1] = "Length";
959  dt[0] = "Inner radius";
960  dtypes.push_back(dt);
961  dt[0] = "Outer radius";
962  dtypes.push_back(dt);
963  dt[0] = "Half length in z";
964  dtypes.push_back(dt);
965  dt[1] = "Angle";
966  dt[0] = "Starting Phi angle";
967  dtypes.push_back(dt);
968  dt[0] = "Delta Phi angle";
969  dtypes.push_back(dt);
970  }
971  if(solidtype == "Trd")
972  {
973  dt[1] = "Length";
974  dt[0] = "Half-length along x at the surface at -dz";
975  dtypes.push_back(dt);
976  dt[0] = "Half-length along x at the surface at +dz";
977  dtypes.push_back(dt);
978  dt[0] = "Half-length along y at the surface at -dz";
979  dtypes.push_back(dt);
980  dt[0] = "Half-length along y at the surface at +dz";
981  dtypes.push_back(dt);
982  dt[0] = "dz: Half-length along z axis";
983  dtypes.push_back(dt);
984  }
985  if(solidtype == "Cons")
986  {
987  dt[1] = "Length";
988  dt[0] = "Inner radius at start";
989  dtypes.push_back(dt);
990  dt[0] = "Outer radius at start";
991  dtypes.push_back(dt);
992  dt[0] = "Inner radius at end";
993  dtypes.push_back(dt);
994  dt[0] = "Outer radius at end";
995  dtypes.push_back(dt);
996  dt[0] = "Half length in z";
997  dtypes.push_back(dt);
998  dt[1] = "Angle";
999  dt[0] = "Starting Phi angle";
1000  dtypes.push_back(dt);
1001  dt[0] = "Delta Phi angle";
1002  dtypes.push_back(dt);
1003  }
1004  if(solidtype == "G4Trap")
1005  {
1006  dt[1] = "Length";
1007  dt[0] = "Half z length ";
1008  dtypes.push_back(dt);
1009  dt[1] = "Angle";
1010  dt[0] = "Polar angle of the line joining the centres of the faces";
1011  dtypes.push_back(dt);
1012  dt[0] = "Azimuthal angle of the line joining the centre of the face";
1013  dtypes.push_back(dt);
1014  dt[1] = "Length";
1015  dt[0] = "Half y length at -pDz";
1016  dtypes.push_back(dt);
1017  dt[0] = "Half x length of the side at y=-pDy1, -pDz";
1018  dtypes.push_back(dt);
1019  dt[0] = "Half x length of the side at y=+pDy1, -pDz";
1020  dtypes.push_back(dt);
1021  dt[1] = "Angle";
1022  dt[0] = "Angle to the y axis from the centre (lower endcap) ";
1023  dtypes.push_back(dt);
1024  dt[1] = "Length";
1025  dt[0] = "Half y length at +pDz";
1026  dtypes.push_back(dt);
1027  dt[0] = "Half x length of the side at y=-pDy2, +pDz";
1028  dtypes.push_back(dt);
1029  dt[0] = "Half x length of the side at y=+pDy2, +pDz";
1030  dtypes.push_back(dt);
1031  dt[1] = "Angle";
1032  dt[0] = "Angle to the y axis from the centre (upper endcap) ";
1033  dtypes.push_back(dt);
1034  }
1035 
1036  if(solidtype == "G4EllipticalTube")
1037  {
1038  dt[1]="Lenght";
1039  dt[0]="Half length x";
1040  dtypes.push_back(dt);
1041  dt[1]="Lenght";
1042  dt[0]="Half length y";
1043  dtypes.push_back(dt);
1044  dt[1]="Lenght";
1045  dt[0]="Half length z";
1046  dtypes.push_back(dt);
1047  }
1048 
1049  if(solidtype == "Hype")
1050  {
1051  dt[1]="Length";
1052  dt[0]="Inner radius";
1053  dtypes.push_back(dt);
1054  dt[0] = "Outer radius";
1055  dtypes.push_back(dt);
1056  dt[1] = "Angle";
1057  dt[0] = "Inner stereo angle";
1058  dtypes.push_back(dt);
1059  dt[0] = "Outer stereo angle";
1060  dtypes.push_back(dt);
1061  dt[1]="Length";
1062  dt[0] = "Half length in z";
1063  dtypes.push_back(dt);
1064  }
1065 
1066  if(solidtype == "Parallelepiped")
1067  {
1068  dt[1]="Length";
1069  dt[0]="Half length in x";
1070  dtypes.push_back(dt);
1071  dt[0] = "Half length in y";
1072  dtypes.push_back(dt);
1073  dt[0] = "Half length in z";
1074  dtypes.push_back(dt);
1075  dt[1] = "Angle";
1076  dt[0] = "Angle formed by the y axis and the plane joining the centre of the faces parallel to the z-x plane at -dy and +dy" ;
1077  dtypes.push_back(dt);
1078  dt[0] = "Polar angle of the line joining the centres of the faces at -dz and +dz in z ";
1079  dtypes.push_back(dt);
1080  dt[0] = "Azimuthal angle of the line joining the centres of the faces at -dz and +dz in z ";
1081  dtypes.push_back(dt);
1082  }
1083 
1084  if(solidtype == "Torus")
1085  {
1086  dt[1]="Length";
1087  dt[0]="Inner radius";
1088  dtypes.push_back(dt);
1089  dt[0] = "Outer radius";
1090  dtypes.push_back(dt);
1091  dt[0] = "Swept radius of torus";
1092  dtypes.push_back(dt);
1093  dt[1] = "Angle";
1094  dt[0] = "Starting Phi angle";
1095  dtypes.push_back(dt);
1096  dt[0] = "Outer stereo angle";
1097  dtypes.push_back(dt);
1098  }
1099 
1100  if(solidtype == "Ellipsoid")
1101  {
1102  dt[1]="Length";
1103  dt[0]="Semiaxis in X";
1104  dtypes.push_back(dt);
1105  dt[0] = "Semiaxis in Y ";
1106  dtypes.push_back(dt);
1107  dt[0] = "Semiaxis in Z ";
1108  dtypes.push_back(dt);
1109  dt[0] = "lower cut plane level, z";
1110  dtypes.push_back(dt);
1111  dt[0] = "upper cut plane level, z";
1112  dtypes.push_back(dt);
1113  }
1114 
1115  return dtypes;
1116 }
1117 
1118 
1119 ostream &operator<<(ostream &stream, detector Detector)
1120 {
1121  cout << endl;
1122  cout << " Detector name: " << Detector.name << " - " << Detector.description << endl;
1123  cout << " Mother: " << Detector.mother << endl;
1124  cout << " Position (cm): " << Detector.pos/cm << endl;
1125  cout << " Rotation: " << Detector.rot << endl;
1126  cout << " Color: " << Detector.VAtts.GetColour() << endl;
1127  cout << " Type: " << Detector.type << endl;
1128  vector< vector<string> > dtypes = dimensionstype( Detector.type.c_str() );
1129 
1130  if(dtypes.size() != Detector.dimensions.size() && Detector.type.find("CopyOf") != 0)
1131  {
1132  for(unsigned int i=0; i<Detector.dimensions.size(); i++)
1133  cout << " Size " << i + 1 << ": " << Detector.dimensions[i] << endl;
1134  }
1135  if(dtypes.size() == Detector.dimensions.size() && Detector.type.find("CopyOf") != 0)
1136  for(unsigned int i=0; i<dtypes.size(); i++)
1137  cout << " " << dtypes[i][0] << ": " << G4BestUnit(Detector.dimensions[i], dtypes[i][1]) << endl;
1138 
1139  cout << " Material: " << Detector.material << endl;
1140  cout << " Magnetic Field: " << Detector.magfield << endl;
1141  cout << " Copy Number: " << Detector.ncopy << endl;
1142  cout << " Activated: " << ( Detector.exist==1 ? "yes" : "no" ) << endl;
1143  cout << " Visible: " << ( Detector.visible==1 ? "yes" : "no" ) << endl;
1144  cout << " Style: " << ( Detector.style==1 ? "solid" : "wireframe" ) << endl;
1145  cout << " Sensitivity: " << Detector.sensitivity << endl;
1146  if(Detector.sensitivity != "no")
1147  cout << " hitType: " << Detector.hitType << endl;
1148 
1149  if(Detector.identity.size())
1150  cout << Detector.identity ;
1151 
1152  cout << endl;
1153 
1154  return stream;
1155 }
1156 
1157 bool detector::operator == (const detector& D) const
1158 {
1159  // Name uniquely identifies a volume
1160  if(D.name == this->name)
1161  return true;
1162  else return false;
1163 }
1164 
1165 
1166 
1167 
1168 
1169 
1170 
1171 
1172 
1173 
1174 
1175 
1176 
1177 
1178 
1179 
1180 
1181 
1182 
1183 
1184 
1185 
1186 
1187 
1188 
1189 
detector::create_logical_volume
int create_logical_volume(map< string, G4Material * > *, goptions)
Creates the Logical Volume.
Definition: detector.cc:755
detector::sensitivity
string sensitivity
Defines the Sensitive Detector. possible choices: "no" "hits collection name".
Definition: detector.h:89
detector::mother
string mother
Mother Volume name.
Definition: detector.h:68
operator<<
ostream & operator<<(ostream &stream, detector Detector)
Definition: detector.cc:1119
std
STL namespace.
detector::pos
G4ThreeVector pos
Position relative to the mother volume, as G4ThreeVector.
Definition: detector.h:71
detector::create_physical_volumes
int create_physical_volumes(goptions, G4LogicalVolume *)
Creates the Physical Volume.
Definition: detector.cc:819
detector::operator==
bool operator==(const detector &D) const
Overloaded "==" operator for detector class.
Definition: detector.cc:1157
detector::VAtts
G4VisAttributes VAtts
Visual Attributes: color, transparency, style (wireframe, solid), visibility.
Definition: detector.h:74
get_number
double get_number(string v, int warn_no_unit)
Returns number with dimension from string, i.e. 100*cm.
Definition: string_utilities.cc:134
detector::create_solid
int create_solid(goptions, map< string, detector > *)
Creates the Solid. If it&#39;s a Copy Placement, retrieve and assigns LogicV.
Definition: detector.cc:47
dimensionstype
vector< vector< string > > dimensionstype(string solidtype)
Returns dimensions nomenclature for different solid type.
Definition: detector.cc:922
detector::ncopy
int ncopy
copy number
Definition: detector.h:82
detector::name
string name
Name of the volume. Since this is the key of the STL map, it has to be unique.
Definition: detector.h:67
goptions::optMap
map< string, aopt > optMap
Options map.
Definition: options.h:75
detector::material
string material
Volume Material name.
Definition: detector.h:79
detector::detector
detector()
Definition: detector.cc:40
detector::visible
int visible
visibility of the detector: 0=invisible 1=visible
Definition: detector.h:86
detector::style
int style
Visual style: 0=wireframe 1=solid.
Definition: detector.h:87
detector::hitType
string hitType
Hit Process routine name. A Hit Process HitProcess derived class must exists with this name...
Definition: detector.h:90
detector::identity
vector< identifier > identity
Vector of identifiers. Example: superlayer manual 1 type manual 2 segment manual 3 strip manual 4...
Definition: detector.h:91
detector::magfield
string magfield
Magnetic Field. The string "no" means that the field is inherited from the mother volume...
Definition: detector.h:80
detector::rot
G4RotationMatrix rot
Rotation Matrix, defined by rotations along x,y,z axis relative to the mother volume.
Definition: detector.h:72
detector::exist
int exist
detector ON/OFF switch
Definition: detector.h:85
detector::type
string type
solid type. This follows the GEANT4 definitions
Definition: detector.h:76
detector::dimensions
vector< double > dimensions
vector of dimensions. Size, units depends on solid type
Definition: detector.h:77
TrimSpaces
string TrimSpaces(string in)
Removes leading and trailing spaces.
Definition: string_utilities.cc:31
detector::create_replicas
int create_replicas(goptions, G4LogicalVolume *, detector)
Creates the Replica Volumes.
Definition: detector.cc:865
detector::description
string description
Volume Description for documentation.
Definition: detector.h:69
Generated on Wed Jul 6 2016 11:13:41 for GEMC by   doxygen 1.8.11