Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParameterisationPolycone.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4ParameterisationPolycone.cc 69784 2013-05-15 09:16:06Z gcosmo $
28 //
29 // class G4ParameterisationPolycone Implementation file
30 //
31 // 26.05.03 - P.Arce, Initial version
32 // 08.04.04 - I.Hrivnacova, Implemented reflection
33 //---------------------------------------------------------------------
34 
36 
37 #include <iomanip>
38 #include "G4ThreeVector.hh"
39 #include "G4RotationMatrix.hh"
40 #include "G4VPhysicalVolume.hh"
41 #include "G4LogicalVolume.hh"
42 #include "G4ReflectedSolid.hh"
43 
44 //-----------------------------------------------------------------------
47  G4double offset, G4VSolid* msolid,
48  DivisionType divType )
49  : G4VDivisionParameterisation( axis, nDiv, width, offset, divType, msolid )
50 {
51  G4Polycone* msol = (G4Polycone*)(msolid);
52  if ((msolid->GetEntityType() != "G4ReflectedSolid") && (msol->IsGeneric()))
53  {
54  std::ostringstream message;
55  message << "Generic construct for G4Polycone NOT supported." << G4endl
56  << "Sorry! Solid: " << msol->GetName() << ".";
57  G4Exception("G4VParameterisationPolycone::G4VParameterisationPolycone()",
58  "GeomDiv0001", FatalException, message);
59  }
60  if (msolid->GetEntityType() == "G4ReflectedSolid")
61  {
62  // Get constituent solid
63  G4VSolid* mConstituentSolid
64  = ((G4ReflectedSolid*)msolid)->GetConstituentMovedSolid();
65  msol = (G4Polycone*)(mConstituentSolid);
66 
67  // Get parameters
68  G4int nofZplanes = msol->GetOriginalParameters()->Num_z_planes;
69  G4double* zValues = msol->GetOriginalParameters()->Z_values;
70  G4double* rminValues = msol->GetOriginalParameters()->Rmin;
71  G4double* rmaxValues = msol->GetOriginalParameters()->Rmax;
72 
73  // Invert z values
74  G4double* zValuesRefl = new double[nofZplanes];
75  for (G4int i=0; i<nofZplanes; i++) zValuesRefl[i] = - zValues[i];
76 
77  G4Polycone* newSolid
78  = new G4Polycone(msol->GetName(),
79  msol->GetStartPhi(),
80  msol->GetEndPhi() - msol->GetStartPhi(),
81  nofZplanes, zValuesRefl, rminValues, rmaxValues);
82 
83  delete [] zValuesRefl;
84 
85  msol = newSolid;
86  fmotherSolid = newSolid;
87  fReflectedSolid = true;
88  fDeleteSolid = true;
89  }
90 }
91 
92 //---------------------------------------------------------------------
94 {
95 }
96 
97 //---------------------------------------------------------------------
100  G4double width, G4double offset,
101  G4VSolid* msolid, DivisionType divType )
102  : G4VParameterisationPolycone( axis, nDiv, width, offset, msolid, divType )
103 {
105  SetType( "DivisionPolyconeRho" );
106 
107  G4Polycone* msol = (G4Polycone*)(fmotherSolid);
108  G4PolyconeHistorical* origparamMother = msol->GetOriginalParameters();
109 
110  if( divType == DivWIDTH )
111  {
112  fnDiv = CalculateNDiv( origparamMother->Rmax[0]
113  - origparamMother->Rmin[0], width, offset );
114  }
115  else if( divType == DivNDIV )
116  {
117  fwidth = CalculateWidth( origparamMother->Rmax[0]
118  - origparamMother->Rmin[0], nDiv, offset );
119  }
120 
121 #ifdef G4DIVDEBUG
122  if( verbose >= -1 )
123  {
124  G4cout << " G4ParameterisationPolyconeRho - # divisions " << fnDiv
125  << " = " << nDiv << G4endl
126  << " Offset " << foffset << " = " << offset << G4endl
127  << " Width " << fwidth << " = " << width << G4endl;
128  }
129 #endif
130 }
131 
132 //---------------------------------------------------------------------
134 {
135 }
136 
137 //---------------------------------------------------------------------
139 {
141 
142  G4Polycone* msol = (G4Polycone*)(fmotherSolid);
143 
145  {
146  std::ostringstream message;
147  message << "In solid " << msol->GetName() << G4endl
148  << "Division along R will be done with a width "
149  << "different for each solid section." << G4endl
150  << "WIDTH will not be used !";
151  G4Exception("G4VParameterisationPolycone::CheckParametersValidity()",
152  "GeomDiv1001", JustWarning, message);
153  }
154  if( foffset != 0. )
155  {
156  std::ostringstream message;
157  message << "In solid " << msol->GetName() << G4endl
158  << "Division along R will be done with a width "
159  << "different for each solid section." << G4endl
160  << "OFFSET will not be used !";
161  G4Exception("G4VParameterisationPolycone::CheckParametersValidity()",
162  "GeomDiv1001", JustWarning, message);
163  }
164 }
165 
166 //------------------------------------------------------------------------
168 {
169  G4Polycone* msol = (G4Polycone*)(fmotherSolid);
170  G4PolyconeHistorical* original_pars = msol->GetOriginalParameters();
171  return original_pars->Rmax[0] - original_pars->Rmin[0];
172 }
173 
174 
175 //---------------------------------------------------------------------
176 void
179 {
180  //----- translation
181  G4ThreeVector origin(0.,0.,0.);
182  //----- set translation
183  physVol->SetTranslation( origin );
184 
185  //----- calculate rotation matrix: unit
186 
187 #ifdef G4DIVDEBUG
188  if( verbose >= 2 )
189  {
190  G4cout << " G4ParameterisationPolyconeRho " << G4endl
191  << " foffset: " << foffset
192  << " - fwidth: " << fwidth << G4endl;
193  }
194 #endif
195 
196  ChangeRotMatrix( physVol );
197 
198 #ifdef G4DIVDEBUG
199  if( verbose >= 2 )
200  {
201  G4cout << std::setprecision(8) << " G4ParameterisationPolyconeRho "
202  << G4endl
203  << " Position: " << origin/mm
204  << " - Width: " << fwidth/deg
205  << " - Axis: " << faxis << G4endl;
206  }
207 #endif
208 }
209 
210 //---------------------------------------------------------------------
211 void
213 ComputeDimensions( G4Polycone& pcone, const G4int copyNo,
214  const G4VPhysicalVolume* ) const
215 {
216  G4Polycone* msol = (G4Polycone*)(fmotherSolid);
217 
218  G4PolyconeHistorical* origparamMother = msol->GetOriginalParameters();
219  G4PolyconeHistorical origparam( *origparamMother );
220  G4int nZplanes = origparamMother->Num_z_planes;
221 
222  G4double width = 0.;
223  for( G4int ii = 0; ii < nZplanes; ii++ )
224  {
225  width = CalculateWidth( origparamMother->Rmax[ii]
226  - origparamMother->Rmin[ii], fnDiv, foffset );
227  origparam.Rmin[ii] = origparamMother->Rmin[ii]+foffset+width*copyNo;
228  origparam.Rmax[ii] = origparamMother->Rmin[ii]+foffset+width*(copyNo+1);
229  }
230 
231  pcone.SetOriginalParameters(&origparam); // copy values & transfer pointers
232  pcone.Reset(); // reset to new solid parameters
233 
234 #ifdef G4DIVDEBUG
235  if( verbose >= -2 )
236  {
237  G4cout << "G4ParameterisationPolyconeRho::ComputeDimensions()" << G4endl
238  << "-- Parametrised pcone copy-number: " << copyNo << G4endl;
239  pcone.DumpInfo();
240  }
241 #endif
242 }
243 
244 //---------------------------------------------------------------------
247  G4double width, G4double offset,
248  G4VSolid* msolid, DivisionType divType )
249  : G4VParameterisationPolycone( axis, nDiv, width, offset, msolid, divType )
250 {
252  SetType( "DivisionPolyconePhi" );
253 
254  G4Polycone* msol = (G4Polycone*)(fmotherSolid);
255  G4double deltaPhi = msol->GetEndPhi() - msol->GetStartPhi();
256 
257  if( divType == DivWIDTH )
258  {
259  fnDiv = CalculateNDiv( deltaPhi, width, offset );
260  }
261  else if( divType == DivNDIV )
262  {
263  fwidth = CalculateWidth( deltaPhi, nDiv, offset );
264  }
265 
266 #ifdef G4DIVDEBUG
267  if( verbose >= 1 )
268  {
269  G4cout << " G4ParameterisationPolyconePhi - # divisions " << fnDiv
270  << " = " << nDiv << G4endl
271  << " Offset " << foffset/deg << " = " << offset/deg << G4endl
272  << " Width " << fwidth/deg << " = " << width/deg << G4endl;
273  }
274 #endif
275 }
276 
277 //---------------------------------------------------------------------
279 {
280 }
281 
282 //------------------------------------------------------------------------
284 {
285  G4Polycone* msol = (G4Polycone*)(fmotherSolid);
286  return msol->GetEndPhi() - msol->GetStartPhi();
287 }
288 
289 //---------------------------------------------------------------------
290 void
292 ComputeTransformation( const G4int copyNo, G4VPhysicalVolume *physVol ) const
293 {
294  //----- translation
295  G4ThreeVector origin(0.,0.,0.);
296  //----- set translation
297  physVol->SetTranslation( origin );
298 
299  //----- calculate rotation matrix (so that all volumes point to the centre)
300  G4double posi = foffset + copyNo*fwidth;
301 
302 #ifdef G4DIVDEBUG
303  if( verbose >= 2 )
304  {
305  G4cout << " G4ParameterisationPolyconePhi - position: " << posi/deg
306  << G4endl
307  << " copyNo: " << copyNo << " - foffset: " << foffset/deg
308  << " - fwidth: " << fwidth/deg << G4endl;
309  }
310 #endif
311 
312  ChangeRotMatrix( physVol, -posi );
313 
314 #ifdef G4DIVDEBUG
315  if( verbose >= 2 )
316  {
317  G4cout << std::setprecision(8) << " G4ParameterisationPolyconePhi "
318  << copyNo << G4endl
319  << " Position: " << origin << " - Width: " << fwidth
320  << " - Axis: " << faxis << G4endl;
321  }
322 #endif
323 }
324 
325 //---------------------------------------------------------------------
326 void
329  const G4VPhysicalVolume* ) const
330 {
331  G4Polycone* msol = (G4Polycone*)(fmotherSolid);
332 
333  G4PolyconeHistorical* origparamMother = msol->GetOriginalParameters();
334  G4PolyconeHistorical origparam( *origparamMother );
335  origparam.Start_angle = origparamMother->Start_angle;
336  origparam.Opening_angle = fwidth;
337 
338  pcone.SetOriginalParameters(&origparam); // copy values & transfer pointers
339  pcone.Reset(); // reset to new solid parameters
340 
341 #ifdef G4DIVDEBUG
342  if( verbose >= 2 )
343  {
344  G4cout << "G4ParameterisationPolyconePhi::ComputeDimensions():" << G4endl;
345  pcone.DumpInfo();
346  }
347 #endif
348 }
349 
350 //---------------------------------------------------------------------
353  G4double width, G4double offset,
354  G4VSolid* msolid, DivisionType divType)
355  : G4VParameterisationPolycone( axis, nDiv, width, offset, msolid, divType ),
356  fNSegment(0),
357  fOrigParamMother(((G4Polycone*)fmotherSolid)->GetOriginalParameters())
358 {
359 
361  SetType( "DivisionPolyconeZ" );
362 
363  if( divType == DivWIDTH )
364  {
365  fnDiv =
366  CalculateNDiv( fOrigParamMother->Z_values[fOrigParamMother->Num_z_planes-1]
367  - fOrigParamMother->Z_values[0] , width, offset );
368  }
369  else if( divType == DivNDIV )
370  {
371  fwidth =
372  CalculateNDiv( fOrigParamMother->Z_values[fOrigParamMother->Num_z_planes-1]
373  - fOrigParamMother->Z_values[0] , nDiv, offset );
374  }
375 
376 #ifdef G4DIVDEBUG
377  if( verbose >= 1 )
378  {
379  G4cout << " G4ParameterisationPolyconeZ - # divisions " << fnDiv << " = "
380  << nDiv << G4endl
381  << " Offset " << foffset << " = " << offset << G4endl
382  << " Width " << fwidth << " = " << width << G4endl;
383  }
384 #endif
385 }
386 
387 //---------------------------------------------------------------------
389 {
390 }
391 
392 //------------------------------------------------------------------------
393 G4double G4ParameterisationPolyconeZ::GetR(G4double z,
394  G4double z1, G4double r1,
395  G4double z2, G4double r2) const
396 {
397  // Linear parameterisation:
398  // r = az + b
399  // a = (r1 - r2)/(z1-z2)
400  // b = r1 - a*z1
401 
402  return (r1-r2)/(z1-z2)*z + ( r1 - (r1-r2)/(z1-z2)*z1 ) ;
403 }
404 
405 //------------------------------------------------------------------------
406 G4double G4ParameterisationPolyconeZ::GetRmin(G4double z, G4int nseg) const
407 {
408 // Get Rmin in the given z position for the given polycone segment
409 
410  return GetR(z,
411  fOrigParamMother->Z_values[nseg],
412  fOrigParamMother->Rmin[nseg],
413  fOrigParamMother->Z_values[nseg+1],
414  fOrigParamMother->Rmin[nseg+1]);
415 }
416 
417 //------------------------------------------------------------------------
418 G4double G4ParameterisationPolyconeZ::GetRmax(G4double z, G4int nseg) const
419 {
420 // Get Rmax in the given z position for the given polycone segment
421 
422  return GetR(z,
423  fOrigParamMother->Z_values[nseg],
424  fOrigParamMother->Rmax[nseg],
425  fOrigParamMother->Z_values[nseg+1],
426  fOrigParamMother->Rmax[nseg+1]);
427 }
428 
429 //------------------------------------------------------------------------
431 {
432  return std::abs (fOrigParamMother->Z_values[fOrigParamMother->Num_z_planes-1]
433  -fOrigParamMother->Z_values[0]);
434 }
435 
436 //---------------------------------------------------------------------
438 {
440 
441  // Division will be following the mother polycone segments
442  if( fDivisionType == DivNDIV ) {
443  if( fnDiv > fOrigParamMother->Num_z_planes-1 ) {
444  std::ostringstream error;
445  error << "Configuration not supported." << G4endl
446  << "Division along Z will be done by splitting in the defined"
447  << G4endl
448  << "Z planes, i.e, the number of division would be: "
449  << fOrigParamMother->Num_z_planes-1
450  << ", instead of: " << fnDiv << " !";
451  G4Exception("G4ParameterisationPolyconeZ::CheckParametersValidity()",
452  "GeomDiv0001", FatalException, error);
453  }
454  }
455 
456  // Division will be done within one polycone segment
457  // with applying given width and offset
459  // Check if divided region does not span over more
460  // than one z segment
461 
462  G4int isegstart = -1; // number of the segment containing start position
463  G4int isegend = -1; // number of the segment containing end position
464 
465  if ( ! fReflectedSolid ) {
466  // The start/end position of the divided region
467  G4double zstart
468  = fOrigParamMother->Z_values[0] + foffset;
469  G4double zend
470  = fOrigParamMother->Z_values[0] + foffset + fnDiv* fwidth;
471 
472  G4int counter = 0;
473  while ( isegend < 0 && counter < fOrigParamMother->Num_z_planes-1 ) {
474  // first segment
475  if ( zstart >= fOrigParamMother->Z_values[counter] &&
476  zstart < fOrigParamMother->Z_values[counter+1] ) {
477  isegstart = counter;
478  }
479  // last segment
480  if ( zend > fOrigParamMother->Z_values[counter] &&
481  zend <= fOrigParamMother->Z_values[counter+1] ) {
482  isegend = counter;
483  }
484  ++counter;
485  }
486  }
487  else {
488  // The start/end position of the divided region
489  G4double zstart
490  = fOrigParamMother->Z_values[0] - foffset;
491  G4double zend
492  = fOrigParamMother->Z_values[0] - ( foffset + fnDiv* fwidth);
493 
494  G4int counter = 0;
495  while ( isegend < 0 && counter < fOrigParamMother->Num_z_planes-1 ) {
496  // first segment
497  if ( zstart <= fOrigParamMother->Z_values[counter] &&
498  zstart > fOrigParamMother->Z_values[counter+1] ) {
499  isegstart = counter;
500  }
501  // last segment
502  if ( zend < fOrigParamMother->Z_values[counter] &&
503  zend >= fOrigParamMother->Z_values[counter+1] ) {
504  isegend = counter;
505  }
506  ++counter;
507  }
508  }
509 
510 
511  if ( isegstart != isegend ) {
512  std::ostringstream message;
513  message << "Condiguration not supported." << G4endl
514  << "Division with user defined width." << G4endl
515  << "Solid " << fmotherSolid->GetName() << G4endl
516  << "Divided region is not between two z planes.";
517  G4Exception("G4ParameterisationPolyconeZ::CheckParametersValidity()",
518  "GeomDiv0001", FatalException, message);
519  }
520 
521  fNSegment = isegstart;
522  }
523 }
524 
525 //---------------------------------------------------------------------
526 void
528 ComputeTransformation( const G4int copyNo, G4VPhysicalVolume* physVol) const
529 {
530  if ( fDivisionType == DivNDIV ) {
531  // The position of the centre of copyNo-th mother polycone segment
532  G4double posi
533  = ( fOrigParamMother->Z_values[copyNo]
534  + fOrigParamMother->Z_values[copyNo+1])/2;
535  physVol->SetTranslation( G4ThreeVector(0, 0, posi) );
536  }
537 
539  // The position of the centre of copyNo-th division
540  G4double posi = fOrigParamMother->Z_values[0];
541 
542  if ( ! fReflectedSolid )
543  posi += foffset + (2*copyNo + 1) * fwidth/2.;
544  else
545  posi -= foffset + (2*copyNo + 1) * fwidth/2.;
546 
547  physVol->SetTranslation( G4ThreeVector(0, 0, posi) );
548  }
549 
550  //----- calculate rotation matrix: unit
551 
552 #ifdef G4DIVDEBUG
553  if( verbose >= 2 )
554  {
555  G4cout << " G4ParameterisationPolyconeZ - position: " << posi << G4endl
556  << " copyNo: " << copyNo << " - foffset: " << foffset/deg
557  << " - fwidth: " << fwidth/deg << G4endl;
558  }
559 #endif
560 
561  ChangeRotMatrix( physVol );
562 
563 #ifdef G4DIVDEBUG
564  if( verbose >= 2 )
565  {
566  G4cout << std::setprecision(8) << " G4ParameterisationPolyconeZ "
567  << copyNo << G4endl
568  << " Position: " << origin << " - Width: " << fwidth
569  << " - Axis: " << faxis << G4endl;
570  }
571 #endif
572 }
573 
574 //---------------------------------------------------------------------
575 void
577 ComputeDimensions( G4Polycone& pcone, const G4int copyNo,
578  const G4VPhysicalVolume* ) const
579 {
580 
581  // Define division solid
582  G4PolyconeHistorical origparam;
583  G4int nz = 2;
584  origparam.Num_z_planes = nz;
585  origparam.Start_angle = fOrigParamMother->Start_angle;
586  origparam.Opening_angle = fOrigParamMother->Opening_angle;
587 
588  // Define division solid z sections
589  origparam.Z_values = new G4double[nz];
590  origparam.Rmin = new G4double[nz];
591  origparam.Rmax = new G4double[nz];
592 
593  if ( fDivisionType == DivNDIV ) {
594  // The position of the centre of copyNo-th mother polycone segment
595  G4double posi = (fOrigParamMother->Z_values[copyNo]
596  + fOrigParamMother->Z_values[copyNo+1])/2;
597 
598  origparam.Z_values[0] = fOrigParamMother->Z_values[copyNo] - posi;
599  origparam.Z_values[1] = fOrigParamMother->Z_values[copyNo+1] - posi;
600  origparam.Rmin[0] = fOrigParamMother->Rmin[copyNo];
601  origparam.Rmin[1] = fOrigParamMother->Rmin[copyNo+1];
602  origparam.Rmax[0] = fOrigParamMother->Rmax[copyNo];
603  origparam.Rmax[1] = fOrigParamMother->Rmax[copyNo+1];
604  }
605 
607  if ( ! fReflectedSolid ) {
608  origparam.Z_values[0] = - fwidth/2.;
609  origparam.Z_values[1] = fwidth/2.;
610 
611  // The position of the centre of copyNo-th division
612  G4double posi
613  = fOrigParamMother->Z_values[0] + foffset + (2*copyNo + 1) * fwidth/2.;
614 
615  // The first and last z sides z values
616  G4double zstart = posi - fwidth/2.;
617  G4double zend = posi + fwidth/2.;
618  origparam.Rmin[0] = GetRmin(zstart, fNSegment);
619  origparam.Rmax[0] = GetRmax(zstart, fNSegment);
620  origparam.Rmin[1] = GetRmin(zend, fNSegment);
621  origparam.Rmax[1] = GetRmax(zend, fNSegment);
622  }
623  else {
624  origparam.Z_values[0] = fwidth/2.;
625  origparam.Z_values[1] = - fwidth/2.;
626 
627  // The position of the centre of copyNo-th division
628  G4double posi
629  = fOrigParamMother->Z_values[0] - ( foffset + (2*copyNo + 1) * fwidth/2.);
630 
631  // The first and last z sides z values
632  G4double zstart = posi + fwidth/2.;
633  G4double zend = posi - fwidth/2.;
634  origparam.Rmin[0] = GetRmin(zstart, fNSegment);
635  origparam.Rmax[0] = GetRmax(zstart, fNSegment);
636  origparam.Rmin[1] = GetRmin(zend, fNSegment);
637  origparam.Rmax[1] = GetRmax(zend, fNSegment);
638  }
639 
640  // It can happen due to rounding errors
641  if ( origparam.Rmin[0] < 0.0 ) origparam.Rmin[0] = 0.0;
642  if ( origparam.Rmin[nz-1] < 0.0 ) origparam.Rmin[1] = 0.0;
643  }
644 
645  pcone.SetOriginalParameters(&origparam); // copy values & transfer pointers
646  pcone.Reset(); // reset to new solid parameters
647 
648 #ifdef G4DIVDEBUG
649  if( verbose >= 2 )
650  {
651  G4cout << "G4ParameterisationPolyconeZ::ComputeDimensions()" << G4endl
652  << "-- Parametrised pcone copy-number: " << copyNo << G4endl;
653  pcone.DumpInfo();
654  }
655 #endif
656 }