9#include <ir/geometry/GeoFarBeamLineCreator.h>
11#include <geometry/Materials.h>
12#include <geometry/CreatorFactory.h>
13#include <geometry/utilities.h>
14#include <framework/gearbox/Unit.h>
15#include <ir/simulation/SensitiveDetector.h>
16#include <simulation/background/BkgSensitiveDetector.h>
19#include <boost/algorithm/string.hpp>
21#include <G4LogicalVolume.hh>
22#include <G4PVPlacement.hh>
28#include <G4Polycone.hh>
30#include <G4IntersectionSolid.hh>
31#include <G4SubtractionSolid.hh>
32#include <G4UserLimits.hh>
43 using namespace geometry;
51 geometry::CreatorFactory<GeoFarBeamLineCreator> GeoFarBeamLineFactory(
"FarBeamLineCreator");
79 map<string, FarBeamLineElement> elements;
91 std::string prep =
"TubeR.";
95 std::vector<double> TubeR_Z(TubeR_N);
96 std::vector<double> TubeR_R(TubeR_N);
97 std::vector<double> TubeR_r(TubeR_N);
99 for (
int i = 0; i < TubeR_N; ++i) {
100 ostringstream ossZID;
103 ostringstream ossRID;
106 ostringstream ossrID;
114 tubeR.
transform = G4Translate3D(0.0, 0.0, 0.0);
117 tubeR.
geo =
new G4Polycone(
"geo_TubeR_name", 0.0, 2 * M_PI, TubeR_N, &(TubeR_Z[0]), &(TubeR_r[0]), &(TubeR_R[0]));
121 elements[
"TubeR"] = tubeR;
133 std::vector<double> TubeL_Z(TubeL_N);
134 std::vector<double> TubeL_R(TubeL_N);
135 std::vector<double> TubeL_r(TubeL_N);
137 for (
int i = 0; i < TubeL_N; ++i) {
138 ostringstream ossZID;
141 ostringstream ossRID;
144 ostringstream ossrID;
152 tubeL.
transform = G4Translate3D(0.0, 0.0, 0.0);
155 tubeL.
geo =
new G4Polycone(
"geo_TubeL_name", 0.0, 2 * M_PI, TubeL_N, &(TubeL_Z[0]), &(TubeL_r[0]), &(TubeL_R[0]));
159 elements[
"TubeL"] = tubeL;
161 std::vector<double> zero_r(N, 0.);
163 std::vector<std::string> straightSections;
165 for (
const auto& name : straightSections) {
175 std::vector<double> Polycone_Z(N);
176 std::vector<double> Polycone_R(N);
177 std::vector<double> Polycone_r(N);
191 polycone.
transform = G4Translate3D(Polycone_X0, Polycone_Y0, Polycone_Z0);
193 if (Polycone_PHIYZ != 0)
200 string geo_polyconexx_name =
"geo_" + name +
"xx_name";
201 string geo_polyconex_name =
"geo_" + name +
"x_name";
202 string geo_polycone_name =
"geo_" + name +
"_name";
204 G4VSolid* geo_polyconexx(NULL), *geo_polycone(NULL);
208 geo_polyconexx =
new G4Polycone(geo_polyconexx_name, 0.0, 2 * M_PI, N, &(Polycone_Z[0]), &(zero_r[0]), &(Polycone_R[0]));
210 geo_polyconexx =
new G4Polycone(geo_polyconexx_name, 0.0, 2 * M_PI, N, &(Polycone_Z[0]), &(Polycone_r[0]), &(Polycone_R[0]));
211 else if (type ==
"pipe")
212 geo_polycone =
new G4Polycone(geo_polycone_name, 0.0, 2 * M_PI, N, &(Polycone_Z[0]), &(zero_r[0]), &(Polycone_R[0]));
214 geo_polycone =
new G4Polycone(geo_polycone_name, 0.0, 2 * M_PI, N, &(Polycone_Z[0]), &(Polycone_r[0]), &(Polycone_R[0]));
218 G4VSolid* geo_polyconex =
new G4SubtractionSolid(geo_polyconex_name, geo_polyconexx, elements[subtract].geo,
219 polycone.
transform.inverse()*elements[subtract].transform);
220 geo_polycone =
new G4IntersectionSolid(geo_polycone_name, geo_polyconex, elements[
intersect].geo,
222 }
else if (subtract !=
"")
223 geo_polycone =
new G4SubtractionSolid(geo_polycone_name, geo_polyconexx, elements[subtract].geo,
224 polycone.
transform.inverse()*elements[subtract].transform);
226 geo_polycone =
new G4IntersectionSolid(geo_polycone_name, geo_polyconexx, elements[
intersect].geo,
229 polycone.
geo = geo_polycone;
234 string logi_polycone_name =
"logi_" + name +
"_name";
235 G4LogicalVolume* logi_polycone =
new G4LogicalVolume(polycone.
geo, mat_polycone, logi_polycone_name);
236 setColor(*logi_polycone,
"#CC0000");
239 polycone.
logi = logi_polycone;
242 string phys_polycone_name =
"phys_" + name +
"_name";
243 new G4PVPlacement(polycone.
transform, logi_polycone, phys_polycone_name, &topVolume,
false, 0);
245 elements[name] = polycone;
248 for (
int i = 0; i < N; ++i)
249 sum += Polycone_r[i];
250 if (type ==
"pipe" && sum != 0) {
254 string nameVac = name +
"Vac";
257 string geo_vacuumxx_name =
"geo_" + nameVac +
"xx_name";
258 string geo_vacuum_name =
"geo_" + nameVac +
"_name";
260 G4VSolid* geo_vacuumxx, *geo_vacuum;
262 geo_vacuumxx =
new G4Polycone(geo_vacuumxx_name, 0.0, 2 * M_PI, N, &(Polycone_Z[0]), &(zero_r[0]), &(Polycone_r[0]));
263 geo_vacuum =
new G4IntersectionSolid(geo_vacuumxx_name, geo_vacuumxx, geo_polycone);
265 vacuum.
geo = geo_vacuum;
270 string logi_vacuum_name =
"logi_" + nameVac +
"_name";
271 G4LogicalVolume* logi_vacuum =
new G4LogicalVolume(vacuum.
geo, mat_vacuum, logi_vacuum_name);
272 if (flag_limitStep) logi_vacuum->SetUserLimits(
new G4UserLimits(stepMax));
276 vacuum.
logi = logi_vacuum;
279 string phys_vacuum_name =
"phys_" + nameVac +
"_name";
281 new G4PVPlacement(0, G4ThreeVector(0, 0, 0), logi_vacuum, phys_vacuum_name, logi_polycone,
false, 0);
283 elements[nameVac] = vacuum;
288 std::vector<std::string> bendingSections;
290 for (
const auto& name : bendingSections) {
307 torus.
transform = G4Translate3D(torus_X0, 0.0, torus_Z0);
314 string geo_torusxx_name =
"geo_" + name +
"xx_name";
315 string geo_torusx_name =
"geo_" + name +
"x_name";
316 string geo_torus_name =
"geo_" + name +
"_name";
318 G4VSolid* geo_torus(NULL);
321 G4VSolid* geo_torusxx(NULL);
323 geo_torusxx =
new G4Torus(geo_torusxx_name, 0, torus_R, torus_RT, torus_SPHI, torus_DPHI);
325 geo_torusxx =
new G4Torus(geo_torusxx_name, torus_r, torus_R, torus_RT, torus_SPHI, torus_DPHI);
327 G4VSolid* geo_torusx =
new G4SubtractionSolid(geo_torusx_name, geo_torusxx, elements[subtract].geo,
328 torus.
transform.inverse()*elements[subtract].transform);
329 geo_torus =
new G4IntersectionSolid(geo_torus_name, geo_torusx, elements[
intersect].geo,
331 }
else if (subtract !=
"")
332 geo_torus =
new G4SubtractionSolid(geo_torus_name, geo_torusxx, elements[subtract].geo,
333 torus.
transform.inverse()*elements[subtract].transform);
335 geo_torus =
new G4IntersectionSolid(geo_torus_name, geo_torusxx, elements[
intersect].geo,
337 }
else if (type ==
"pipe")
338 geo_torus =
new G4Torus(geo_torus_name, 0, torus_R, torus_RT, torus_SPHI, torus_DPHI);
340 geo_torus =
new G4Torus(geo_torus_name, torus_r, torus_R, torus_RT, torus_SPHI, torus_DPHI);
342 torus.
geo = geo_torus;
347 string logi_torus_name =
"logi_" + name +
"_name";
348 G4LogicalVolume* logi_torus =
new G4LogicalVolume(torus.
geo, mat_torus, logi_torus_name);
352 torus.
logi = logi_torus;
355 string phys_torus_name =
"phys_" + name +
"_name";
356 new G4PVPlacement(torus.
transform, logi_torus, phys_torus_name, &topVolume,
false, 0);
358 elements[name] = torus;
360 if (type ==
"pipe" && torus_r != 0) {
364 string nameVac = name +
"Vac";
367 string geo_vacuumxx_name =
"geo_" + nameVac +
"xx_name";
368 string geo_vacuum_name =
"geo_" + nameVac +
"_name";
370 G4VSolid* geo_vacuumxx, *geo_vacuum;
372 geo_vacuumxx =
new G4Torus(geo_vacuumxx_name, 0.0, torus_r, torus_RT, torus_SPHI, torus_DPHI);
373 geo_vacuum =
new G4IntersectionSolid(geo_vacuum_name, geo_vacuumxx, geo_torus);
375 vacuum.
geo = geo_vacuum;
380 string logi_vacuum_name =
"logi_" + nameVac +
"_name";
381 G4LogicalVolume* logi_vacuum =
new G4LogicalVolume(vacuum.
geo, mat_vacuum, logi_vacuum_name);
382 if (flag_limitStep) logi_vacuum->SetUserLimits(
new G4UserLimits(stepMax));
386 vacuum.
logi = logi_vacuum;
389 string phys_vacuum_name =
"phys_" + nameVac +
"_name";
391 new G4PVPlacement(0, G4ThreeVector(0, 0, 0), logi_vacuum, phys_vacuum_name, logi_torus,
false, 0);
393 elements[nameVac] = vacuum;
401 std::vector<std::string> shields;
403 for (
const auto& name : shields) {
427 shield.
transform = G4Translate3D(shield_X0, shield_Y0, shield_Z0);
430 G4Transform3D transform_shield_hole = G4Translate3D(shield_hole_dX, shield_hole_dY, shield_hole_dZ);
433 string geo_shieldx_name =
"geo_" + name +
"x_name";
434 string geo_shield_hole_name =
"geo_" + name +
"_hole_name";
435 string geo_shield_name =
"geo_" + name +
"_name";
437 if (shield_hole_W == 0 || shield_hole_H == 0 || shield_hole_L == 0) {
438 G4Box* geo_shield =
new G4Box(geo_shield_name, shield_W / 2.0, shield_H / 2.0, shield_L / 2.0);
440 shield.
geo = geo_shield;
442 G4Box* geo_shieldx =
new G4Box(geo_shieldx_name, shield_W / 2.0, shield_H / 2.0, shield_L / 2.0);
443 G4Box* geo_shield_hole =
new G4Box(geo_shield_hole_name, shield_hole_W / 2.0, shield_hole_H / 2.0, shield_hole_L / 2.0);
444 G4SubtractionSolid* geo_shield =
new G4SubtractionSolid(geo_shield_name, geo_shieldx, geo_shield_hole,
445 transform_shield_hole);
447 shield.
geo = geo_shield;
453 string logi_shield_name =
"logi_" + name +
"_name";
454 G4LogicalVolume* logi_shield =
new G4LogicalVolume(shield.
geo, mat_shield, logi_shield_name);
456 shield.
logi = logi_shield;
461 string phys_shield_name =
"phys_" + name +
"_name";
462 new G4PVPlacement(shield.
transform, shield.
logi, phys_shield_name, &topVolume,
false, 0);
464 elements[name] = shield;
472 G4Tubs* geo_Tube =
new G4Tubs(
"geo_Tube_name", 3995 * CLHEP::mm, 4000 * CLHEP::mm, 29 * CLHEP::m, 0. * CLHEP::deg, 360.*CLHEP::deg);
474 G4LogicalVolume* logi_Tube =
new G4LogicalVolume(geo_Tube, mat_Tube,
"logi_Tube_name");
479 bool radiation_study =
false;
481 if (radiation_study && elements.count(
"GateShieldL")) {
482 new G4PVPlacement(elements[
"GateShieldL"].transform, logi_Tube,
"phys_Tube_name", &topVolume,
false, 0);
491 if (radiation_study) {
493 if (elements.count(
"PolyShieldL"))
495 if (elements.count(
"PolyShieldR"))
499 if (elements.count(
"ConcreteShieldL"))
501 if (elements.count(
"ConcreteShieldR"))
505 if (elements.count(
"GateShieldL"))
507 if (elements.count(
"GateShieldR"))
518 std::vector<std::string> collimators;
520 for (
const auto& name : collimators) {
528 string motherVolumeVacuum = motherVolume +
"Vac";
533 G4Translate3D translation;
534 elements[motherVolumeVacuum].transform.getDecomposition(scale, rotation, translation);
535 double zz = rotation.zz();
548 B2WARNING(
"Collimator " << name <<
" displacement d1 is set to " << collimator_d1 <<
"mm (must be negative)");
549 B2WARNING(
"Collimator " << name <<
" displacement d2 is set to " << collimator_d2 <<
"mm (must be positive)");
559 double head_dz = collimator_headH / 2.0;
560 if (type ==
"vertical") {
561 head_dx1 = collimator_th / 2.0;
562 head_dx2 = collimator_th / 2.0;
563 head_dy1 = ((collimator_maxW - collimator_minW) * collimator_headH / collimator_fullH + collimator_minW) / 2.0;
564 head_dy2 = collimator_minW / 2.0;
566 head_dx1 = ((collimator_maxW - collimator_minW) * collimator_headH / collimator_fullH + collimator_minW) / 2.0;
567 head_dx2 = collimator_minW / 2.0;
568 head_dy1 = collimator_th / 2.0;
569 head_dy2 = collimator_th / 2.0;
577 G4Transform3D transform_head1 = G4Translate3D(0.0, 0.0, collimator_Z);
578 G4Transform3D transform_head2 = G4Translate3D(0.0, 0.0, collimator_Z);
581 if (type ==
"vertical") {
582 transform_head1 = transform_head1 * G4Translate3D(0.0, -head_dz + collimator_d1, 0.0);
583 transform_head1 = transform_head1 * G4RotateX3D(-M_PI / 2 /
Unit::rad);
585 transform_head2 = transform_head2 * G4Translate3D(0.0, head_dz + collimator_d2, 0.0);
586 transform_head2 = transform_head2 * G4RotateX3D(M_PI / 2 /
Unit::rad);
589 transform_head1 = transform_head1 * G4Translate3D(-head_dz + collimator_d1, 0.0, 0.0);
590 transform_head1 = transform_head1 * G4RotateY3D(M_PI / 2 /
Unit::rad);
592 transform_head2 = transform_head2 * G4Translate3D(head_dz + collimator_d2, 0.0, 0.0);
593 transform_head2 = transform_head2 * G4RotateY3D(-M_PI / 2 /
Unit::rad);
595 transform_head1 = transform_head1 * G4Translate3D(head_dz - collimator_d1, 0.0, 0.0);
596 transform_head1 = transform_head1 * G4RotateY3D(-M_PI / 2 /
Unit::rad);
598 transform_head2 = transform_head2 * G4Translate3D(-head_dz - collimator_d2, 0.0, 0.0);
599 transform_head2 = transform_head2 * G4RotateY3D(M_PI / 2 /
Unit::rad);
603 collimator_head1.
transform = transform_head1;
604 collimator_head2.
transform = transform_head2;
607 string geo_headx_name =
"geo_" + name +
"_headx_name";
609 string geo_head1_name =
"geo_" + name +
"_head1_name";
610 string geo_head2_name =
"geo_" + name +
"_head2_name";
612 G4VSolid* geo_headx =
new G4Trd(geo_headx_name, head_dx1, head_dx2, head_dy1, head_dy2, head_dz);
614 G4VSolid* geo_head1 =
new G4IntersectionSolid(geo_head1_name, geo_headx, elements[motherVolumeVacuum].geo,
616 G4VSolid* geo_head2 =
new G4IntersectionSolid(geo_head2_name, geo_headx, elements[motherVolumeVacuum].geo,
619 collimator_head1.
geo = geo_head1;
620 collimator_head2.
geo = geo_head2;
625 string logi_head1_name =
"logi_" + name +
"_head1_name";
626 string logi_head2_name =
"logi_" + name +
"_head2_name";
627 G4LogicalVolume* logi_head1 =
new G4LogicalVolume(geo_head1, mat_head, logi_head1_name);
628 G4LogicalVolume* logi_head2 =
new G4LogicalVolume(geo_head2, mat_head, logi_head2_name);
635 double volume_head1 = logi_head1->GetSolid()->GetCubicVolume();
636 double volume_head2 = logi_head2->GetSolid()->GetCubicVolume();
638 collimator_head1.
logi = logi_head1;
639 collimator_head2.
logi = logi_head2;
642 string phys_head1_name =
"phys_" + name +
"_head1" +
"_name";
643 string phys_head2_name =
"phys_" + name +
"_head2" +
"_name";
644 if (volume_head1 != 0)
645 new G4PVPlacement(collimator_head1.
transform, logi_head1, phys_head1_name, elements[motherVolumeVacuum].logi,
false, 0);
646 if (volume_head2 != 0)
647 new G4PVPlacement(collimator_head2.
transform, logi_head2, phys_head2_name, elements[motherVolumeVacuum].logi,
false, 0);
650 collimator_head1.
transform = collimator_head1.
transform * elements[motherVolumeVacuum].transform;
651 collimator_head2.
transform = collimator_head2.
transform * elements[motherVolumeVacuum].transform;
653 string name_head1 = name +
"_head1";
654 string name_head2 = name +
"_head2";
655 elements[name_head1] = collimator_head1;
656 elements[name_head2] = collimator_head2;
666 double body_dz = (collimator_fullH - collimator_headH) / 2.0;
667 if (type ==
"vertical") {
668 body_dx1 = collimator_th / 2.0;
669 body_dx2 = collimator_th / 2.0;
670 body_dy1 = collimator_maxW / 2.0;
671 body_dy2 = ((collimator_maxW - collimator_minW) * collimator_headH / collimator_fullH + collimator_minW) / 2.0;
673 body_dx1 = collimator_maxW / 2.0;
674 body_dx2 = ((collimator_maxW - collimator_minW) * collimator_headH / collimator_fullH + collimator_minW) / 2.0;
675 body_dy1 = collimator_th / 2.0;
676 body_dy2 = collimator_th / 2.0;
684 if (type ==
"vertical") {
685 collimator_body1.
transform = G4Translate3D(0.0, -head_dz - body_dz, 0.0) * transform_head1;
686 collimator_body2.
transform = G4Translate3D(0.0, head_dz + body_dz, 0.0) * transform_head2;
689 collimator_body1.
transform = G4Translate3D(-head_dz - body_dz, 0.0, 0.0) * transform_head1;
690 collimator_body2.
transform = G4Translate3D(head_dz + body_dz, 0.0, 0.0) * transform_head2;
692 collimator_body1.
transform = G4Translate3D(head_dz + body_dz, 0.0, 0.0) * transform_head1;
693 collimator_body2.
transform = G4Translate3D(-head_dz - body_dz, 0.0, 0.0) * transform_head2;
698 string geo_bodyx_name =
"geo_" + name +
"_bodyx_name";
700 string geo_body1_name =
"geo_" + name +
"_body1_name";
701 string geo_body2_name =
"geo_" + name +
"_body2_name";
703 G4VSolid* geo_bodyx =
new G4Trd(geo_bodyx_name, body_dx1, body_dx2, body_dy1, body_dy2, body_dz);
705 G4VSolid* geo_body1 =
new G4IntersectionSolid(geo_body1_name, geo_bodyx, elements[motherVolumeVacuum].geo,
707 G4VSolid* geo_body2 =
new G4IntersectionSolid(geo_body2_name, geo_bodyx, elements[motherVolumeVacuum].geo,
710 collimator_body1.
geo = geo_body1;
711 collimator_body2.
geo = geo_body2;
716 string logi_body1_name =
"logi_" + name +
"_body1_name";
717 string logi_body2_name =
"logi_" + name +
"_body2_name";
718 G4LogicalVolume* logi_body1 =
new G4LogicalVolume(geo_body1, mat_body, logi_body1_name);
719 G4LogicalVolume* logi_body2 =
new G4LogicalVolume(geo_body2, mat_body, logi_body2_name);
726 double volume_body1 = logi_body1->GetSolid()->GetCubicVolume();
727 double volume_body2 = logi_body2->GetSolid()->GetCubicVolume();
729 collimator_body1.
logi = logi_body1;
730 collimator_body2.
logi = logi_body2;
733 string phys_body1_name =
"phys_" + name +
"_body1" +
"_name";
734 string phys_body2_name =
"phys_" + name +
"_body2" +
"_name";
735 if (volume_body1 != 0)
736 new G4PVPlacement(collimator_body1.
transform, logi_body1, phys_body1_name, elements[motherVolumeVacuum].logi,
false, 0);
737 if (volume_body2 != 0)
738 new G4PVPlacement(collimator_body2.
transform, logi_body2, phys_body2_name, elements[motherVolumeVacuum].logi,
false, 0);
741 collimator_body1.
transform = collimator_body1.
transform * elements[motherVolumeVacuum].transform;
742 collimator_body2.
transform = collimator_body2.
transform * elements[motherVolumeVacuum].transform;
744 string name_body1 = name +
"_body1";
745 string name_body2 = name +
"_body2";
746 elements[name_body1] = collimator_body1;
747 elements[name_body2] = collimator_body2;
The Class for BeamBackground Sensitive Detector.
const std::string & getParameterStr(const std::string &name) const
Get string parameter.
double getParameter(const std::string &name) const
Get parameter value.
static const double mm
[millimeters]
static const double rad
Standard of [angle].
static const double cm
Standard units with the value = 1.
static G4Material * get(const std::string &name)
Find given material.
FarBeamLineGeo m_config
geometry parameters object
GeoFarBeamLineCreator()
Constructor of the GeoFarBeamLineCreator class.
virtual ~GeoFarBeamLineCreator()
The destructor of the GeoFarBeamLineCreator class.
void createGeometry(G4LogicalVolume &topVolume, geometry::GeometryTypes type)
Create detector geometry.
SensitiveDetector * m_sensitive
Sensitive detector.
The IR Sensitive Detector class.
int intersect(const TRGCDCLpar &lp1, const TRGCDCLpar &lp2, CLHEP::HepVector &v1, CLHEP::HepVector &v2)
intersection
void setVisibility(G4LogicalVolume &volume, bool visible)
Helper function to quickly set the visibility of a given volume.
void setColor(G4LogicalVolume &volume, const std::string &color)
Set the color of a logical volume.
GeometryTypes
Flag indiciating the type of geometry to be used.
Abstract base class for different kinds of events.
The struct for FarBeamLineElement.
G4LogicalVolume * logi
Logical volume.
G4VSolid * geo
Solid volume.
G4Transform3D transform
Transformation.