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 // Code developed by 27 // Silvia Pozzi (1), silvia.pozzi@iss.it 28 // Barbara Caccia (1), barbara.caccia@iss.it 29 // Carlo Mancini Terracciano (2), carlo.mancin 30 // (1) Istituto Superiore di Sanita' and INFN 31 // (2) Univ. La Sapienza and INFN Roma, Italy 32 33 #include "DetectorConstruction.hh" 34 #include "DetectorMessenger.hh" 35 36 #include <G4LogicalVolume.hh> 37 #include <G4PVPlacement.hh> 38 #include <G4PVReplica.hh> 39 #include <G4ProductionCuts.hh> 40 #include <G4NistManager.hh> 41 #include <G4SystemOfUnits.hh> 42 #include <G4VisAttributes.hh> 43 #include <G4Box.hh> 44 #include <G4VSolid.hh> 45 #include <G4Cons.hh> 46 #include <G4Tubs.hh> 47 #include <G4BooleanSolid.hh> 48 #include <G4SubtractionSolid.hh> 49 #include <G4SDManager.hh> 50 #include <G4GlobalMagFieldMessenger.hh> 51 #include <G4UnionSolid.hh> 52 #include <G4SubtractionSolid.hh> 53 #include <G4RunManager.hh> 54 #include <G4GeometryManager.hh> 55 #include <G4PhysicalVolumeStore.hh> 56 #include <G4LogicalVolumeStore.hh> 57 #include <G4SolidStore.hh> 58 #include <G4Transform3D.hh> 59 #include <CLHEP/Vector/Rotation.h> 60 61 62 G4VPhysicalVolume* DetectorConstruction::Const 63 { 64 detectorMessenger = new DetectorMessenger(th 65 66 //// Saturn 43 67 G4Material* air = G4NistManager::Instance()- 68 69 G4double worldSizeX = 3.5 * m; 70 G4double worldSizeY = 3.5 * m; 71 G4double worldSizeZ = 3.5 * m; 72 73 G4VSolid* worldBox = new G4Box("world", worl 74 75 G4LogicalVolume* worldLog = new G4LogicalVol 76 77 G4VisAttributes* visAttr = new G4VisAttribut 78 visAttr->SetVisibility(false); 79 worldLog->SetVisAttributes(visAttr); 80 81 world_phys = new G4PVPlacement(nullptr, {}, 82 83 ConstructMaterials(); 84 ConstructAccelerator(); 85 ConstructPhantom(); 86 87 return world_phys; 88 } 89 90 void DetectorConstruction::ConstructPhantom() 91 { 92 G4Material* water = G4NistManager::Instance( 93 94 G4ThreeVector phantomHalfDim = {phantomSideD 95 96 G4VSolid* Phantom = new G4Box("Phantom", pha 97 G4LogicalVolume* Phantom_log = new G4Logical 98 99 G4VSolid* spessore = new G4Box("Spessore", p 100 G4LogicalVolume* spessore_log = new G4Logica 101 new G4PVPlacement(nullptr, {0.,0., 1.25*mm}, 102 103 // axes origin on the face of the phantom 104 G4ThreeVector phantomPos = {0., 0., phantomH 105 106 phantom_phys = new G4PVPlacement(nullptr, ph 107 "phantom_phys", Phantom_log, world_ph 108 109 PhantomSegmentation(Phantom_log); 110 111 // Region for cuts 112 G4Region *regVol= new G4Region("fullWaterPha 113 G4ProductionCuts* cuts = new G4ProductionCut 114 cuts->SetProductionCut(0.1*mm); 115 regVol->SetProductionCuts(cuts); 116 117 Phantom_log -> SetRegion(regVol); 118 regVol->AddRootLogicalVolume(Phantom_log); 119 120 //--------- Visualization attributes ------- 121 /// Phantom 122 G4VisAttributes* simpleH2OVisAtt= new G4VisA 123 simpleH2OVisAtt->SetVisibility(true); 124 simpleH2OVisAtt->SetForceSolid(false); 125 Phantom_log->SetVisAttributes(simpleH2OVisAt 126 127 } 128 129 G4Material* DetectorConstruction::GetMaterial( 130 { 131 G4Material* material = nullptr; 132 133 if (materialName == "kapton") 134 { 135 material = mat_Kapton; 136 } 137 else if (materialName == "steel") 138 { 139 material = mat_Ssteel; 140 } 141 else if (materialName == "xc10") 142 { 143 material = mat_XC10; 144 } 145 else if (materialName == "wnicu") 146 { 147 material = mat_WNICU; 148 } 149 else 150 { 151 G4cerr << "* DetectorConstruction::GetMa 152 materialName << " not found" << G4en 153 } 154 155 return material; 156 } 157 158 void DetectorConstruction::ConstructMaterials( 159 { 160 G4NistManager* nist = G4NistManager::Insta 161 162 G4Material* el_H = nist -> FindOrBuildMater 163 G4Material* el_C = nist -> FindOrBuildMater 164 G4Material* el_N = nist -> FindOrBuildMater 165 G4Material* el_O = nist -> FindOrBuildMater 166 G4Material* el_Si = nist -> FindOrBuildMater 167 G4Material* el_Cr = nist -> FindOrBuildMater 168 G4Material* el_Mn = nist -> FindOrBuildMater 169 G4Material* el_Fe = nist -> FindOrBuildMater 170 G4Material* el_Ni = nist -> FindOrBuildMater 171 G4Material* el_Cu = nist -> FindOrBuildMater 172 G4Material* el_W = nist -> FindOrBuildMater 173 174 mat_Ssteel = new G4Material("StainlessSteel" 175 mat_Ssteel -> AddMaterial(el_Fe, 0.6898); 176 mat_Ssteel -> AddMaterial(el_Cr, 0.18); 177 mat_Ssteel -> AddMaterial(el_Ni, 0.10); 178 mat_Ssteel -> AddMaterial(el_Mn, 0.02); 179 mat_Ssteel -> AddMaterial(el_Si, 0.01); 180 mat_Ssteel -> AddMaterial(el_C, 0.0002); 181 182 mat_XC10 = new G4Material("CARBON_STEEL", 7. 183 mat_XC10 -> AddMaterial(el_Fe, 0.993); 184 mat_XC10 -> AddMaterial(el_Mn, 0.006); 185 mat_XC10 -> AddMaterial(el_C, 0.001); 186 187 mat_WNICU = new G4Material("Denal(WNICU)", 1 188 mat_WNICU -> AddMaterial(el_W, 0.905); 189 mat_WNICU -> AddMaterial(el_Ni, 0.07); 190 mat_WNICU -> AddMaterial(el_Cu, 0.025); 191 192 mat_Kapton = new G4Material("Kapton", 1.42*g 193 mat_Kapton -> AddMaterial(el_C, 69.1133 *per 194 mat_Kapton -> AddMaterial(el_O, 20.9235 *per 195 mat_Kapton -> AddMaterial(el_N, 7.3270 *perC 196 mat_Kapton -> AddMaterial(el_H, 2.6362 *perC 197 198 } 199 200 void DetectorConstruction::ConstructAccelerato 201 { 202 G4bool checkOverlap = true; 203 204 G4Material* Vacuum = G4NistManager::Instance 205 206 accHalfSize = {600.*mm, 600.*mm, 600.*mm}; 207 G4ThreeVector initialCentre = {0.*mm, 0.*mm, 208 209 G4Box* accWorldBox = new G4Box("accWorld", a 210 G4LogicalVolume* accWorldLV = new G4LogicalV 211 accWorld_phys = new G4PVPlacement(nullptr, i 212 213 G4VisAttributes* simpleAlSVisAtt = new G4Vis 214 simpleAlSVisAtt -> SetVisibility(false); 215 accWorldLV -> SetVisAttributes(simpleAlSVisA 216 217 ConstructTarget(); 218 ConstructPrimaryCollimator(); 219 ConstructVacuumWindow(); 220 ConstructFlatteningFilter(); 221 ConstructIonizationChamber(); 222 ConstructJawsX(); 223 ConstructJawsY(); 224 } 225 226 void DetectorConstruction::ConstructTarget() 227 { 228 G4bool checkOverlap = true; 229 230 G4ThreeVector targetCentre, boxHalfSide, tub 231 G4double tubeRadius, tubeHalfLengthZ; 232 targetCentre.set(0, 0, -3.5*mm); 233 boxHalfSide.set(5.*mm, 5.*mm, 7.5*mm); 234 235 tubeRadius = 3.*mm; 236 tubeHalfLengthZ = 5.5*mm; // 11 mm tot 237 238 G4Material* diskMaterial = G4NistManager::In 239 G4Material* targetMaterial = G4NistManager:: 240 241 // Physical volumes 242 G4Box* box = new G4Box("targetBox", boxHalfS 243 G4Tubs* tube = new G4Tubs("targetTube", 0.*m 244 245 G4SubtractionSolid* BoxMinusTube = new G4Sub 246 G4LogicalVolume* BoxMinusTubeLV = new G4Logi 247 new G4PVPlacement(0, targetCentre, "TargetPV 248 249 G4Tubs* disk = new G4Tubs("targetDisk", 0.*m 250 G4LogicalVolume* diskLV = new G4LogicalVolum 251 new G4PVPlacement(0, {0.,0.,-15.*mm}, "diskP 252 253 // Visualization 254 G4VisAttributes* simpleAlSVisAtt; 255 simpleAlSVisAtt = new G4VisAttributes(G4Colo 256 simpleAlSVisAtt -> SetVisibility(true); 257 simpleAlSVisAtt -> SetForceSolid(true); 258 BoxMinusTubeLV -> SetVisAttributes(simpleAl 259 diskLV -> SetVisAttributes(simpleAlSVisAtt); 260 261 // Region for cuts 262 G4Region *regVol; 263 regVol = new G4Region("TargetR"); 264 G4ProductionCuts* cuts = new G4ProductionCut 265 cuts -> SetProductionCut(1.*mm); 266 regVol -> SetProductionCuts(cuts); 267 BoxMinusTubeLV -> SetRegion(regVol); 268 regVol -> AddRootLogicalVolume(BoxMinusTubeL 269 diskLV -> SetRegion(regVol); 270 regVol -> AddRootLogicalVolume(diskLV); 271 } 272 273 void DetectorConstruction::ConstructPrimaryCol 274 { 275 G4bool checkOverlap = true; 276 277 G4double bigTube1HalfLengthZ, bigTube1Radius 278 G4double bigTube2HalfLengthZ, bigTube2Radius 279 G4double bigTube3HalfLengthZ, bigTube3Radius 280 281 G4double tube1Radius, tube1HalfLengthZ, tube 282 283 G4double cone1RadiusMin, cone1RadiusMax, con 284 G4double cone2RadiusMin, cone2RadiusMax, con 285 G4double cone3RadiusMin, cone3RadiusMax, con 286 287 // upper principal tube = bigTube1 288 bigTube1Radius = 141./2.*mm; 289 bigTube1HalfLengthZ = 79./2. * mm; 290 bigTube1Z = 5. * mm + bigTube1HalfLe 291 tube1Radius = 34./2. * mm; 292 tube1HalfLengthZ = 15.86/2. * mm; 293 cone1RadiusMin = 34./2.*mm; // upper cone 294 cone1RadiusMax = 54./2.*mm; 295 cone1HalfLengthZ = (bigTube1HalfLengthZ - t 296 297 // middle principal tube = bigTube2 298 bigTube2Radius = 141./2. * mm; 299 bigTube2HalfLengthZ = 57.5/2. * mm; 300 bigTube2Z = bigTube1Z + bigTube1Half 301 tube2Radius = 53./2. * mm; 302 tube2HalfLengthZ = 1.35/2. * mm; 303 cone2RadiusMin = 53./2.*mm; // middle con 304 cone2RadiusMax = 81./2.*mm; 305 cone2HalfLengthZ = bigTube2HalfLengthZ - tu 306 307 // lower principal tube = bigTube3 308 bigTube3Radius = 241./2.*mm; 309 bigTube3HalfLengthZ = 40.50/2.*mm; 310 bigTube3Z = bigTube2Z + bigTube2HalfLe 311 cone3RadiusMin = 88.06/2.*mm; // lower co 312 cone3RadiusMax = 109.76/2.*mm; 313 cone3HalfLengthZ = 40.50/2.*mm; 314 315 G4Material* upperTubeMaterial = GetMaterial( 316 G4Material* middleTubeMaterial = GetMaterial 317 G4Material* lowerTubeMaterial = G4NistManage 318 319 // Physical volumes 320 G4Tubs* bigTube1 = new G4Tubs("PriCollBigTub 321 G4Tubs* tube1 = new G4Tubs("PriCollTube1", 0 322 G4Cons* cone1 = new G4Cons("PriCollCone1", 0 323 324 G4Tubs* bigTube2 = new G4Tubs("PriCollBigTub 325 G4Tubs* tube2 = new G4Tubs("PriCollTube2", 0 326 G4Cons* cone2 = new G4Cons("PriCollCone2", 0 327 328 G4Tubs* bigTube3 = new G4Tubs("PriCollBigTub 329 G4Cons* cone3 = new G4Cons("PriCollCone3", 0 330 331 G4UnionSolid* tube1AndCone1 = new G4UnionSol 332 tube1, cone1, nullptr, {0., 0., (tube1 333 G4SubtractionSolid* bigTube1NotTube1AndCone1 334 "PriColltube1NotTube1AndCone1", bigTube1 335 {0., 0., (-bigTube1HalfLengthZ + tube1Ha 336 G4LogicalVolume* bigTube1NotTube1AndCone1LV 337 bigTube1NotTube1AndCone1, upperTubeMater 338 339 G4UnionSolid* tube2AndCone2 = new G4UnionSol 340 tube2, cone2, nullptr, {0., 0., (tube2 341 G4SubtractionSolid* bigTube2NotTube2AndCone2 342 "PriCollBigTube2NotTube2Cone2", bigTube2 343 {0., 0., (-bigTube2HalfLengthZ + tube2Ha 344 G4LogicalVolume* bigTube2NotTube2AndCone2LV 345 bigTube2NotTube2AndCone2, middleTubeMate 346 347 G4SubtractionSolid* bigTube3NotCone3 = new G 348 G4LogicalVolume* bigTube3NotCone3LV = new G 349 lowerTubeMaterial, "PriCollBigTube3NotCo 350 351 new G4PVPlacement(nullptr, {0, 0, bigTube1Z} 352 bigTube1NotTube1AndCone1LV, accWorld_phy 353 new G4PVPlacement(nullptr, {0, 0, bigTube2Z} 354 bigTube2NotTube2AndCone2LV, accWorld_phy 355 new G4PVPlacement(nullptr, {0, 0, bigTube3Z} 356 bigTube3NotCone3LV, accWorld_phys, false 357 358 // Visualization 359 G4VisAttributes* simpleAlSVisAtt1 = new G4Vi 360 G4VisAttributes* simpleAlSVisAtt2 = new G4Vi 361 G4VisAttributes* simpleAlSVisAtt3 = new G4Vi 362 simpleAlSVisAtt1 -> SetVisibility(true); 363 simpleAlSVisAtt1 -> SetForceSolid(false); 364 simpleAlSVisAtt1 -> SetForceWireframe(true); 365 bigTube1NotTube1AndCone1LV -> SetVisAttribut 366 simpleAlSVisAtt2 -> SetVisibility(true); 367 simpleAlSVisAtt2 -> SetForceSolid(false); 368 simpleAlSVisAtt2 -> SetForceWireframe(true); 369 bigTube2NotTube2AndCone2LV -> SetVisAttribut 370 simpleAlSVisAtt3 -> SetVisibility(true); 371 simpleAlSVisAtt3 -> SetForceSolid(false); 372 simpleAlSVisAtt3 -> SetForceWireframe(true); 373 bigTube3NotCone3LV -> SetVisAttributes(simpl 374 375 // Region for cuts 376 G4Region* regVol = new G4Region("primaryColl 377 G4ProductionCuts* cuts = new G4ProductionCut 378 cuts->SetProductionCut(1.*cm); 379 regVol->SetProductionCuts(cuts); 380 bigTube1NotTube1AndCone1LV -> SetRegion(regV 381 regVol -> AddRootLogicalVolume(bigTube1NotTu 382 bigTube2NotTube2AndCone2LV -> SetRegion(regV 383 regVol -> AddRootLogicalVolume(bigTube2NotTu 384 bigTube3NotCone3LV -> SetRegion(regVol); 385 regVol -> AddRootLogicalVolume(bigTube3NotCo 386 } 387 388 void DetectorConstruction::ConstructFlattening 389 { 390 G4double tube1HalfLengthZ; //tube1Radius, 391 G4double cone1RadiusMin, cone1RadiusMax, con 392 G4double cone2RadiusMin, cone2RadiusMax, con 393 G4double cone3RadiusMin, cone3RadiusMax, con 394 395 tubeFFRadius = 108./2.*mm; // top principa 396 tube1HalfLengthZ = 7.5/2.*mm; 397 398 cone1RadiusMin = 54./2.*mm; // top cone 399 cone1RadiusMax = 76./2.*mm; 400 cone1HalfLengthZ = 13.7/2.*mm; 401 402 cone2RadiusMin = 8./2.*mm; // mid cone 403 cone2RadiusMax = 54./2.*mm; 404 cone2HalfLengthZ = (44.3-13.7)/2.*mm; 405 406 cone3RadiusMin = 0.000001*mm; // bottom co 407 cone3RadiusMax = 8./2.*mm; 408 cone3HalfLengthZ = (46.8-44.3)/2.*mm; 409 tubeFFFirstFaceZ = 149.50*mm; 410 411 ffZ = tubeFFFirstFaceZ + tube1HalfLengthZ; / 412 413 G4Material* ffMaterial = GetMaterial("steel" 414 415 // Physical volumes 416 G4Tubs* tube1 = new G4Tubs("ffTube1", 0.*mm, 417 G4Cons* cone1 = new G4Cons("ffCone1", 0., co 418 G4Cons* cone2 = new G4Cons("ffCone2", 0., co 419 G4Cons* cone3 = new G4Cons("ffCone3", 0., co 420 421 G4double pos = (cone1HalfLengthZ - tube1Half 422 G4UnionSolid *tube1AndCone1 = new G4UnionSol 423 cone1, 0, {0., 0., pos}); 424 425 pos += cone1HalfLengthZ + cone2HalfLengthZ; 426 G4UnionSolid* tubeCone1AndCone2 = new G4Unio 427 tube1AndCone1, cone2, 0, {0., 0., pos}); 428 429 pos += cone2HalfLengthZ + cone3HalfLengthZ; 430 G4UnionSolid* tubeCone12AndCone3 = new G4Uni 431 tubeCone1AndCone2, cone3, 0, {0., 0., po 432 433 G4LogicalVolume* ffLV = new G4LogicalVolume( 434 435 new G4PVPlacement(0, {0,0,ffZ}, "ffPV", ffLV 436 437 // Visualization 438 G4VisAttributes* simpleAlSVisAtt= new G4VisA 439 simpleAlSVisAtt->SetVisibility(true); 440 simpleAlSVisAtt->SetForceSolid(true); 441 ffLV->SetVisAttributes(simpleAlSVisAtt); 442 443 // Region for cuts 444 G4Region* regVol; 445 regVol= new G4Region("ffR"); 446 G4ProductionCuts* cuts = new G4ProductionCut 447 cuts->SetProductionCut(0.1*cm); 448 regVol->SetProductionCuts(cuts); 449 ffLV->SetRegion(regVol); 450 regVol->AddRootLogicalVolume(ffLV); 451 } 452 453 void DetectorConstruction::ConstructIonization 454 { 455 G4bool checkOverlap = true; 456 457 G4double tubeRadius, tubeA1Z, tubeA2Z, tubeA 458 G4double kaptonThickness, AlThickness1, AlTh 459 460 kaptonThickness = 0.025/2. *mm; 461 AlThickness1 = 1.6e-5/2. *cm; 462 AlThickness8 = 8e-6/2. *cm; 463 464 tubeRadius = 110./2.*mm; // upper principal 465 tubeA1Z = (202.5 + AlThickness8)*mm; 466 tubeA2Z = (204.5 + AlThickness1)*mm; 467 tubeA3Z = (206.5 + AlThickness8)*mm; 468 tubeA4Z = (207.5 + AlThickness8)*mm; 469 tubeA5Z = (209.5 + AlThickness1)*mm; 470 tubeA6Z = (211.5 + AlThickness8)*mm; 471 472 G4double pos1 = (AlThickness1 + kaptonThickn 473 G4double pos8 = (AlThickness8 + kaptonThickn 474 475 G4Material* el_Al = G4NistManager::Instance( 476 G4Material* Kapton = GetMaterial("kapton"); 477 478 // Physical volumes 479 G4Tubs* tubeAl8 = new G4Tubs("ICTube1A", 0.* 480 G4Tubs* tubeK = new G4Tubs("ICTube1B", 0.*mm 481 G4Tubs* tubeAl1 = new G4Tubs("ICTube2A", 0.* 482 483 G4LogicalVolume* IC_Al8_LV = new G4LogicalVo 484 G4LogicalVolume* IC_Al1_LV = new G4LogicalVo 485 G4LogicalVolume* IC_K_LV = new G4LogicalVo 486 487 new G4PVPlacement(0, {0,0,tubeA1Z}, "IC_AL8_ 488 new G4PVPlacement(0, {0,0,tubeA1Z + pos8}, " 489 490 new G4PVPlacement(0, {0,0,tubeA2Z}, "IC_AL1_ 491 new G4PVPlacement(0, {0,0,tubeA2Z + pos1}, " 492 493 new G4PVPlacement(0, {0,0,tubeA3Z}, "IC_AL8_ 494 new G4PVPlacement(0, {0,0,tubeA3Z + pos8}, " 495 496 new G4PVPlacement(0, {0,0,tubeA4Z}, "IC_AL8_ 497 new G4PVPlacement(0, {0,0,tubeA4Z + pos8}, " 498 499 new G4PVPlacement(0, {0,0,tubeA5Z}, "IC_AL1_ 500 new G4PVPlacement(0, {0,0,tubeA5Z + pos1}, " 501 502 new G4PVPlacement(0, {0,0,tubeA6Z}, "IC_AL8_ 503 new G4PVPlacement(0, {0,0,tubeA6Z + pos8}, " 504 505 // Visualization 506 G4VisAttributes* simpleAlSVisAttAl1 = new G4 507 simpleAlSVisAttAl1 -> SetVisibility(true); 508 simpleAlSVisAttAl1 -> SetForceSolid(true); 509 IC_Al1_LV -> SetVisAttributes(simpleAlSVisAt 510 511 G4VisAttributes* simpleAlSVisAttAl8 = new G4 512 simpleAlSVisAttAl8 -> SetVisibility(true); 513 simpleAlSVisAttAl8 -> SetForceSolid(true); 514 IC_Al8_LV -> SetVisAttributes(simpleAlSVisAt 515 516 G4VisAttributes* simpleAlSVisAttK = new G4Vi 517 simpleAlSVisAttK -> SetVisibility(true); 518 simpleAlSVisAttK -> SetForceSolid(true); 519 IC_K_LV -> SetVisAttributes(simpleAlSVisAttK 520 521 // Region for cuts 522 G4Region * regVol = new G4Region("IonChamber 523 G4ProductionCuts* cuts = new G4ProductionCut 524 cuts -> SetProductionCut(0.1*mm); 525 regVol -> SetProductionCuts(cuts); 526 527 IC_Al1_LV -> SetRegion(regVol); 528 regVol -> AddRootLogicalVolume(IC_Al1_LV); 529 IC_Al8_LV -> SetRegion(regVol); 530 regVol -> AddRootLogicalVolume(IC_Al8_LV); 531 IC_K_LV -> SetRegion(regVol); 532 regVol -> AddRootLogicalVolume(IC_K_LV); 533 } 534 535 void DetectorConstruction::ConstructVacuumWind 536 { 537 G4double tubeRadius, tubeHalfLengthZ, tubeZ; 538 539 tubeRadius = 116.53/2.*mm; // upper princ 540 tubeHalfLengthZ = 2./2.*mm; 541 tubeZ = 215.75*mm + tubeHalfLengthZ; 542 543 G4Material* el_Al= G4NistManager::Instance() 544 545 // Physical volumes 546 G4Tubs* tube = new G4Tubs("VacWindowTube", 0 547 G4LogicalVolume* tubeLV = new G4LogicalVolum 548 549 new G4PVPlacement(0, {0,0,tubeZ}, "VacWindow 550 551 // Visualization 552 G4VisAttributes* simpleAlSVisAtt = new G4Vis 553 simpleAlSVisAtt -> SetVisibility(true); 554 simpleAlSVisAtt -> SetForceSolid(true); 555 tubeLV -> SetVisAttributes(simpleAlSVisAtt); 556 557 // Region for cuts 558 G4Region* regVol= new G4Region("VacWindowR") 559 G4ProductionCuts* cuts = new G4ProductionCut 560 cuts->SetProductionCut(0.1*cm); 561 regVol->SetProductionCuts(cuts); 562 tubeLV -> SetRegion(regVol); 563 regVol -> AddRootLogicalVolume(tubeLV); 564 } 565 void DetectorConstruction::ConstructJawsX() 566 { 567 G4bool checkOverlap = false; 568 569 G4double boxSide = 101.*mm; 570 G4double box1HalfLengthZ = 3.*mm; 571 G4double box2HalfLengthZ = 35.*mm; 572 G4double box3HalfLengthZ = 35.*mm; 573 G4double box4HalfLengthZ = 27.*mm; 574 G4double box5HalfLengthZ = 10.*mm; 575 G4double box6HalfLengthZ = 5.*mm; 576 577 G4Material* el_Pb = G4NistManager::Instance( 578 G4Material* XC10 = GetMaterial("xc10"); 579 G4Material* WNICU = GetMaterial("wnicu"); 580 581 G4ThreeVector boxJawXSide = {boxSide, boxSid 582 583 G4Box* boxJaw1X = new G4Box("Jaw1Xbox", boxJ 584 G4LogicalVolume* boxJaw1XLV = new G4LogicalV 585 G4NistManager::Instance()->FindOrBuildMa 586 G4Box* boxJaw2X = new G4Box("Jaw2Xbox", boxJ 587 G4LogicalVolume* boxJaw2XLV = new G4LogicalV 588 G4NistManager::Instance()->FindOrBuild 589 590 jaw1XInitialPos.setX(boxJawXSide.getX()/2.); 591 jaw1XInitialPos.setY(0.); 592 jaw1XInitialPos.setZ(275.5*mm + boxJawXSide. 593 594 G4RotationMatrix* cRotationJaw1X = new G4Rot 595 cRotationJaw1X -> rotateY(std::fabs(std::ata 596 597 G4ThreeVector position = *cRotationJaw1X * j 598 *cRotationJaw1X = cRotationJaw1X -> inverse( 599 600 boxJaw1X_phys = 601 new G4PVPlacement (cRotationJaw1X, positio 602 "Jaw1XPV", boxJaw1XLV, accWorld_phys, 603 604 jaw2XInitialPos.setX(-jaw1XInitialPos.getX() 605 jaw2XInitialPos.setY(jaw1XInitialPos.getY()) 606 jaw2XInitialPos.setZ(jaw1XInitialPos.getZ()) 607 608 G4RotationMatrix* cRotationJaw2X = new G4Rot 609 cRotationJaw2X -> rotateY(-std::fabs(std::at 610 611 position = *cRotationJaw2X * jaw2XInitialPos 612 *cRotationJaw2X = cRotationJaw2X -> inverse( 613 614 boxJaw2X_phys = 615 new G4PVPlacement (cRotationJaw2X, pos 616 "Jaw2XPV", boxJaw2XLV, accWorld_ph 617 618 // Physical volumes 619 G4Box* box1 = new G4Box("Jaws1XBox1", boxSid 620 G4Box* box2 = new G4Box("Jaws1XBox2", boxSid 621 G4Box* box3 = new G4Box("Jaws1XBox3", boxSid 622 G4Box* box4 = new G4Box("Jaws1XBox4", boxSid 623 G4Box* box5 = new G4Box("Jaws1XBox5", boxSid 624 G4Box* box6 = new G4Box("Jaws1XBox6", boxSid 625 G4LogicalVolume* box1LV1 = new G4LogicalVolu 626 G4LogicalVolume* box2LV1 = new G4LogicalVolu 627 G4LogicalVolume* box3LV1 = new G4LogicalVolu 628 G4LogicalVolume* box4LV1 = new G4LogicalVolu 629 G4LogicalVolume* box5LV1 = new G4LogicalVolu 630 G4LogicalVolume* box6LV1 = new G4LogicalVolu 631 632 G4LogicalVolume* box1LV2 = new G4LogicalVolu 633 G4LogicalVolume* box2LV2 = new G4LogicalVolu 634 G4LogicalVolume* box3LV2 = new G4LogicalVolu 635 G4LogicalVolume* box4LV2 = new G4LogicalVolu 636 G4LogicalVolume* box5LV2 = new G4LogicalVolu 637 G4LogicalVolume* box6LV2 = new G4LogicalVolu 638 639 G4ThreeVector centre = {0., 0., 0.}; 640 G4double zCentreCurrentBox = -boxJawXSide.ge 641 642 centre.setZ(zCentreCurrentBox); 643 new G4PVPlacement(nullptr, centre, "Jaws1XPV 644 new G4PVPlacement(nullptr, centre, "Jaws2XPV 645 646 zCentreCurrentBox += (box1HalfLengthZ+box2Ha 647 centre.setZ(zCentreCurrentBox); 648 new G4PVPlacement(nullptr, centre, "Jaws1XPV 649 new G4PVPlacement(nullptr, centre, "Jaws2XPV 650 651 zCentreCurrentBox += (box2HalfLengthZ+box3Ha 652 centre.setZ(zCentreCurrentBox); 653 new G4PVPlacement(nullptr, centre, "Jaws1XPV 654 new G4PVPlacement(nullptr, centre, "Jaws2XPV 655 656 zCentreCurrentBox += (box3HalfLengthZ+box4Ha 657 centre.setZ(zCentreCurrentBox); 658 new G4PVPlacement(nullptr, centre, "Jaws1XPV 659 new G4PVPlacement(nullptr, centre, "Jaws2XPV 660 661 zCentreCurrentBox += (box4HalfLengthZ+box5Ha 662 centre.setZ(zCentreCurrentBox); 663 new G4PVPlacement(nullptr, centre, "Jaws1XPV 664 new G4PVPlacement(nullptr, centre, "Jaws2XPV 665 666 zCentreCurrentBox += (box5HalfLengthZ+box6Ha 667 centre.setZ(zCentreCurrentBox); 668 new G4PVPlacement(nullptr, centre, "Jaws1XPV 669 new G4PVPlacement(nullptr, centre, "Jaws2XPV 670 671 // Visualization 672 G4VisAttributes* visAtt = new G4VisAttribute 673 visAtt -> SetVisibility(false); 674 visAtt -> SetForceSolid(false); 675 676 G4VisAttributes* simpleAlSVisAttPb = new G4V 677 simpleAlSVisAttPb -> SetVisibility(true); 678 simpleAlSVisAttPb -> SetForceSolid(true); 679 680 G4VisAttributes* simpleAlSVisAttXC10 = new G 681 simpleAlSVisAttXC10 -> SetVisibility(true); 682 simpleAlSVisAttXC10 ->SetForceSolid(true); 683 684 G4VisAttributes* simpleAlSVisAttWNICU = new 685 simpleAlSVisAttWNICU ->SetVisibility(true); 686 simpleAlSVisAttWNICU ->SetForceSolid(true); 687 688 box1LV1 -> SetVisAttributes(simpleAlSVisAttX 689 box2LV1 -> SetVisAttributes(simpleAlSVisAttP 690 box3LV1 -> SetVisAttributes(simpleAlSVisAttW 691 box4LV1 -> SetVisAttributes(simpleAlSVisAttP 692 box5LV1 -> SetVisAttributes(simpleAlSVisAttP 693 box6LV1 -> SetVisAttributes(simpleAlSVisAttW 694 695 box1LV2 -> SetVisAttributes(simpleAlSVisAttX 696 box2LV2 -> SetVisAttributes(simpleAlSVisAttP 697 box3LV2 -> SetVisAttributes(simpleAlSVisAttW 698 box4LV2 -> SetVisAttributes(simpleAlSVisAttP 699 box5LV2 -> SetVisAttributes(simpleAlSVisAttP 700 box6LV2 -> SetVisAttributes(simpleAlSVisAttW 701 702 boxJaw1XLV -> SetVisAttributes(visAtt); 703 boxJaw2XLV -> SetVisAttributes(visAtt); 704 705 // Region for cuts 706 G4Region* regVol = new G4Region("JawsXR"); 707 G4ProductionCuts* cuts = new G4ProductionCut 708 cuts -> SetProductionCut(2.*cm); 709 regVol -> SetProductionCuts(cuts); 710 box1LV1 -> SetRegion(regVol); 711 regVol -> AddRootLogicalVolume(box1LV1); 712 box2LV1 -> SetRegion(regVol); 713 regVol -> AddRootLogicalVolume(box2LV1); 714 box3LV1 -> SetRegion(regVol); 715 regVol -> AddRootLogicalVolume(box3LV1); 716 box4LV1 -> SetRegion(regVol); 717 regVol -> AddRootLogicalVolume(box4LV1); 718 box5LV1 -> SetRegion(regVol); 719 regVol -> AddRootLogicalVolume(box5LV1); 720 box6LV1 -> SetRegion(regVol); 721 regVol -> AddRootLogicalVolume(box6LV1); 722 723 box1LV2 -> SetRegion(regVol); 724 regVol -> AddRootLogicalVolume(box1LV2); 725 box2LV2 -> SetRegion(regVol); 726 regVol -> AddRootLogicalVolume(box2LV2); 727 box3LV2 -> SetRegion(regVol); 728 regVol -> AddRootLogicalVolume(box3LV2); 729 box4LV2 -> SetRegion(regVol); 730 regVol -> AddRootLogicalVolume(box4LV2); 731 box5LV2 -> SetRegion(regVol); 732 regVol -> AddRootLogicalVolume(box5LV2); 733 box6LV2 -> SetRegion(regVol); 734 regVol -> AddRootLogicalVolume(box6LV2); 735 } 736 737 void DetectorConstruction::ConstructJawsY() 738 { 739 G4bool checkOverlap = false; 740 741 G4double boxSide = 101.*mm; 742 G4double box1HalfLengthZ = 15.*mm; 743 G4double box2HalfLengthZ = 31.*mm; 744 G4double box3HalfLengthZ = 35.*mm; 745 G4double box4HalfLengthZ = 21.*mm; 746 747 G4Material* el_Pb = G4NistManager::Instance( 748 G4Material* XC10 = GetMaterial("xc10"); 749 G4Material* WNICU = GetMaterial("wnicu"); 750 751 G4ThreeVector boxJawYSide = {boxSide, boxSid 752 753 G4Box* boxJaw1Y = new G4Box("Jaw1Ybox", boxJ 754 G4LogicalVolume* boxJaw1YLV = new G4LogicalV 755 G4NistManager::Instance()->FindOrBuildMa 756 G4Box* boxJaw2Y = new G4Box("Jaw2Ybox", boxJ 757 G4LogicalVolume* boxJaw2YLV = new G4Logica 758 G4NistManager::Instance()->FindOrBuild 759 760 jaw1YInitialPos.setX(0.); 761 jaw1YInitialPos.setY(boxJawYSide.getX()/2.); 762 jaw1YInitialPos.setZ(248.5*mm + boxJawYSide. 763 764 G4RotationMatrix* cRotationJaw1Y = new G4Rot 765 cRotationJaw1Y -> rotateX(-std::fabs(std::at 766 767 G4ThreeVector position = *cRotationJaw1Y * j 768 *cRotationJaw1Y = cRotationJaw1Y -> inverse( 769 770 boxJaw1Y_phys = 771 new G4PVPlacement (cRotationJaw1Y, posit 772 "Jaw1YPV", boxJaw1YLV, accWorld_phys 773 774 jaw2YInitialPos.setX(jaw1YInitialPos.getX()) 775 jaw2YInitialPos.setY(-jaw1YInitialPos.getY() 776 jaw2YInitialPos.setZ(jaw1YInitialPos.getZ()) 777 778 G4RotationMatrix* cRotationJaw2Y = new G4Rot 779 cRotationJaw2Y -> rotateX(std::fabs(std::ata 780 781 position = *cRotationJaw2Y * jaw2YInitialPos 782 *cRotationJaw2Y = cRotationJaw2Y -> inverse( 783 784 boxJaw2Y_phys = 785 new G4PVPlacement (cRotationJaw2Y, posit 786 "Jaw2YPV", boxJaw2YLV, accWorld_phys 787 788 // Physical volumes 789 G4Box *box1 = new G4Box("Jaws1YBox1", boxSid 790 G4Box *box2 = new G4Box("Jaws1YBox2", boxSid 791 G4Box *box3 = new G4Box("Jaws1YBox3", boxSid 792 G4Box *box4 = new G4Box("Jaws1YBox4", boxSid 793 G4LogicalVolume* box1LV1 = new G4LogicalVolu 794 G4LogicalVolume* box2LV1 = new G4LogicalVolu 795 G4LogicalVolume* box3LV1 = new G4LogicalVolu 796 G4LogicalVolume* box4LV1 = new G4LogicalVolu 797 798 G4LogicalVolume* box1LV2 = new G4LogicalVolu 799 G4LogicalVolume* box2LV2 = new G4LogicalVolu 800 G4LogicalVolume* box3LV2 = new G4LogicalVolu 801 G4LogicalVolume* box4LV2 = new G4LogicalVolu 802 803 G4ThreeVector centre = {0., 0., 0.}; 804 G4double zCentreCurrentBox = -boxJawYSide.ge 805 806 centre.setZ(zCentreCurrentBox); 807 new G4PVPlacement(nullptr, centre, "Jaws1YPV 808 new G4PVPlacement(nullptr, centre, "Jaws2YPV 809 810 zCentreCurrentBox += (box1HalfLengthZ+box2Ha 811 centre.setZ(zCentreCurrentBox); 812 new G4PVPlacement(nullptr, centre, "Jaws1YPV 813 new G4PVPlacement(nullptr, centre, "Jaws2YPV 814 815 zCentreCurrentBox += (box2HalfLengthZ+box3Ha 816 centre.setZ(zCentreCurrentBox); 817 new G4PVPlacement(nullptr, centre, "Jaws1YPV 818 new G4PVPlacement(nullptr, centre, "Jaws2YPV 819 820 zCentreCurrentBox += (box3HalfLengthZ+box4Ha 821 centre.setZ(zCentreCurrentBox); 822 new G4PVPlacement(nullptr, centre, "Jaws1YPV 823 new G4PVPlacement(nullptr, centre, "Jaws2YPV 824 825 // Visualization 826 G4VisAttributes* visAtt = new G4VisAttribute 827 visAtt -> SetVisibility(false); 828 visAtt -> SetForceSolid(false); 829 830 G4VisAttributes* simpleAlSVisAttPb = new G4V 831 simpleAlSVisAttPb -> SetVisibility(true); 832 simpleAlSVisAttPb -> SetForceSolid(true); 833 834 G4VisAttributes* simpleAlSVisAttXC10 = new G 835 simpleAlSVisAttXC10 -> SetVisibility(true); 836 simpleAlSVisAttXC10 -> SetForceSolid(true); 837 838 G4VisAttributes* simpleAlSVisAttWNICU = new 839 simpleAlSVisAttWNICU -> SetVisibility(true); 840 simpleAlSVisAttWNICU -> SetForceSolid(true); 841 842 box1LV1 -> SetVisAttributes(simpleAlSVisAttW 843 box2LV1 -> SetVisAttributes(simpleAlSVisAttP 844 box3LV1 -> SetVisAttributes(simpleAlSVisAttW 845 box4LV1 -> SetVisAttributes(simpleAlSVisAttP 846 847 box1LV2 -> SetVisAttributes(simpleAlSVisAttW 848 box2LV2 -> SetVisAttributes(simpleAlSVisAttP 849 box3LV2 -> SetVisAttributes(simpleAlSVisAttW 850 box4LV2 -> SetVisAttributes(simpleAlSVisAttP 851 852 boxJaw1YLV -> SetVisAttributes(visAtt); 853 boxJaw2YLV -> SetVisAttributes(visAtt); 854 855 // Region for cuts 856 G4Region* regVol = new G4Region("JawsYR"); 857 G4ProductionCuts* cuts = new G4ProductionCut 858 cuts -> SetProductionCut(2.*cm); 859 regVol -> SetProductionCuts(cuts); 860 861 box1LV1 -> SetRegion(regVol); 862 regVol -> AddRootLogicalVolume(box1LV1); 863 box2LV1 -> SetRegion(regVol); 864 regVol -> AddRootLogicalVolume(box2LV1); 865 box3LV1 -> SetRegion(regVol); 866 regVol -> AddRootLogicalVolume(box3LV1); 867 box4LV1 -> SetRegion(regVol); 868 regVol -> AddRootLogicalVolume(box4LV1); 869 870 box1LV2 -> SetRegion(regVol); 871 regVol -> AddRootLogicalVolume(box1LV2); 872 box2LV2 -> SetRegion(regVol); 873 regVol -> AddRootLogicalVolume(box2LV2); 874 box3LV2 -> SetRegion(regVol); 875 regVol -> AddRootLogicalVolume(box3LV2); 876 box4LV2 -> SetRegion(regVol); 877 regVol -> AddRootLogicalVolume(box4LV2); 878 } 879 880 void DetectorConstruction::PhantomSegmentation 881 { 882 883 G4ThreeVector phantomHalfDim = {phantomSideD 884 885 nSideCells = CheckPhantomSegmentation(phanto 886 nDepthCells = CheckPhantomSegmentation(phant 887 888 G4Material* water = G4NistManager::Instance( 889 890 G4String yRepName("RepY"); 891 G4VSolid* solYRep = new G4Box(yRepName, phan 892 phantomHalfDim.getY()/nSideCells, phanto 893 G4LogicalVolume* logYRep = new G4LogicalVolu 894 new G4PVReplica(yRepName, logYRep, Phantom_l 895 2.*phantomHalfDim.getY()/nSideCells); 896 897 G4String xRepName("RepX"); 898 G4VSolid* solXRep = new G4Box(xRepName, phan 899 phantomHalfDim.getY()/nSideCells, phanto 900 G4LogicalVolume* logXRep = new G4LogicalVolu 901 new G4PVReplica(xRepName, logXRep, logYRep, 902 2.*phantomHalfDim.getX()/nSideCells); 903 904 G4String zVoxName("phantomSens"); 905 G4VSolid* solVoxel = new G4Box(zVoxName, pha 906 phantomHalfDim.getY()/nSideCells, phanto 907 908 LVPhantomSens = new G4LogicalVolume(solVoxel 909 new G4PVReplica(zVoxName, LVPhantomSens, log 910 2.*phantomHalfDim.getZ()/nDepthCells); 911 912 } 913 914 G4int DetectorConstruction::CheckPhantomSegmen 915 { 916 // I want an odd number of voxels 917 // check whether nCells is integer 918 if ( nCells != (int)nCells ) 919 { 920 nCells = std::ceil(nCells); 921 } 922 // check whether nCells is even 923 if ( std::fmod(nCells,2) == 0) 924 { 925 nCells += 1; 926 } 927 928 return nCells; 929 } 930 931 G4double DetectorConstruction::GetAccOriginPos 932 { 933 return sourceToSkinDistance; 934 } 935 936 G4double DetectorConstruction::GetFFilterRadiu 937 { 938 return tubeFFRadius; 939 } 940 941 G4double DetectorConstruction::GetFFilterZ() 942 { 943 return tubeFFFirstFaceZ; 944 } 945 946 G4double DetectorConstruction::GetVoxelDepthDi 947 { 948 return voxelDepthDim; 949 } 950 951 G4int DetectorConstruction::GetNumberSideCells 952 { 953 return nSideCells; 954 } 955 956 G4int DetectorConstruction::GetNumberDepthCell 957 { 958 return nDepthCells; 959 } 960 961 G4int DetectorConstruction::GetPhantomDepth() 962 { 963 return phantomSideDim; 964 } 965 966 void DetectorConstruction::SetJaws(G4double va 967 { 968 jawAperture = value; 969 } 970 971 void DetectorConstruction::SetTargetPosition(G 972 { 973 sourceToSkinDistance = value; 974 } 975 976 void DetectorConstruction::SetPhantomSide(G4do 977 { 978 phantomSideDim = value; 979 } 980 981 void DetectorConstruction::SetVoxelSide(G4doub 982 { 983 voxelSideDim = value; 984 } 985 986 void DetectorConstruction::SetVoxelDepth(G4dou 987 { 988 voxelDepthDim = value; 989 } 990 991 void DetectorConstruction::UpdateGeometry(G4St 992 { 993 994 if ( string == "fieldSide" ) 995 { 996 // 997 G4double halfFieldSide = value/2; 998 G4ThreeVector position; 999 G4RotationMatrix* rMatrix1X = new G4Rotati 1000 1001 rMatrix1X -> rotateY(std::fabs(std::atan( 1002 position = *rMatrix1X * jaw1XInitialPos; 1003 boxJaw1X_phys -> SetTranslation(position) 1004 *rMatrix1X = rMatrix1X -> inverse(); 1005 boxJaw1X_phys -> SetRotation(rMatrix1X); 1006 1007 G4RotationMatrix* rMatrix2X = new G4Rotat 1008 rMatrix2X -> rotateY(-std::fabs(std::atan 1009 1010 position = *rMatrix2X * jaw2XInitialPos; 1011 boxJaw2X_phys -> SetTranslation(position) 1012 *rMatrix2X = rMatrix2X -> inverse(); 1013 boxJaw2X_phys -> SetRotation(rMatrix2X); 1014 1015 G4RotationMatrix* rMatrix1Y = new G4Rotat 1016 rMatrix1Y -> rotateX(-std::fabs(std::atan 1017 position = *rMatrix1Y * jaw1YInitialPos; 1018 boxJaw1Y_phys -> SetTranslation(position) 1019 *rMatrix1Y = rMatrix1Y -> inverse(); 1020 boxJaw1Y_phys -> SetRotation(rMatrix1Y); 1021 1022 G4RotationMatrix* rMatrix2Y = new G4Rotat 1023 rMatrix2Y -> rotateX(std::fabs(std::atan( 1024 1025 position = *rMatrix2Y * jaw2YInitialPos; 1026 boxJaw2Y_phys -> SetTranslation(position) 1027 *rMatrix2Y = rMatrix2Y -> inverse(); 1028 boxJaw2Y_phys -> SetRotation(rMatrix2Y); 1029 } 1030 else if ( string == "sourceToSkinDistance") 1031 { 1032 accWorld_phys -> SetTranslation({0., 0., 1033 } 1034 else if (string == "phantomSide") 1035 { 1036 G4Box* box = (G4Box *) phantom_phys -> Ge 1037 1038 // change phantom side 1039 box -> SetXHalfLength(phantomSideDim/2.); 1040 box -> SetYHalfLength(phantomSideDim/2.); 1041 box -> SetZHalfLength(phantomSideDim/2.); 1042 phantom_phys -> SetTranslation({0., 0., p 1043 1044 // clear daughters and rebuild them 1045 phantom_phys -> GetLogicalVolume() -> Cle 1046 PhantomSegmentation(phantom_phys -> GetLo 1047 } 1048 else if (string == "voxelSide") 1049 { 1050 phantom_phys -> GetLogicalVolume() -> Cle 1051 PhantomSegmentation(phantom_phys -> GetLo 1052 } 1053 else if (string == "voxelDepth") 1054 { 1055 phantom_phys -> GetLogicalVolume() -> Cle 1056 PhantomSegmentation(phantom_phys -> GetLo 1057 } 1058 else 1059 G4cerr << "*** DetectorMessenger::UpdateG 1060 1061 G4RunManager::GetRunManager()->GeometryHasB 1062 } 1063 1064 1065