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 // G4ParameterisationTrd[X/Y/Z] implementation 27 // 28 // 26.05.03 - P.Arce, Initial version 29 // 08.04.04 - I.Hrivnacova, Implemented reflec 30 // 21.04.10 - M.Asai, Added gaps 31 // ------------------------------------------- 32 33 #include "G4ParameterisationTrd.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 "G4Trd.hh" 42 #include "G4Trap.hh" 43 44 //-------------------------------------------- 45 G4VParameterisationTrd:: 46 G4VParameterisationTrd( EAxis axis, G4int nDiv 47 G4double offset, G4VSo 48 DivisionType divType ) 49 : G4VDivisionParameterisation( axis, nDiv, 50 { 51 auto msol = (G4Trd*)(msolid); 52 if (msolid->GetEntityType() == "G4ReflectedS 53 { 54 // Get constituent solid 55 G4VSolid* mConstituentSolid 56 = ((G4ReflectedSolid*)msolid)->GetConst 57 msol = (G4Trd*)(mConstituentSolid); 58 59 // Create a new solid with inversed parame 60 auto newSolid 61 = new G4Trd(msol->GetName(), 62 msol->GetXHalfLength2(), mso 63 msol->GetYHalfLength2(), mso 64 msol->GetZHalfLength()); 65 msol = newSolid; 66 fmotherSolid = newSolid; 67 fReflectedSolid = true; 68 fDeleteSolid = true; 69 } 70 } 71 72 //-------------------------------------------- 73 G4VParameterisationTrd::~G4VParameterisationTr 74 75 //-------------------------------------------- 76 G4ParameterisationTrdX:: 77 G4ParameterisationTrdX( EAxis axis, G4int nDiv 78 G4double width, G4doub 79 G4VSolid* msolid, Divi 80 : G4VParameterisationTrd( axis, nDiv, width 81 { 82 CheckParametersValidity(); 83 SetType( "DivisionTrdX" ); 84 85 auto msol = (G4Trd*)(fmotherSolid); 86 if( divType == DivWIDTH ) 87 { 88 fnDiv = CalculateNDiv( msol->GetXHalfLengt 89 width, offset ); 90 } 91 else if( divType == DivNDIV ) 92 { 93 fwidth = CalculateWidth( msol->GetXHalfLen 94 nDiv, offset ); 95 } 96 97 #ifdef G4DIVDEBUG 98 if( verbose >= 1 ) 99 { 100 G4cout << " G4ParameterisationTrdX - ## di 101 << nDiv << G4endl 102 << " Offset " << foffset << " = " < 103 << " Width " << fwidth << " = " << 104 } 105 #endif 106 107 G4double mpDx1 = msol->GetXHalfLength1(); 108 G4double mpDx2 = msol->GetXHalfLength2(); 109 if( std::fabs(mpDx1 - mpDx2) > kCarTolerance 110 { 111 bDivInTrap = true; 112 } 113 } 114 115 //-------------------------------------------- 116 G4ParameterisationTrdX::~G4ParameterisationTrd 117 = default; 118 119 //-------------------------------------------- 120 G4double G4ParameterisationTrdX::GetMaxParamet 121 { 122 auto msol = (G4Trd*)(fmotherSolid); 123 return (msol->GetXHalfLength1()+msol->GetXHa 124 } 125 126 //-------------------------------------------- 127 void 128 G4ParameterisationTrdX:: 129 ComputeTransformation( const G4int copyNo, 130 G4VPhysicalVolume *phys 131 { 132 auto msol = (G4Trd*)(fmotherSolid ); 133 G4double mdx = ( msol->GetXHalfLength1() + m 134 //----- translation 135 G4ThreeVector origin(0.,0.,0.); 136 G4double posi; 137 posi = -mdx + foffset + (copyNo+0.5)*fwidth; 138 if( faxis == kXAxis ) 139 { 140 origin.setX( posi ); 141 } 142 else 143 { 144 std::ostringstream message; 145 message << "Only axes along X are allowed 146 G4Exception("G4ParameterisationTrdX::Compu 147 "GeomDiv0002", FatalException, 148 } 149 150 #ifdef G4DIVDEBUG 151 if( verbose >= 2 ) 152 { 153 G4cout << std::setprecision(8) 154 << " G4ParameterisationTrdX::Comput 155 << copyNo << G4endl 156 << " Position: " << origin << " - A 157 } 158 #endif 159 160 //----- set translation 161 physVol->SetTranslation( origin ); 162 } 163 164 //-------------------------------------------- 165 void 166 G4ParameterisationTrdX:: 167 ComputeDimensions( G4Trd& trd, [[maybe_unused] 168 { 169 auto msol = (G4Trd*)(fmotherSolid); 170 G4double pDy1 = msol->GetYHalfLength1(); 171 G4double pDy2 = msol->GetYHalfLength2(); 172 G4double pDz = msol->GetZHalfLength(); 173 G4double pDx = fwidth/2. - fhgap; 174 175 trd.SetAllParameters ( pDx, pDx, pDy1, pDy2, 176 177 #ifdef G4DIVDEBUG 178 if( verbose >= 2 ) 179 { 180 G4cout << " G4ParameterisationTrdX::Comput 181 << copyNo << G4endl; 182 trd.DumpInfo(); 183 } 184 #endif 185 } 186 187 //-------------------------------------------- 188 void 189 G4ParameterisationTrdX::ComputeDimensions( G4T 190 con 191 { 192 auto msol = (G4Trd*)(fmotherSolid); 193 G4double pDy1 = msol->GetYHalfLength1(); 194 G4double pDy2 = msol->GetYHalfLength2(); 195 G4double pDz = msol->GetZHalfLength(); 196 //fwidth is at Z=0 197 G4double pDx1 = msol->GetXHalfLength1(); 198 G4double pDx2 = msol->GetXHalfLength2(); 199 // G4double pDxAVE = (pDx1+pDx2)/2.; 200 G4double xChangeRatio = (pDx2-pDx1) / (pDx2+ 201 G4double fWidChange = xChangeRatio*fwidth; 202 G4double fWid1 = fwidth - fWidChange; 203 G4double fWid2 = fwidth + fWidChange; 204 G4double fOffsetChange = xChangeRatio*foffse 205 G4double fOffset1 = foffset - fOffsetChange; 206 G4double fOffset2 = foffset + fOffsetChange; 207 G4double cxy1 = -pDx1+fOffset1 + (copyNo+0.5 208 G4double cxy2 = -pDx2+fOffset2 + (copyNo+0.5 209 G4double alp = std::atan( (cxy2-cxy1)/(pDz*2 210 211 pDx1 = fwidth/2. - fWidChange/2.; 212 pDx2 = fwidth/2. + fWidChange/2.; 213 214 215 trap.SetAllParameters ( pDz, 216 alp, 217 0., 218 pDy1, 219 pDx1, 220 pDx1, 221 0., 222 pDy2, 223 pDx2 - fhgap, 224 pDx2 - fhgap * pDx2/pDx1, 225 0.); 226 227 #ifdef G4DIVDEBUG 228 if( verbose >= 2 ) 229 { 230 G4cout << " G4ParameterisationTrdX::Comput 231 << copyNo << G4endl; 232 trap.DumpInfo(); 233 } 234 #endif 235 } 236 237 238 //-------------------------------------------- 239 G4ParameterisationTrdY:: 240 G4ParameterisationTrdY( EAxis axis, G4int nDiv 241 G4double width, G4doub 242 G4VSolid* msolid, Divi 243 : G4VParameterisationTrd( axis, nDiv, width, 244 { 245 CheckParametersValidity(); 246 SetType( "DivisionTrdY" ); 247 248 auto msol = (G4Trd*)(fmotherSolid); 249 if( divType == DivWIDTH ) 250 { 251 fnDiv = CalculateNDiv( 2*msol->GetYHalfLen 252 width, offset ); 253 } 254 else if( divType == DivNDIV ) 255 { 256 fwidth = CalculateWidth( 2*msol->GetYHalfL 257 nDiv, offset ); 258 } 259 260 #ifdef G4DIVDEBUG 261 if( verbose >= 1 ) 262 { 263 G4cout << " G4ParameterisationTrdY no div 264 << " = " << nDiv << G4endl 265 << " Offset " << foffset << " = " 266 << " width " << fwidth << " = " << 267 } 268 #endif 269 } 270 271 //-------------------------------------------- 272 G4ParameterisationTrdY::~G4ParameterisationTrd 273 274 //-------------------------------------------- 275 G4double G4ParameterisationTrdY::GetMaxParamet 276 { 277 auto msol = (G4Trd*)(fmotherSolid); 278 return (msol->GetYHalfLength1()+msol->GetYHa 279 } 280 281 //-------------------------------------------- 282 void 283 G4ParameterisationTrdY:: 284 ComputeTransformation( const G4int copyNo, G4V 285 { 286 auto msol = (G4Trd*)(fmotherSolid ); 287 G4double mdy = ( msol->GetYHalfLength1() + 288 289 //----- translation 290 G4ThreeVector origin(0.,0.,0.); 291 G4double posi = -mdy + foffset + (copyNo+0.5 292 293 if( faxis == kYAxis ) 294 { 295 origin.setY( posi ); 296 } 297 else 298 { 299 std::ostringstream message; 300 message << "Only axes along Y are allowed 301 G4Exception("G4ParameterisationTrdY::Compu 302 "GeomDiv0002", FatalException, 303 } 304 305 #ifdef G4DIVDEBUG 306 if( verbose >= 2 ) 307 { 308 G4cout << std::setprecision(8) 309 << " G4ParameterisationTrdY::Comput 310 << " pos " << origin << " rot mat " 311 } 312 #endif 313 314 //----- set translation 315 physVol->SetTranslation( origin ); 316 } 317 318 //-------------------------------------------- 319 void 320 G4ParameterisationTrdY:: 321 ComputeDimensions(G4Trd& trd, const G4int, con 322 { 323 //---- The division along Y of a Trd will re 324 //--- if Y at -Z and +Z are equal, else use 325 auto msol = (G4Trd*)(fmotherSolid); 326 327 G4double pDx1 = msol->GetXHalfLength1(); 328 G4double pDx2 = msol->GetXHalfLength2(); 329 G4double pDz = msol->GetZHalfLength(); 330 G4double pDy = fwidth/2. - fhgap; 331 332 trd.SetAllParameters ( pDx1, pDx2, pDy, pDy, 333 334 #ifdef G4DIVDEBUG 335 if( verbose >= 2 ) 336 { 337 G4cout << " G4ParameterisationTrdY::Comput 338 trd.DumpInfo(); 339 } 340 #endif 341 } 342 343 //-------------------------------------------- 344 void 345 G4ParameterisationTrdY::ComputeDimensions( G4T 346 con 347 { 348 auto msol = (G4Trd*)(fmotherSolid); 349 G4double pDx1 = msol->GetXHalfLength1(); 350 G4double pDx2 = msol->GetXHalfLength2(); 351 G4double pDz = msol->GetZHalfLength(); 352 //fwidth is at Z=0 353 G4double pDy1 = msol->GetYHalfLength1(); 354 G4double pDy2 = msol->GetYHalfLength2(); 355 // G4double pDxAVE = (pDy1+pDy2)/2.; 356 G4double yChangeRatio = (pDy2-pDy1) / (pDy2+ 357 G4double fWidChange = yChangeRatio*fwidth; 358 G4double fWid1 = fwidth - fWidChange; 359 G4double fWid2 = fwidth + fWidChange; 360 G4double fOffsetChange = yChangeRatio*foffse 361 G4double fOffset1 = foffset - fOffsetChange; 362 G4double fOffset2 = foffset + fOffsetChange; 363 G4double cyx1 = -pDy1+fOffset1 + (copyNo+0.5 364 G4double cyx2 = -pDy2+fOffset2 + (copyNo+0.5 365 G4double alp = std::atan( (cyx2-cyx1)/(pDz*2 366 367 pDy1 = fwidth/2. - fWidChange/2.; 368 pDy2 = fwidth/2. + fWidChange/2.; 369 370 371 trap.SetAllParameters ( pDz, 372 alp, 373 90*CLHEP::deg, 374 pDy1, 375 pDx1, 376 pDx1, 377 0., 378 pDy2, 379 pDx2 - fhgap, 380 pDx2 - fhgap * pDx2/pDx1, 381 0.); 382 383 #ifdef G4DIVDEBUG 384 if( verbose >= 2 ) 385 { 386 G4cout << " G4ParameterisationTrdY::Comput 387 << copyNo << G4endl; 388 trap.DumpInfo(); 389 } 390 #endif 391 } 392 393 //-------------------------------------------- 394 G4ParameterisationTrdZ:: 395 G4ParameterisationTrdZ( EAxis axis, G4int nDiv 396 G4double width, G4doub 397 G4VSolid* msolid, Divi 398 : G4VParameterisationTrd( axis, nDiv, width, 399 { 400 CheckParametersValidity(); 401 SetType( "DivTrdZ" ); 402 403 auto msol = (G4Trd*)(fmotherSolid); 404 if( divType == DivWIDTH ) 405 { 406 fnDiv = CalculateNDiv( 2*msol->GetZHalfLen 407 width, offset ); 408 } 409 else if( divType == DivNDIV ) 410 { 411 fwidth = CalculateWidth( 2*msol->GetZHalfL 412 nDiv, offset ); 413 } 414 415 #ifdef G4DIVDEBUG 416 if( verbose >= 1 ) 417 { 418 G4cout << " G4ParameterisationTrdZ no divi 419 << " = " << nDiv << G4endl 420 << " Offset " << foffset << " = " < 421 << " Width " << fwidth << " = " << 422 } 423 #endif 424 } 425 426 //-------------------------------------------- 427 G4ParameterisationTrdZ::~G4ParameterisationTrd 428 429 //-------------------------------------------- 430 G4double G4ParameterisationTrdZ::GetMaxParamet 431 { 432 auto msol = (G4Trd*)(fmotherSolid); 433 return 2*msol->GetZHalfLength(); 434 } 435 436 //-------------------------------------------- 437 void 438 G4ParameterisationTrdZ:: 439 ComputeTransformation(const G4int copyNo, G4VP 440 { 441 auto msol = (G4Trd*)(fmotherSolid ); 442 G4double mdz = msol->GetZHalfLength(); 443 444 //----- translation 445 G4ThreeVector origin(0.,0.,0.); 446 G4double posi = -mdz + OffsetZ() + (copyNo+0 447 if( faxis == kZAxis ) 448 { 449 origin.setZ( posi ); 450 } 451 else 452 { 453 std::ostringstream message; 454 message << "Only axes along Z are allowed 455 G4Exception("G4ParameterisationTrdZ::Compu 456 "GeomDiv0002", FatalException, 457 } 458 459 #ifdef G4DIVDEBUG 460 if( verbose >= 1 ) 461 { 462 G4cout << std::setprecision(8) << " G4Para 463 << copyNo << G4endl 464 << " Position: " << origin << " - O 465 << " - Width: " << fwidth << " Axis 466 } 467 #endif 468 469 //----- set translation 470 physVol->SetTranslation( origin ); 471 } 472 473 //-------------------------------------------- 474 void 475 G4ParameterisationTrdZ:: 476 ComputeDimensions(G4Trd& trd, const G4int copy 477 const G4VPhysicalVolume*) co 478 { 479 //---- The division along Z of a Trd will re 480 auto msol = (G4Trd*)(fmotherSolid); 481 482 G4double pDx1 = msol->GetXHalfLength1(); 483 G4double DDx = (msol->GetXHalfLength2() - ms 484 G4double pDy1 = msol->GetYHalfLength1(); 485 G4double DDy = (msol->GetYHalfLength2() - ms 486 G4double pDz = fwidth/2. - fhgap; 487 G4double zLength = 2*msol->GetZHalfLength(); 488 489 trd.SetAllParameters( pDx1+DDx*(OffsetZ()+co 490 pDx1+DDx*(OffsetZ()+(c 491 pDy1+DDy*(OffsetZ()+co 492 pDy1+DDy*(OffsetZ()+(c 493 pDz ); 494 495 #ifdef G4DIVDEBUG 496 if( verbose >= 1 ) 497 { 498 G4cout << " G4ParameterisationTrdZ::Comput 499 << " - Mother TRD " << G4endl; 500 msol->DumpInfo(); 501 G4cout << " - Parameterised TRD: " 502 << copyNo << G4endl; 503 trd.DumpInfo(); 504 } 505 #endif 506 } 507