| 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 // G4ParameterisationPolycone[Rho/Phi/Z] imple
27 //
28 // 26.05.03 - P.Arce, Initial version
29 // 08.04.04 - I.Hrivnacova, Implemented reflec
30 //--------------------------------------------
31
32 #include "G4ParameterisationPolycone.hh"
33
34 #include <iomanip>
35 #include "G4ThreeVector.hh"
36 #include "G4RotationMatrix.hh"
37 #include "G4VPhysicalVolume.hh"
38 #include "G4LogicalVolume.hh"
39 #include "G4ReflectedSolid.hh"
40
41 //--------------------------------------------
42 G4VParameterisationPolycone::
43 G4VParameterisationPolycone( EAxis axis, G4int
44 G4double offset,
45 DivisionType divT
46 : G4VDivisionParameterisation( axis, nDiv,
47 {
48 /*#ifdef G4MULTITHREADED
49 std::ostringstream message;
50 message << "Divisions for G4Polycone curren
51 << G4endl
52 << "Sorry! Solid: " << msolid->GetN
53 G4Exception("G4VParameterisationPolycone::G
54 "GeomDiv0001", FatalException,
55 #endif */
56 auto msol = (G4Polycone*)(msolid);
57 if (msolid->GetEntityType() == "G4ReflectedS
58 {
59 // Get constituent solid
60 //
61 G4VSolid* mConstituentSolid
62 = ((G4ReflectedSolid*)msolid)->GetConst
63 msol = (G4Polycone*)(mConstituentSolid);
64
65 // Get parameters
66 //
67 G4int nofZplanes = msol->GetOriginalPara
68 G4double* zValues = msol->GetOriginalPara
69 G4double* rminValues = msol->GetOriginalP
70 G4double* rmaxValues = msol->GetOriginalP
71
72 // Invert z values
73 //
74 auto zValuesRefl = new G4double[nofZplanes
75 for (G4int i=0; i<nofZplanes; ++i) { zVal
76
77 auto newSolid
78 = new G4Polycone(msol->GetName(),
79 msol->GetStartPhi(),
80 msol->GetEndPhi() - mso
81 nofZplanes, zValuesRefl
82
83 delete [] zValuesRefl;
84
85 msol = newSolid;
86 fmotherSolid = newSolid;
87 fReflectedSolid = true;
88 fDeleteSolid = true;
89 }
90 }
91
92 //--------------------------------------------
93 G4VParameterisationPolycone::~G4VParameterisat
94
95 //--------------------------------------------
96 G4ParameterisationPolyconeRho::
97 G4ParameterisationPolyconeRho( EAxis axis, G4i
98 G4double width,
99 G4VSolid* msoli
100 : G4VParameterisationPolycone( axis, nDiv,
101 {
102 CheckParametersValidity();
103 SetType( "DivisionPolyconeRho" );
104
105 auto msol = (G4Polycone*)(fmotherSolid);
106 G4PolyconeHistorical* origparamMother = msol
107
108 if( divType == DivWIDTH )
109 {
110 fnDiv = CalculateNDiv( origparamMother->Rm
111 - origparamMother->Rm
112 }
113 else if( divType == DivNDIV )
114 {
115 fwidth = CalculateWidth( origparamMother->
116 - origparamMother->
117 }
118
119 #ifdef G4DIVDEBUG
120 if( verbose >= -1 )
121 {
122 G4cout << " G4ParameterisationPolyconeRho
123 << " = " << nDiv << G4endl
124 << " Offset " << foffset << " = " <
125 << " Width " << fwidth << " = " <<
126 }
127 #endif
128 }
129
130 //--------------------------------------------
131 G4ParameterisationPolyconeRho::~G4Parameterisa
132
133 //--------------------------------------------
134 void G4ParameterisationPolyconeRho::CheckParam
135 {
136 G4VDivisionParameterisation::CheckParameters
137
138 auto msol = (G4Polycone*)(fmotherSolid);
139
140 if( fDivisionType == DivNDIVandWIDTH || fDiv
141 {
142 std::ostringstream message;
143 message << "In solid " << msol->GetName()
144 << "Division along R will be done
145 << "different for each solid secti
146 << "WIDTH will not be used !";
147 G4Exception("G4VParameterisationPolycone::
148 "GeomDiv1001", JustWarning, me
149 }
150 if( foffset != 0. )
151 {
152 std::ostringstream message;
153 message << "In solid " << msol->GetName()
154 << "Division along R will be done
155 << "different for each solid secti
156 << "OFFSET will not be used !";
157 G4Exception("G4VParameterisationPolycone::
158 "GeomDiv1001", JustWarning, me
159 }
160 }
161
162 //--------------------------------------------
163 G4double G4ParameterisationPolyconeRho::GetMax
164 {
165 auto msol = (G4Polycone*)(fmotherSolid);
166 G4PolyconeHistorical* original_pars = msol->
167 return original_pars->Rmax[0] - original_par
168 }
169
170
171 //--------------------------------------------
172 void
173 G4ParameterisationPolyconeRho::
174 ComputeTransformation( const G4int, G4VPhysica
175 {
176 //----- translation
177 G4ThreeVector origin(0.,0.,0.);
178 //----- set translation
179 physVol->SetTranslation( origin );
180
181 //----- calculate rotation matrix: unit
182
183 #ifdef G4DIVDEBUG
184 if( verbose >= 2 )
185 {
186 G4cout << " G4ParameterisationPolyconeRho
187 << " foffset: " << foffset
188 << " - fwidth: " << fwidth << G4end
189 }
190 #endif
191
192 ChangeRotMatrix( physVol );
193
194 #ifdef G4DIVDEBUG
195 if( verbose >= 2 )
196 {
197 G4cout << std::setprecision(8) << " G4Para
198 << G4endl
199 << " Position: (0,0,0)"
200 << " - Width: " << fwidth/CLHEP::de
201 << " - Axis: " << faxis << G4endl;
202 }
203 #endif
204 }
205
206 //--------------------------------------------
207 void
208 G4ParameterisationPolyconeRho::
209 ComputeDimensions( G4Polycone& pcone, const G4
210 const G4VPhysicalVolume* )
211 {
212 auto msol = (G4Polycone*)(fmotherSolid);
213
214 G4PolyconeHistorical* origparamMother = msol
215 G4PolyconeHistorical origparam( *origparamMo
216 G4int nZplanes = origparamMother->Num_z_plan
217
218 G4double width = 0.;
219 for( G4int ii = 0; ii < nZplanes; ++ii )
220 {
221 width = CalculateWidth( origparamMother->R
222 - origparamMother->R
223 origparam.Rmin[ii] = origparamMother->Rmin
224 origparam.Rmax[ii] = origparamMother->Rmin
225 }
226
227 pcone.SetOriginalParameters(&origparam); //
228 pcone.Reset(); //
229
230 #ifdef G4DIVDEBUG
231 if( verbose >= -2 )
232 {
233 G4cout << "G4ParameterisationPolyconeRho::
234 << "-- Parametrised pcone copy-numb
235 pcone.DumpInfo();
236 }
237 #endif
238 }
239
240 //--------------------------------------------
241 G4ParameterisationPolyconePhi::
242 G4ParameterisationPolyconePhi( EAxis axis, G4i
243 G4double width,
244 G4VSolid* msoli
245 : G4VParameterisationPolycone( axis, nDiv,
246 {
247 CheckParametersValidity();
248 SetType( "DivisionPolyconePhi" );
249
250 auto msol = (G4Polycone*)(fmotherSolid);
251 G4double deltaPhi = msol->GetEndPhi() - msol
252
253 if( divType == DivWIDTH )
254 {
255 fnDiv = CalculateNDiv( deltaPhi, width, of
256 }
257 else if( divType == DivNDIV )
258 {
259 fwidth = CalculateWidth( deltaPhi, nDiv, o
260 }
261
262 #ifdef G4DIVDEBUG
263 if( verbose >= 1 )
264 {
265 G4cout << " G4ParameterisationPolyconePhi
266 << " = " << nDiv << G4endl
267 << " Offset " << foffset/CLHEP::deg
268 << " Width " << fwidth/CLHEP::deg <
269 }
270 #endif
271 }
272
273 //--------------------------------------------
274 G4ParameterisationPolyconePhi::~G4Parameterisa
275
276 //--------------------------------------------
277 G4double G4ParameterisationPolyconePhi::GetMax
278 {
279 auto msol = (G4Polycone*)(fmotherSolid);
280 return msol->GetEndPhi() - msol->GetStartPhi
281 }
282
283 //--------------------------------------------
284 void
285 G4ParameterisationPolyconePhi::
286 ComputeTransformation( const G4int copyNo, G4V
287 {
288 //----- translation
289 G4ThreeVector origin(0.,0.,0.);
290 //----- set translation
291 physVol->SetTranslation( origin );
292
293 //----- calculate rotation matrix (so that a
294 G4double posi = foffset + copyNo*fwidth;
295
296 #ifdef G4DIVDEBUG
297 if( verbose >= 2 )
298 {
299 G4cout << " G4ParameterisationPolyconePhi
300 << G4endl
301 << " copyNo: " << copyNo << " - fof
302 << " - fwidth: " << fwidth/CLHEP::d
303 }
304 #endif
305
306 ChangeRotMatrix( physVol, -posi );
307
308 #ifdef G4DIVDEBUG
309 if( verbose >= 2 )
310 {
311 G4cout << std::setprecision(8) << " G4Para
312 << copyNo << G4endl
313 << " Position: (0,0,0) - Width: " << fwid
314 << " - Axis: " << faxis << G4endl;
315 }
316 #endif
317 }
318
319 //--------------------------------------------
320 void
321 G4ParameterisationPolyconePhi::
322 ComputeDimensions( G4Polycone& pcone, const G4
323 const G4VPhysicalVolume* )
324 {
325 auto msol = (G4Polycone*)(fmotherSolid);
326
327 G4PolyconeHistorical* origparamMother = msol
328 G4PolyconeHistorical origparam( *origparamMo
329 origparam.Start_angle = origparamMother->Sta
330 origparam.Opening_angle = fwidth;
331
332 pcone.SetOriginalParameters(&origparam); //
333 pcone.Reset(); //
334
335 #ifdef G4DIVDEBUG
336 if( verbose >= 2 )
337 {
338 G4cout << "G4ParameterisationPolyconePhi::
339 pcone.DumpInfo();
340 }
341 #endif
342 }
343
344 //--------------------------------------------
345 G4ParameterisationPolyconeZ::
346 G4ParameterisationPolyconeZ( EAxis axis, G4int
347 G4double width, G
348 G4VSolid* msolid,
349 : G4VParameterisationPolycone( axis, nDiv, w
350 fOrigParamMother(((G4Polycone*)fmotherSoli
351 {
352
353 CheckParametersValidity();
354 SetType( "DivisionPolyconeZ" );
355
356 if( divType == DivWIDTH )
357 {
358 fnDiv =
359 CalculateNDiv( fOrigParamMother->Z_value
360 - fOrigParamMother->Z_val
361 }
362 else if( divType == DivNDIV )
363 {
364 fwidth =
365 CalculateNDiv( fOrigParamMother->Z_value
366 - fOrigParamMother->Z_val
367 }
368
369 #ifdef G4DIVDEBUG
370 if( verbose >= 1 )
371 {
372 G4cout << " G4ParameterisationPolyconeZ -
373 << nDiv << G4endl
374 << " Offset " << foffset << " = " <
375 << " Width " << fwidth << " = " <<
376 }
377 #endif
378 }
379
380 //--------------------------------------------
381 G4ParameterisationPolyconeZ::~G4Parameterisati
382
383 //--------------------------------------------
384 G4double G4ParameterisationPolyconeZ::GetR(G4d
385 G4d
386 G4d
387 {
388 // Linear parameterisation:
389 // r = az + b
390 // a = (r1 - r2)/(z1-z2)
391 // b = r1 - a*z1
392
393 return (r1-r2)/(z1-z2)*z + ( r1 - (r1-r2)/(z
394 }
395
396 //--------------------------------------------
397 G4double G4ParameterisationPolyconeZ::GetRmin(
398 {
399 // Get Rmin in the given z position for the gi
400
401 return GetR(z,
402 fOrigParamMother->Z_values[nseg]
403 fOrigParamMother->Rmin[nseg],
404 fOrigParamMother->Z_values[nseg+
405 fOrigParamMother->Rmin[nseg+1]);
406 }
407
408 //--------------------------------------------
409 G4double G4ParameterisationPolyconeZ::GetRmax(
410 {
411 // Get Rmax in the given z position for the gi
412
413 return GetR(z,
414 fOrigParamMother->Z_values[nseg]
415 fOrigParamMother->Rmax[nseg],
416 fOrigParamMother->Z_values[nseg+
417 fOrigParamMother->Rmax[nseg+1]);
418 }
419
420 //--------------------------------------------
421 G4double G4ParameterisationPolyconeZ::GetMaxPa
422 {
423 return std::abs (fOrigParamMother->Z_values[
424 -fOrigParamMother->Z_values[
425 }
426
427 //--------------------------------------------
428 void G4ParameterisationPolyconeZ::CheckParamet
429 {
430 G4VDivisionParameterisation::CheckParameters
431
432 // Division will be following the mother pol
433 //
434 if( fDivisionType == DivNDIV )
435 {
436 if( fnDiv > fOrigParamMother->Num_z_planes
437 {
438 std::ostringstream error;
439 error << "Configuration not supported."
440 << "Division along Z will be done
441 << G4endl
442 << "Z planes, i.e, the number of
443 << fOrigParamMother->Num_z_planes
444 << ", instead of: " << fnDiv << "
445 G4Exception("G4ParameterisationPolyconeZ
446 "GeomDiv0001", FatalExceptio
447 }
448 }
449
450 // Division will be done within one polycone
451 // with applying given width and offset
452 //
453 if( fDivisionType == DivNDIVandWIDTH || fDiv
454 {
455 // Check if divided region does not span o
456 // than one z segment
457
458 G4int isegstart = -1; // number of the se
459 G4int isegend = -1; // number of the se
460
461 if ( !fReflectedSolid )
462 {
463 // The start/end position of the divided
464 //
465 G4double zstart
466 = fOrigParamMother->Z_values[0] + foff
467 G4double zend
468 = fOrigParamMother->Z_values[0] + foff
469
470 G4int counter = 0;
471 while ( isegend < 0 && counter < fOrigPa
472 {
473 // first segment
474 if ( zstart >= fOrigParamMother->Z_val
475 zstart < fOrigParamMother->Z_val
476 {
477 isegstart = counter;
478 }
479 // last segment
480 if ( zend > fOrigParamMother->Z_value
481 zend <= fOrigParamMother->Z_value
482 {
483 isegend = counter;
484 }
485 ++counter;
486 } // Loop checking, 06.08.2015, G.Cosmo
487 }
488 else
489 {
490 // The start/end position of the divided
491 //
492 G4double zstart
493 = fOrigParamMother->Z_values[0] - foff
494 G4double zend
495 = fOrigParamMother->Z_values[0] - ( fo
496
497 G4int counter = 0;
498 while ( isegend < 0 && counter < fOrigPa
499 {
500 // first segment
501 if ( zstart <= fOrigParamMother->Z_val
502 zstart > fOrigParamMother->Z_val
503 {
504 isegstart = counter;
505 }
506 // last segment
507 if ( zend < fOrigParamMother->Z_value
508 zend >= fOrigParamMother->Z_value
509 {
510 isegend = counter;
511 }
512 ++counter;
513 } // Loop checking, 06.08.2015, G.Cosmo
514 }
515
516
517 if ( isegstart != isegend )
518 {
519 std::ostringstream message;
520 message << "Condiguration not supported.
521 << "Division with user defined w
522 << "Solid " << fmotherSolid->Get
523 << "Divided region is not betwee
524 G4Exception("G4ParameterisationPolyconeZ
525 "GeomDiv0001", FatalExceptio
526 }
527
528 fNSegment = isegstart;
529 }
530 }
531
532 //--------------------------------------------
533 void
534 G4ParameterisationPolyconeZ::
535 ComputeTransformation( const G4int copyNo, G4V
536 {
537 G4double posi = 0.;
538 if ( fDivisionType == DivNDIV )
539 {
540 // The position of the centre of copyNo-th
541 //
542 posi = ( fOrigParamMother->Z_values[copyNo
543 + fOrigParamMother->Z_valu
544 physVol->SetTranslation( G4ThreeVector(0,
545 }
546
547 if ( fDivisionType == DivNDIVandWIDTH || fDi
548 {
549 // The position of the centre of copyNo-th
550 //
551 posi = fOrigParamMother->Z_values[0];
552
553 if ( !fReflectedSolid )
554 posi += foffset + (2*copyNo + 1) * fwidt
555 else
556 posi -= foffset + (2*copyNo + 1) * fwidt
557
558 physVol->SetTranslation( G4ThreeVector(0,
559 }
560
561 //----- calculate rotation matrix: unit
562
563 #ifdef G4DIVDEBUG
564 if( verbose >= 2 )
565 {
566 G4cout << " G4ParameterisationPolyconeZ -
567 << " copyNo: " << copyNo << " - fof
568 << " - fwidth: " << fwidth/CLHEP::d
569 }
570 #endif
571
572 ChangeRotMatrix( physVol );
573
574 #ifdef G4DIVDEBUG
575 if( verbose >= 2 )
576 {
577 G4cout << std::setprecision(8) << " G4Para
578 << copyNo << G4endl
579 << " Position: (0,0,0) - Width: " <
580 << " - Axis: " << faxis << G4endl;
581 }
582 #endif
583 }
584
585 //--------------------------------------------
586 void
587 G4ParameterisationPolyconeZ::
588 ComputeDimensions( G4Polycone& pcone, const G4
589 const G4VPhysicalVolume* )
590 {
591
592 // Define division solid
593 //
594 G4PolyconeHistorical origparam;
595 G4int nz = 2;
596 origparam.Num_z_planes = nz;
597 origparam.Start_angle = fOrigParamMother->St
598 origparam.Opening_angle = fOrigParamMother->
599
600 // Define division solid z sections
601 //
602 origparam.Z_values = new G4double[nz];
603 origparam.Rmin = new G4double[nz];
604 origparam.Rmax = new G4double[nz];
605
606 if ( fDivisionType == DivNDIV )
607 {
608 // The position of the centre of copyNo-th
609 G4double posi = (fOrigParamMother->Z_value
610 + fOrigParamMother->Z_value
611
612 origparam.Z_values[0] = fOrigParamMother->
613 origparam.Z_values[1] = fOrigParamMother->
614 origparam.Rmin[0] = fOrigParamMother->Rmin
615 origparam.Rmin[1] = fOrigParamMother->Rmin
616 origparam.Rmax[0] = fOrigParamMother->Rmax
617 origparam.Rmax[1] = fOrigParamMother->Rmax
618 }
619
620 if ( fDivisionType == DivNDIVandWIDTH || fDi
621 {
622 if ( !fReflectedSolid )
623 {
624 origparam.Z_values[0] = - fwidth/2.;
625 origparam.Z_values[1] = fwidth/2.;
626
627 // The position of the centre of copyNo-
628 //
629 G4double posi = fOrigParamMother->Z_valu
630 + foffset + (2*copyNo + 1)
631
632 // The first and last z sides z values
633 //
634 G4double zstart = posi - fwidth/2.;
635 G4double zend = posi + fwidth/2.;
636 origparam.Rmin[0] = GetRmin(zstart, fNSe
637 origparam.Rmax[0] = GetRmax(zstart, fNSe
638 origparam.Rmin[1] = GetRmin(zend, fNSegm
639 origparam.Rmax[1] = GetRmax(zend, fNSegm
640 }
641 else
642 {
643 origparam.Z_values[0] = fwidth/2.;
644 origparam.Z_values[1] = - fwidth/2.;
645
646 // The position of the centre of copyNo-
647 //
648 G4double posi = fOrigParamMother->Z_valu
649 - ( foffset + (2*copyNo +
650
651 // The first and last z sides z values
652 //
653 G4double zstart = posi + fwidth/2.;
654 G4double zend = posi - fwidth/2.;
655 origparam.Rmin[0] = GetRmin(zstart, fNSe
656 origparam.Rmax[0] = GetRmax(zstart, fNSe
657 origparam.Rmin[1] = GetRmin(zend, fNSegm
658 origparam.Rmax[1] = GetRmax(zend, fNSegm
659 }
660
661 // It can happen due to rounding errors
662 //
663 if ( origparam.Rmin[0] < 0.0 ) origpara
664 if ( origparam.Rmin[nz-1] < 0.0 ) origpara
665 }
666
667 pcone.SetOriginalParameters(&origparam); //
668 pcone.Reset(); //
669
670 #ifdef G4DIVDEBUG
671 if( verbose >= 2 )
672 {
673 G4cout << "G4ParameterisationPolyconeZ::Co
674 << "-- Parametrised pcone copy-numb
675 pcone.DumpInfo();
676 }
677 #endif
678 }
679