Geant4 Cross Reference |
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 // G4ParameterisationCons[Rho/Phi/Z] implementation 27 // 28 // 26.05.03 - P.Arce, Initial version 29 // 08.04.04 - I.Hrivnacova, Implemented reflection 30 // 21.04.10 - M.Asai, Added gaps 31 // -------------------------------------------------------------------- 32 33 #include "G4ParameterisationCons.hh" 34 35 #include <iomanip> 36 #include "G4ThreeVector.hh" 37 #include "G4RotationMatrix.hh" 38 #include "G4VPhysicalVolume.hh" 39 #include "G4LogicalVolume.hh" 40 #include "G4ReflectedSolid.hh" 41 #include "G4Cons.hh" 42 43 //-------------------------------------------------------------------------- 44 G4VParameterisationCons:: 45 G4VParameterisationCons( EAxis axis, G4int nDiv, G4double width, 46 G4double offset, G4VSolid* msolid, 47 DivisionType divType ) 48 : G4VDivisionParameterisation( axis, nDiv, width, offset, divType, msolid ) 49 { 50 auto msol = (G4Cons*)(msolid); 51 if (msolid->GetEntityType() == "G4ReflectedSolid") 52 { 53 // Get constituent solid 54 G4VSolid* mConstituentSolid 55 = ((G4ReflectedSolid*)msolid)->GetConstituentMovedSolid(); 56 msol = (G4Cons*)(mConstituentSolid); 57 58 // Create a new solid with inversed parameters 59 auto newSolid 60 = new G4Cons(msol->GetName(), 61 msol->GetInnerRadiusPlusZ(), msol->GetOuterRadiusPlusZ(), 62 msol->GetInnerRadiusMinusZ(), msol->GetOuterRadiusMinusZ(), 63 msol->GetZHalfLength(), 64 msol->GetStartPhiAngle(), msol->GetDeltaPhiAngle()); 65 msol = newSolid; 66 fmotherSolid = newSolid; 67 fReflectedSolid = true; 68 fDeleteSolid = true; 69 } 70 } 71 72 //------------------------------------------------------------------------ 73 G4VParameterisationCons::~G4VParameterisationCons() = default; 74 75 //-------------------------------------------------------------------------- 76 G4ParameterisationConsRho:: 77 G4ParameterisationConsRho( EAxis axis, G4int nDiv, 78 G4double width, G4double offset, 79 G4VSolid* msolid, DivisionType divType ) 80 : G4VParameterisationCons( axis, nDiv, width, offset, msolid, divType ) 81 { 82 CheckParametersValidity(); 83 SetType( "DivisionConsRho" ); 84 85 auto msol = (G4Cons*)(fmotherSolid); 86 if( msol->GetInnerRadiusPlusZ() == 0. ) 87 { 88 std::ostringstream message; 89 message << "OuterRadiusMinusZ = 0" << G4endl 90 << "Width is calculated as that of OuterRadiusMinusZ !"; 91 G4Exception("G4ParameterisationConsRho::G4ParameterisationConsRho()", 92 "GeomDiv1001", JustWarning, message); 93 } 94 95 if( divType == DivWIDTH ) 96 { 97 fnDiv = CalculateNDiv( msol->GetOuterRadiusMinusZ() 98 - msol->GetInnerRadiusMinusZ(), width, offset ); 99 } 100 else if( divType == DivNDIV ) 101 { 102 auto mconsol = (G4Cons*)(msolid); 103 fwidth = CalculateWidth( mconsol->GetOuterRadiusMinusZ() 104 - mconsol->GetInnerRadiusMinusZ(), nDiv, offset ); 105 } 106 107 #ifdef G4DIVDEBUG 108 if( verbose >= 1 ) 109 { 110 G4cout << " G4ParameterisationConsRho - no divisions " << fnDiv << " = " 111 << nDiv << G4endl 112 << " Offset " << foffset << " = " << offset 113 << " - Width " << fwidth << " = " << width << G4endl; 114 } 115 #endif 116 } 117 118 //-------------------------------------------------------------------------- 119 G4ParameterisationConsRho::~G4ParameterisationConsRho() = default; 120 121 //------------------------------------------------------------------------ 122 G4double G4ParameterisationConsRho::GetMaxParameter() const 123 { 124 auto msol = (G4Cons*)(fmotherSolid); 125 return msol->GetOuterRadiusMinusZ() - msol->GetInnerRadiusMinusZ(); 126 } 127 128 //-------------------------------------------------------------------------- 129 void 130 G4ParameterisationConsRho:: 131 ComputeTransformation( const G4int, G4VPhysicalVolume *physVol ) const 132 { 133 //----- translation 134 G4ThreeVector origin(0.,0.,0.); 135 //----- set translation 136 physVol->SetTranslation( origin ); 137 138 //----- calculate rotation matrix: unit 139 140 #ifdef G4DIVDEBUG 141 if( verbose >= 2 ) 142 { 143 G4cout << " G4ParameterisationConsRho " << G4endl 144 << " Offset: " << foffset 145 << " - Width: " << fwidth << G4endl; 146 } 147 #endif 148 149 ChangeRotMatrix( physVol ); 150 151 #ifdef G4DIVDEBUG 152 if( verbose >= 2 ) 153 { 154 G4cout << std::setprecision(8) << " G4ParameterisationConsRho" << G4endl 155 << " Position: " << origin << " - Width: " << fwidth 156 << " - Axis: " << faxis << G4endl; 157 } 158 #endif 159 } 160 161 //-------------------------------------------------------------------------- 162 void 163 G4ParameterisationConsRho:: 164 ComputeDimensions( G4Cons& cons, const G4int copyNo, 165 const G4VPhysicalVolume* ) const 166 { 167 auto msol = (G4Cons*)(fmotherSolid); 168 169 G4double pRMin1 = msol->GetInnerRadiusMinusZ() + foffset + fwidth*copyNo; 170 G4double pRMax1 = msol->GetInnerRadiusMinusZ() + foffset + fwidth*(copyNo+1); 171 172 //width at Z Plus 173 //- G4double fwidthPlus = 174 // fwidth * ( msol->GetOuterRadiusPlusZ()/ msol->GetInnerRadiusPlusZ()) 175 //- / ( msol->GetOuterRadiusMinusZ() - msol->GetInnerRadiusMinusZ()); 176 G4double fwidthPlus = CalculateWidth( msol->GetOuterRadiusPlusZ() 177 - msol->GetInnerRadiusPlusZ(), fnDiv, foffset ); 178 G4double pRMin2 = msol->GetInnerRadiusPlusZ() 179 + foffset + fwidthPlus * copyNo; 180 G4double pRMax2 = msol->GetInnerRadiusPlusZ() 181 + foffset + fwidthPlus * (copyNo+1); 182 G4double pDz = msol->GetZHalfLength(); 183 184 G4double d_half_gap = fhgap * pRMax2 / pRMax1; 185 //- already rotated double pSR = foffset + copyNo*fwidth; 186 G4double pSPhi = msol->GetStartPhiAngle(); 187 G4double pDPhi = msol->GetDeltaPhiAngle();; 188 189 cons.SetInnerRadiusMinusZ( pRMin1 + fhgap ); 190 cons.SetOuterRadiusMinusZ( pRMax1 - fhgap ); 191 cons.SetInnerRadiusPlusZ( pRMin2 + d_half_gap ); 192 cons.SetOuterRadiusPlusZ( pRMax2 - d_half_gap ); 193 cons.SetZHalfLength( pDz ); 194 cons.SetStartPhiAngle( pSPhi, false ); 195 cons.SetDeltaPhiAngle( pDPhi ); 196 197 #ifdef G4DIVDEBUG 198 if( verbose >= 2 ) 199 { 200 G4cout << " G4ParameterisationConsRho::ComputeDimensions()" << G4endl 201 << " pRMin: " << pRMin1 << " - pRMax: " << pRMax1 << G4endl; 202 if( verbose >= 4 ) cons.DumpInfo(); 203 } 204 #endif 205 } 206 207 //-------------------------------------------------------------------------- 208 G4ParameterisationConsPhi:: 209 G4ParameterisationConsPhi( EAxis axis, G4int nDiv, 210 G4double width, G4double offset, 211 G4VSolid* msolid, DivisionType divType ) 212 : G4VParameterisationCons( axis, nDiv, width, offset, msolid, divType ) 213 { 214 CheckParametersValidity(); 215 SetType( "DivisionConsPhi" ); 216 217 auto msol = (G4Cons*)(fmotherSolid); 218 if( divType == DivWIDTH ) 219 { 220 fnDiv = CalculateNDiv( msol->GetDeltaPhiAngle(), width, offset ); 221 } 222 else if( divType == DivNDIV ) 223 { 224 fwidth = CalculateWidth( msol->GetDeltaPhiAngle(), nDiv, offset ); 225 } 226 227 #ifdef G4DIVDEBUG 228 if( verbose >= 1 ) 229 { 230 G4cout << " G4ParameterisationConsPhi no divisions " << fnDiv << " = " 231 << nDiv << G4endl 232 << " Offset " << foffset << " = " << offset << G4endl 233 << " Width " << fwidth << " = " << width << G4endl; 234 } 235 #endif 236 } 237 238 //-------------------------------------------------------------------------- 239 G4ParameterisationConsPhi::~G4ParameterisationConsPhi() = default; 240 241 //------------------------------------------------------------------------ 242 G4double G4ParameterisationConsPhi::GetMaxParameter() const 243 { 244 auto msol = (G4Cons*)(fmotherSolid); 245 return msol->GetDeltaPhiAngle(); 246 } 247 248 //-------------------------------------------------------------------------- 249 void 250 G4ParameterisationConsPhi:: 251 ComputeTransformation( const G4int copyNo, G4VPhysicalVolume *physVol ) const 252 { 253 //----- translation 254 G4ThreeVector origin(0.,0.,0.); 255 //----- set translation 256 physVol->SetTranslation( origin ); 257 258 //----- calculate rotation matrix (so that all volumes point to the centre) 259 G4double posi = foffset + copyNo*fwidth; 260 261 #ifdef G4DIVDEBUG 262 if( verbose >= 2 ) 263 { 264 G4cout << " G4ParameterisationConsPhi - position: " << posi/CLHEP::deg << G4endl 265 << " Origin: " << origin << " copyNo: " << copyNo 266 << " - foffset: " << foffset/CLHEP::deg 267 << " - fwidth: " << fwidth/CLHEP::deg << G4endl 268 << " - Axis: " << faxis << G4endl; 269 } 270 #endif 271 272 ChangeRotMatrix( physVol, -posi ); 273 } 274 275 //-------------------------------------------------------------------------- 276 void 277 G4ParameterisationConsPhi:: 278 ComputeDimensions( G4Cons& cons, const G4int, 279 const G4VPhysicalVolume* ) const 280 { 281 auto msol = (G4Cons*)(fmotherSolid); 282 283 G4double pRMin1 = msol->GetInnerRadiusMinusZ(); 284 G4double pRMax1 = msol->GetOuterRadiusMinusZ(); 285 G4double pRMin2 = msol->GetInnerRadiusPlusZ(); 286 G4double pRMax2 = msol->GetOuterRadiusPlusZ(); 287 G4double pDz = msol->GetZHalfLength(); 288 289 //- already rotated double pSPhi = foffset + copyNo*fwidth; 290 G4double pSPhi = foffset + msol->GetStartPhiAngle() + fhgap; 291 G4double pDPhi = fwidth - 2.*fhgap; 292 293 cons.SetInnerRadiusMinusZ( pRMin1 ); 294 cons.SetOuterRadiusMinusZ( pRMax1 ); 295 cons.SetInnerRadiusPlusZ( pRMin2 ); 296 cons.SetOuterRadiusPlusZ( pRMax2 ); 297 cons.SetZHalfLength( pDz ); 298 cons.SetStartPhiAngle( pSPhi, false ); 299 cons.SetDeltaPhiAngle( pDPhi ); 300 301 #ifdef G4DIVDEBUG 302 if( verbose >= 2 ) 303 { 304 G4cout << " G4ParameterisationConsPhi::ComputeDimensions" << G4endl 305 << " pSPhi: " << pSPhi << " - pDPhi: " << pDPhi << G4endl; 306 if( verbose >= 4 ) cons.DumpInfo(); 307 } 308 #endif 309 } 310 311 //-------------------------------------------------------------------------- 312 G4ParameterisationConsZ:: 313 G4ParameterisationConsZ( EAxis axis, G4int nDiv, 314 G4double width, G4double offset, 315 G4VSolid* msolid, DivisionType divType ) 316 : G4VParameterisationCons( axis, nDiv, width, offset, msolid, divType ) 317 { 318 CheckParametersValidity(); 319 SetType( "DivisionConsZ" ); 320 321 auto msol = (G4Cons*)(fmotherSolid); 322 if( divType == DivWIDTH ) 323 { 324 fnDiv = CalculateNDiv( 2*msol->GetZHalfLength(), width, offset ); 325 } 326 else if( divType == DivNDIV ) 327 { 328 fwidth = CalculateWidth( 2*msol->GetZHalfLength(), nDiv, offset ); 329 } 330 331 #ifdef G4DIVDEBUG 332 if( verbose >= 1 ) 333 { 334 G4cout << " G4ParameterisationConsZ: # divisions " << fnDiv << " = " 335 << nDiv << G4endl 336 << " Offset " << foffset << " = " << offset << G4endl 337 << " Width " << fwidth << " = " << width << G4endl 338 << " - Axis: " << faxis << G4endl; 339 } 340 #endif 341 } 342 343 //-------------------------------------------------------------------------- 344 G4ParameterisationConsZ::~G4ParameterisationConsZ() = default; 345 346 //------------------------------------------------------------------------ 347 G4double G4ParameterisationConsZ::GetMaxParameter() const 348 { 349 auto msol = (G4Cons*)(fmotherSolid); 350 return 2*msol->GetZHalfLength(); 351 } 352 353 //-------------------------------------------------------------------------- 354 void 355 G4ParameterisationConsZ:: 356 ComputeTransformation( const G4int copyNo, G4VPhysicalVolume* physVol ) const 357 { 358 //----- set translation: along Z axis 359 auto motherCons = (G4Cons*)(GetMotherSolid()); 360 G4double posi = - motherCons->GetZHalfLength() + OffsetZ() 361 + fwidth/2 + copyNo*fwidth; 362 G4ThreeVector origin(0.,0.,posi); 363 physVol->SetTranslation( origin ); 364 365 //----- calculate rotation matrix: unit 366 367 #ifdef G4DIVDEBUG 368 if( verbose >= 2 ) 369 { 370 G4cout << " G4ParameterisationConsZ::ComputeTransformation()" << G4endl 371 << " Origin: " << origin << " - copyNo: " << copyNo << G4endl 372 << " foffset: " << foffset << " - fwidth: " << fwidth 373 << G4endl; 374 } 375 #endif 376 377 ChangeRotMatrix( physVol ); 378 } 379 380 381 //-------------------------------------------------------------------------- 382 void 383 G4ParameterisationConsZ:: 384 ComputeDimensions( G4Cons& cons, const G4int copyNo, 385 const G4VPhysicalVolume* ) const 386 { 387 auto msol = (G4Cons*)(fmotherSolid); 388 389 G4double mHalfLength = msol->GetZHalfLength() - fhgap; 390 G4double aRInner = (msol->GetInnerRadiusPlusZ() 391 - msol->GetInnerRadiusMinusZ()) / (2*mHalfLength); 392 G4double bRInner = (msol->GetInnerRadiusPlusZ() 393 + msol->GetInnerRadiusMinusZ()) / 2; 394 G4double aROuter = (msol->GetOuterRadiusPlusZ() 395 - msol->GetOuterRadiusMinusZ()) / (2*mHalfLength); 396 G4double bROuter = (msol->GetOuterRadiusPlusZ() 397 + msol->GetOuterRadiusMinusZ()) / 2; 398 G4double xMinusZ = -mHalfLength + OffsetZ() + fwidth*copyNo + fhgap; 399 G4double xPlusZ = -mHalfLength + OffsetZ() + fwidth*(copyNo+1) - fhgap; 400 cons.SetInnerRadiusMinusZ( aRInner * xMinusZ + bRInner ); 401 cons.SetOuterRadiusMinusZ( aROuter * xMinusZ + bROuter ); 402 cons.SetInnerRadiusPlusZ( aRInner * xPlusZ + bRInner ); 403 cons.SetOuterRadiusPlusZ( aROuter * xPlusZ + bROuter ); 404 405 G4double pDz = fwidth / 2. - fhgap; 406 G4double pSPhi = msol->GetStartPhiAngle(); 407 G4double pDPhi = msol->GetDeltaPhiAngle(); 408 409 cons.SetZHalfLength( pDz ); 410 cons.SetStartPhiAngle( pSPhi, false ); 411 cons.SetDeltaPhiAngle( pDPhi ); 412 413 #ifdef G4DIVDEBUG 414 if( verbose >= 2 ) 415 { 416 G4cout << " G4ParameterisationConsZ::ComputeDimensions()" << G4endl 417 << " pDz: " << pDz << G4endl; 418 if( verbose >= 4 ) cons.DumpInfo(); 419 } 420 #endif 421 422 } 423