| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4ParameterisationPolyhedra[Rho/Phi/Z] impl
27 //
28 // 14.10.03 - P.Arce, Initial version
29 // 08.04.04 - I.Hrivnacova, Implemented reflec
30 // -------------------------------------------
31
32 #include "G4ParameterisationPolyhedra.hh"
33
34 #include <iomanip>
35 #include "G4PhysicalConstants.hh"
36 #include "G4ThreeVector.hh"
37 #include "G4GeometryTolerance.hh"
38 #include "G4RotationMatrix.hh"
39 #include "G4VPhysicalVolume.hh"
40 #include "G4LogicalVolume.hh"
41 #include "G4ReflectedSolid.hh"
42 #include "G4Polyhedra.hh"
43
44 //--------------------------------------------
45 G4VParameterisationPolyhedra::
46 G4VParameterisationPolyhedra( EAxis axis, G4in
47 G4double offset,
48 DivisionType div
49 : G4VDivisionParameterisation( axis, nDiv,
50 {
51 std::ostringstream message;
52 /* #ifdef G4MULTITHREADED
53 message << "Divisions for G4Polyhedra curren
54 << G4endl
55 << "Sorry! Solid: " << msolid->GetNa
56 G4Exception("G4VParameterisationPolyhedra::G
57 "GeomDiv0001", FatalException, m
58 #endif */
59
60 auto msol = (G4Polyhedra*)(msolid);
61 if ((msolid->GetEntityType() != "G4Reflected
62 {
63 message << "Generic construct for G4Polyhe
64 << "Sorry! Solid: " << msol->GetNa
65 G4Exception("G4VParameterisationPolyhedra:
66 "GeomDiv0001", FatalException,
67 }
68 if (msolid->GetEntityType() == "G4ReflectedS
69 {
70 // Get constituent solid
71 //
72 G4VSolid* mConstituentSolid
73 = ((G4ReflectedSolid*)msolid)->GetConst
74 msol = (G4Polyhedra*)(mConstituentSolid);
75
76 // Get parameters
77 //
78 G4int nofSides = msol->GetOriginalParamete
79 G4int nofZplanes = msol->GetOriginalParame
80 G4double* zValues = msol->GetOriginalParam
81 G4double* rminValues = msol->GetOriginalPa
82 G4double* rmaxValues = msol->GetOriginalPa
83
84 // Invert z values, convert radius paramet
85 //
86 auto rminValues2 = new G4double[nofZplanes
87 auto rmaxValues2 = new G4double[nofZplanes
88 auto zValuesRefl = new G4double[nofZplanes
89 for (G4int i=0; i<nofZplanes; ++i)
90 {
91 rminValues2[i] = rminValues[i] * Convert
92 rmaxValues2[i] = rmaxValues[i] * Convert
93 zValuesRefl[i] = - zValues[i];
94 }
95
96 auto newSolid
97 = new G4Polyhedra(msol->GetName(),
98 msol->GetStartPhi(),
99 msol->GetEndPhi() - ms
100 nofSides,
101 nofZplanes, zValuesRef
102
103 delete [] rminValues2;
104 delete [] rmaxValues2;
105 delete [] zValuesRefl;
106
107 msol = newSolid;
108 fmotherSolid = newSolid;
109 fReflectedSolid = true;
110 fDeleteSolid = true;
111 }
112 }
113
114 //--------------------------------------------
115 G4VParameterisationPolyhedra::~G4VParameterisa
116
117 //--------------------------------------------
118 G4double
119 G4VParameterisationPolyhedra::
120 ConvertRadiusFactor(const G4Polyhedra& phedra)
121 {
122 G4double phiTotal = phedra.GetEndPhi() - phe
123 G4int nofSides = phedra.GetOriginalParameter
124
125 if ( (phiTotal <=0) || (phiTotal >
126 2*pi+G4GeometryTolerance::GetInstance(
127 { phiTotal = 2*pi; }
128
129 return std::cos(0.5*phiTotal/nofSides);
130 }
131
132 //--------------------------------------------
133 G4ParameterisationPolyhedraRho::
134 G4ParameterisationPolyhedraRho( EAxis axis, G4
135 G4double width,
136 G4VSolid* msoli
137 : G4VParameterisationPolyhedra( axis, nDiv,
138 {
139
140 CheckParametersValidity();
141 SetType( "DivisionPolyhedraRho" );
142
143 auto msol = (G4Polyhedra*)(fmotherSolid);
144 G4PolyhedraHistorical* original_pars = msol-
145
146 if( divType == DivWIDTH )
147 {
148 fnDiv = CalculateNDiv( original_pars->Rmax
149 - original_pars->Rmin
150 }
151 else if( divType == DivNDIV )
152 {
153 fwidth = CalculateWidth( original_pars->Rm
154 - original_pars->Rm
155 }
156
157 #ifdef G4DIVDEBUG
158 if( verbose >= 1 )
159 {
160 G4cout << " G4ParameterisationPolyhedraRho
161 << " = " << nDiv << G4endl
162 << " Offset " << foffset << " = " <
163 << " Width " << fwidth << " = " <<
164 }
165 #endif
166 }
167
168 //--------------------------------------------
169 G4ParameterisationPolyhedraRho::~G4Parameteris
170
171 //--------------------------------------------
172 void G4ParameterisationPolyhedraRho::CheckPara
173 {
174 G4VDivisionParameterisation::CheckParameters
175
176 auto msol = (G4Polyhedra*)(fmotherSolid);
177
178 if( fDivisionType == DivNDIVandWIDTH || fDiv
179 {
180 std::ostringstream message;
181 message << "In solid " << msol->GetName()
182 << "Division along R will be done
183 << "different for each solid secti
184 << "WIDTH will not be used !";
185 G4Exception("G4ParameterisationPolyhedraRh
186 "GeomDiv1001", JustWarning, me
187 }
188 if( foffset != 0. )
189 {
190 std::ostringstream message;
191 message << "In solid " << msol->GetName()
192 << "Division along R will be done
193 << "different for each solid secti
194 << "OFFSET will not be used !";
195 G4Exception("G4ParameterisationPolyhedraRh
196 "GeomDiv1001", JustWarning, me
197 }
198 }
199
200 //--------------------------------------------
201 G4double G4ParameterisationPolyhedraRho::GetMa
202 {
203 auto msol = (G4Polyhedra*)(fmotherSolid);
204 G4PolyhedraHistorical* original_pars = msol-
205 return original_pars->Rmax[0] - original_par
206 }
207
208 //--------------------------------------------
209 void
210 G4ParameterisationPolyhedraRho::
211 ComputeTransformation( const G4int, G4VPhysica
212 {
213 //----- translation
214 G4ThreeVector origin(0.,0.,0.);
215
216 //----- set translation
217 physVol->SetTranslation( origin );
218
219 //----- calculate rotation matrix: unit
220
221 #ifdef G4DIVDEBUG
222 if( verbose >= 2 )
223 {
224 G4cout << " G4ParameterisationPolyhedraRho
225 << " foffset: " << foffset/CLHEP::d
226 << " - fwidth: " << fwidth/CLHEP::d
227 }
228 #endif
229
230 ChangeRotMatrix( physVol );
231
232 #ifdef G4DIVDEBUG
233 if( verbose >= 2 )
234 {
235 G4cout << std::setprecision(8) << " G4Para
236 << G4endl
237 << " Position: " << origin
238 << " - Width: " << fwidth
239 << " - Axis: " << faxis << G4endl;
240 }
241 #endif
242 }
243
244 //--------------------------------------------
245 void
246 G4ParameterisationPolyhedraRho::
247 ComputeDimensions( G4Polyhedra& phedra, const
248 const G4VPhysicalVolume* )
249 {
250 auto msol = (G4Polyhedra*)(fmotherSolid);
251
252 G4PolyhedraHistorical* origparamMother = mso
253 G4PolyhedraHistorical origparam( *origparamM
254 G4int nZplanes = origparamMother->Num_z_plan
255
256 G4double width = 0.;
257 for( G4int ii = 0; ii < nZplanes; ++ii )
258 {
259 width = CalculateWidth( origparamMother->R
260 - origparamMother->R
261 origparam.Rmin[ii] = origparamMother->Rmin
262 origparam.Rmax[ii] = origparamMother->Rmin
263 }
264
265 phedra.SetOriginalParameters(&origparam); //
266 phedra.Reset(); //
267
268 #ifdef G4DIVDEBUG
269 if( verbose >= -2 )
270 {
271 G4cout << "G4ParameterisationPolyhedraRho:
272 << "-- Parametrised phedra copy-num
273 phedra.DumpInfo();
274 }
275 #endif
276 }
277
278 //--------------------------------------------
279 G4ParameterisationPolyhedraPhi::
280 G4ParameterisationPolyhedraPhi( EAxis axis, G4
281 G4double width,
282 G4VSolid* msoli
283 : G4VParameterisationPolyhedra( axis, nDiv,
284 {
285 CheckParametersValidity();
286 SetType( "DivisionPolyhedraPhi" );
287
288 auto msol = (G4Polyhedra*)(fmotherSolid);
289 G4double deltaPhi = msol->GetEndPhi() - msol
290
291 if( divType == DivWIDTH )
292 {
293 fnDiv = msol->GetNumSide();
294 }
295
296 fwidth = CalculateWidth( deltaPhi, fnDiv, 0.
297
298 #ifdef G4DIVDEBUG
299 if( verbose >= 1 )
300 {
301 G4cout << " G4ParameterisationPolyhedraPhi
302 << " = " << nDiv << G4endl
303 << " Offset " << foffset << " = " <
304 << " Width " << fwidth << " = " <<
305 }
306 #endif
307 }
308
309 //--------------------------------------------
310 G4ParameterisationPolyhedraPhi::~G4Parameteris
311
312 //--------------------------------------------
313 G4double G4ParameterisationPolyhedraPhi::GetMa
314 {
315 auto msol = (G4Polyhedra*)(fmotherSolid);
316 return msol->GetEndPhi() - msol->GetStartPhi
317 }
318
319 //--------------------------------------------
320 void G4ParameterisationPolyhedraPhi::CheckPara
321 {
322 G4VDivisionParameterisation::CheckParameters
323
324 auto msol = (G4Polyhedra*)(fmotherSolid);
325
326 if( fDivisionType == DivNDIVandWIDTH || fDiv
327 {
328 std::ostringstream message;
329 message << "In solid " << msol->GetName()
330 << " Division along PHI will be do
331 << "in the defined numSide." << G4
332 << "WIDTH will not be used !";
333 G4Exception("G4ParameterisationPolyhedraPh
334 "GeomDiv1001", JustWarning, me
335 }
336 if( foffset != 0. )
337 {
338 std::ostringstream message;
339 message << "In solid " << msol->GetName()
340 << "Division along PHI will be don
341 << "in the defined numSide." << G4
342 << "OFFSET will not be used !";
343 G4Exception("G4ParameterisationPolyhedraPh
344 "GeomDiv1001", JustWarning, me
345 }
346
347 G4PolyhedraHistorical* origparamMother = mso
348
349 if( origparamMother->numSide != fnDiv && fD
350 {
351 std::ostringstream message;
352 message << "Configuration not supported."
353 << "Division along PHI will be don
354 << G4endl
355 << "numSide, i.e, the number of di
356 << origparamMother->numSide << " i
357 G4Exception("G4ParameterisationPolyhedraPh
358 "GeomDiv0001", FatalException,
359 }
360 }
361
362 //--------------------------------------------
363 void
364 G4ParameterisationPolyhedraPhi::
365 ComputeTransformation( const G4int copyNo, G4V
366 {
367 //----- translation
368 G4ThreeVector origin(0.,0.,0.);
369 //----- set translation
370 physVol->SetTranslation( origin );
371
372 //----- calculate rotation matrix (so that a
373 G4double posi = copyNo*fwidth;
374
375 #ifdef G4DIVDEBUG
376 if( verbose >= 2 )
377 {
378 G4cout << " G4ParameterisationPolyhedraPhi
379 << G4endl
380 << " copyNo: " << copyNo
381 << " - fwidth: " << fwidth/CLHEP::d
382 }
383 #endif
384
385 ChangeRotMatrix( physVol, -posi );
386
387 #ifdef G4DIVDEBUG
388 if( verbose >= 2 )
389 {
390 G4cout << std::setprecision(8) << " G4Para
391 << copyNo << G4endl
392 << " Position: " << origin << " - W
393 << " - Axis: " << faxis << G4endl;
394 }
395 #endif
396 }
397
398 //--------------------------------------------
399 void
400 G4ParameterisationPolyhedraPhi::
401 ComputeDimensions( G4Polyhedra& phedra, const
402 const G4VPhysicalVolume* )
403 {
404 auto msol = (G4Polyhedra*)(fmotherSolid);
405
406 G4PolyhedraHistorical* origparamMother = mso
407 G4PolyhedraHistorical origparam( *origparamM
408
409 origparam.numSide = 1;
410 origparam.Start_angle = origparamMother->Sta
411 origparam.Opening_angle = fwidth;
412
413 phedra.SetOriginalParameters(&origparam); /
414 phedra.Reset(); /
415
416 #ifdef G4DIVDEBUG
417 if( verbose >= 2 )
418 {
419 G4cout << "G4ParameterisationPolyhedraPhi:
420 phedra.DumpInfo();
421 }
422 #endif
423 }
424
425 //--------------------------------------------
426 G4ParameterisationPolyhedraZ::
427 G4ParameterisationPolyhedraZ( EAxis axis, G4in
428 G4double width, G
429 G4VSolid* msolid,
430 : G4VParameterisationPolyhedra( axis, nDiv,
431 fOrigParamMother(((G4Polyhedra*)fmotherSo
432 {
433 CheckParametersValidity();
434 SetType( "DivisionPolyhedraZ" );
435
436 if( divType == DivWIDTH )
437 {
438 fnDiv =
439 CalculateNDiv(fOrigParamMother->Z_values
440 - fOrigParamMother->Z_valu
441 }
442 else if( divType == DivNDIV )
443 {
444 fwidth =
445 CalculateNDiv(fOrigParamMother->Z_values
446 - fOrigParamMother->Z_valu
447 }
448
449 #ifdef G4DIVDEBUG
450 if( verbose >= 1 )
451 {
452 G4cout << " G4ParameterisationPolyhedraZ -
453 << nDiv << G4endl
454 << " Offset " << foffset << " = " <
455 << " Width " << fwidth << " = " <<
456 }
457 #endif
458 }
459
460 //--------------------------------------------
461 G4ParameterisationPolyhedraZ::~G4Parameterisat
462
463 //--------------------------------------------
464 G4double G4ParameterisationPolyhedraZ::GetR(G4
465 G4d
466 G4d
467 {
468 // Linear parameterisation:
469 // r = az + b
470 // a = (r1 - r2)/(z1-z2)
471 // b = r1 - a*z1
472
473 return (r1-r2)/(z1-z2)*z + ( r1 - (r1-r2)/(z
474 }
475
476 //--------------------------------------------
477 G4double G4ParameterisationPolyhedraZ::GetRmin
478 {
479 // Get Rmin in the given z position for the gi
480
481 return GetR(z,
482 fOrigParamMother->Z_values[nseg]
483 fOrigParamMother->Rmin[nseg],
484 fOrigParamMother->Z_values[nseg+
485 fOrigParamMother->Rmin[nseg+1]);
486 }
487
488 //--------------------------------------------
489 G4double G4ParameterisationPolyhedraZ::GetRmax
490 {
491 // Get Rmax in the given z position for the gi
492
493 return GetR(z,
494 fOrigParamMother->Z_values[nseg]
495 fOrigParamMother->Rmax[nseg],
496 fOrigParamMother->Z_values[nseg+
497 fOrigParamMother->Rmax[nseg+1]);
498 }
499
500 //--------------------------------------------
501 G4double G4ParameterisationPolyhedraZ::GetMaxP
502 {
503 return std::abs(fOrigParamMother->Z_values[f
504 -fOrigParamMother->Z_values[0
505 }
506
507 //--------------------------------------------
508 void G4ParameterisationPolyhedraZ::CheckParame
509 {
510 G4VDivisionParameterisation::CheckParameters
511
512 // Division will be following the mother pol
513 //
514 if( fDivisionType == DivNDIV )
515 {
516 if( fOrigParamMother->Num_z_planes-1 != fn
517 {
518 std::ostringstream message;
519 message << "Configuration not supported.
520 << "Division along Z will be don
521 << G4endl
522 << "Z planes, i.e, the number of
523 << fOrigParamMother->Num_z_plane
524 << fnDiv << " !";
525 G4Exception("G4ParameterisationPolyhedra
526 "GeomDiv0001", FatalExceptio
527 }
528 }
529
530 // Division will be done within one polyhedr
531 // with applying given width and offset
532 //
533 if( fDivisionType == DivNDIVandWIDTH || fDiv
534 {
535 // Check if divided region does not span o
536 // than one z segment
537
538 G4int isegstart = -1; // number of the se
539 G4int isegend = -1; // number of the se
540
541 if ( !fReflectedSolid )
542 {
543 // The start/end position of the divided
544 //
545 G4double zstart = fOrigParamMother->Z_va
546 G4double zend = fOrigParamMother->Z_valu
547 + foffset + fnDiv*fwidth;
548
549 G4int counter = 0;
550 while ( isegend < 0 && counter < fOrigPa
551 {
552 // first segment
553 if ( zstart >= fOrigParamMother->Z_val
554 zstart < fOrigParamMother->Z_val
555 {
556 isegstart = counter;
557 }
558 // last segment
559 if ( zend > fOrigParamMother->Z_value
560 zend <= fOrigParamMother->Z_value
561 {
562 isegend = counter;
563 }
564 ++counter;
565 } // Loop checking, 06.08.2015, G.Cosmo
566 }
567 else
568 {
569 // The start/end position of the divided
570 //
571 G4double zstart = fOrigParamMother->Z_va
572 G4double zend = fOrigParamMother->Z_valu
573 - (foffset + fnDiv* fwidth
574
575 G4int counter = 0;
576 while ( isegend < 0 && counter < fOrigPa
577 {
578 // first segment
579 if ( zstart <= fOrigParamMother->Z_val
580 zstart > fOrigParamMother->Z_val
581 {
582 isegstart = counter;
583 }
584 // last segment
585 if ( zend < fOrigParamMother->Z_value
586 zend >= fOrigParamMother->Z_value
587 {
588 isegend = counter;
589 }
590 ++counter;
591 } // Loop checking, 06.08.2015, G.Cosmo
592 }
593
594 if ( isegstart != isegend )
595 {
596 std::ostringstream message;
597 message << "Configuration not supported.
598 << "Division with user defined w
599 << "Solid " << fmotherSolid->Get
600 << "Divided region is not betwee
601 G4Exception("G4ParameterisationPolyhedra
602 "GeomDiv0001", FatalExceptio
603 }
604
605 fNSegment = isegstart;
606 }
607 }
608
609 //--------------------------------------------
610 void
611 G4ParameterisationPolyhedraZ::
612 ComputeTransformation( const G4int copyNo, G4V
613 {
614 G4double posi;
615 if ( fDivisionType == DivNDIV )
616 {
617 // The position of the centre of copyNo-th
618
619 posi = ( fOrigParamMother->Z_values[copyNo
620 + fOrigParamMother->Z_valu
621 physVol->SetTranslation( G4ThreeVector(0,
622 }
623
624 if ( fDivisionType == DivNDIVandWIDTH || fDi
625 {
626 // The position of the centre of copyNo-th
627
628 posi = fOrigParamMother->Z_values[0];
629
630 if ( !fReflectedSolid )
631 posi += foffset + (2*copyNo + 1) * fwidt
632 else
633 posi -= foffset + (2*copyNo + 1) * fwidt
634
635 physVol->SetTranslation( G4ThreeVector(0,
636 }
637
638 //----- calculate rotation matrix: unit
639
640 #ifdef G4DIVDEBUG
641 if( verbose >= 2 )
642 {
643 G4cout << " G4ParameterisationPolyhedraZ -
644 << " copyNo: " << copyNo << " - fof
645 << " - fwidth: " << fwidth/CLHEP::d
646 }
647 #endif
648
649 ChangeRotMatrix( physVol );
650
651 #ifdef G4DIVDEBUG
652 if( verbose >= 2 )
653 {
654 G4cout << std::setprecision(8) << " G4Para
655 << copyNo << G4endl
656 << " Position: (0,0,0) - Width: " <
657 << " - Axis: " << faxis << G4endl;
658 }
659 #endif
660 }
661
662 //--------------------------------------------
663 void
664 G4ParameterisationPolyhedraZ::
665 ComputeDimensions( G4Polyhedra& phedra, const
666 const G4VPhysicalVolume* )
667 {
668 // Define division solid
669 //
670 G4PolyhedraHistorical origparam;
671 G4int nz = 2;
672 origparam.Num_z_planes = nz;
673 origparam.numSide = fOrigParamMother->numSid
674 origparam.Start_angle = fOrigParamMother->St
675 origparam.Opening_angle = fOrigParamMother->
676
677 // Define division solid z sections
678 //
679 origparam.Z_values = new G4double[nz];
680 origparam.Rmin = new G4double[nz];
681 origparam.Rmax = new G4double[nz];
682 origparam.Z_values[0] = - fwidth/2.;
683 origparam.Z_values[1] = fwidth/2.;
684
685 if ( fDivisionType == DivNDIV )
686 {
687 // The position of the centre of copyNo-th
688 //
689 G4double posi = ( fOrigParamMother->Z_valu
690 + fOrigParamMother->Z_valu
691
692 origparam.Z_values[0] = fOrigParamMother->
693 origparam.Z_values[1] = fOrigParamMother->
694 origparam.Rmin[0] = fOrigParamMother->Rmin
695 origparam.Rmin[1] = fOrigParamMother->Rmin
696 origparam.Rmax[0] = fOrigParamMother->Rmax
697 origparam.Rmax[1] = fOrigParamMother->Rmax
698 }
699
700 if ( fDivisionType == DivNDIVandWIDTH || fDi
701 {
702 if ( !fReflectedSolid )
703 {
704 origparam.Z_values[0] = -fwidth/2.;
705 origparam.Z_values[1] = fwidth/2.;
706
707 // The position of the centre of copyNo-
708 //
709 G4double posi = fOrigParamMother->Z_valu
710 + foffset + (2*copyNo + 1)
711
712 // The first and last z sides z values
713 G4double zstart = posi - fwidth/2.;
714 G4double zend = posi + fwidth/2.;
715 origparam.Rmin[0] = GetRmin(zstart, fNSe
716 origparam.Rmax[0] = GetRmax(zstart, fNSe
717 origparam.Rmin[1] = GetRmin(zend, fNSegm
718 origparam.Rmax[1] = GetRmax(zend, fNSegm
719 }
720 else
721 {
722 origparam.Z_values[0] = fwidth/2.;
723 origparam.Z_values[1] = -fwidth/2.;
724
725 // The position of the centre of copyNo-
726 //
727 G4double posi = fOrigParamMother->Z_valu
728 - ( foffset + (2*copyNo +
729
730 // The first and last z sides z values
731 //
732 G4double zstart = posi + fwidth/2.;
733 G4double zend = posi - fwidth/2.;
734 origparam.Rmin[0] = GetRmin(zstart, fNSe
735 origparam.Rmax[0] = GetRmax(zstart, fNSe
736 origparam.Rmin[1] = GetRmin(zend, fNSegm
737 origparam.Rmax[1] = GetRmax(zend, fNSegm
738 }
739
740 // It can happen due to rounding errors
741 //
742 if ( origparam.Rmin[0] < 0.0 ) origpara
743 if ( origparam.Rmin[nz-1] < 0.0 ) origpara
744 }
745
746 phedra.SetOriginalParameters(&origparam); /
747 phedra.Reset(); /
748
749 #ifdef G4DIVDEBUG
750 if( verbose >= 2 )
751 {
752 G4cout << "G4ParameterisationPolyhedraZ::C
753 << "-- Parametrised phedra copy-num
754 phedra.DumpInfo();
755 }
756 #endif
757 }
758