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