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/format.hpp>
20#include <boost/foreach.hpp>
21#include <boost/algorithm/string.hpp>
23#include <G4LogicalVolume.hh>
24#include <G4PVPlacement.hh>
30#include <G4Polycone.hh>
32#include <G4IntersectionSolid.hh>
33#include <G4SubtractionSolid.hh>
34#include <G4UserLimits.hh>
46 using namespace geometry;
54 geometry::CreatorFactory<GeoFarBeamLineCreator> GeoFarBeamLineFactory(
"FarBeamLineCreator");
82 map<string, FarBeamLineElement> elements;
94 std::string prep =
"TubeR.";
98 std::vector<double> TubeR_Z(TubeR_N);
99 std::vector<double> TubeR_R(TubeR_N);
100 std::vector<double> TubeR_r(TubeR_N);
102 for (
int i = 0; i < TubeR_N; ++i) {
103 ostringstream ossZID;
106 ostringstream ossRID;
109 ostringstream ossrID;
117 tubeR.
transform = G4Translate3D(0.0, 0.0, 0.0);
120 tubeR.
geo =
new G4Polycone(
"geo_TubeR_name", 0.0, 2 * M_PI, TubeR_N, &(TubeR_Z[0]), &(TubeR_r[0]), &(TubeR_R[0]));
124 elements[
"TubeR"] = tubeR;
136 std::vector<double> TubeL_Z(TubeL_N);
137 std::vector<double> TubeL_R(TubeL_N);
138 std::vector<double> TubeL_r(TubeL_N);
140 for (
int i = 0; i < TubeL_N; ++i) {
141 ostringstream ossZID;
144 ostringstream ossRID;
147 ostringstream ossrID;
155 tubeL.
transform = G4Translate3D(0.0, 0.0, 0.0);
158 tubeL.
geo =
new G4Polycone(
"geo_TubeL_name", 0.0, 2 * M_PI, TubeL_N, &(TubeL_Z[0]), &(TubeL_r[0]), &(TubeL_R[0]));
162 elements[
"TubeL"] = tubeL;
164 std::vector<double> zero_r(N, 0.);
166 std::vector<std::string> straightSections;
168 for (
const auto& name : straightSections) {
178 std::vector<double> Polycone_Z(N);
179 std::vector<double> Polycone_R(N);
180 std::vector<double> Polycone_r(N);
194 polycone.
transform = G4Translate3D(Polycone_X0, Polycone_Y0, Polycone_Z0);
196 if (Polycone_PHIYZ != 0)
203 string geo_polyconexx_name =
"geo_" + name +
"xx_name";
204 string geo_polyconex_name =
"geo_" + name +
"x_name";
205 string geo_polycone_name =
"geo_" + name +
"_name";
207 G4VSolid* geo_polyconexx(NULL), *geo_polycone(NULL);
211 geo_polyconexx =
new G4Polycone(geo_polyconexx_name, 0.0, 2 * M_PI, N, &(Polycone_Z[0]), &(zero_r[0]), &(Polycone_R[0]));
213 geo_polyconexx =
new G4Polycone(geo_polyconexx_name, 0.0, 2 * M_PI, N, &(Polycone_Z[0]), &(Polycone_r[0]), &(Polycone_R[0]));
214 else if (type ==
"pipe")
215 geo_polycone =
new G4Polycone(geo_polycone_name, 0.0, 2 * M_PI, N, &(Polycone_Z[0]), &(zero_r[0]), &(Polycone_R[0]));
217 geo_polycone =
new G4Polycone(geo_polycone_name, 0.0, 2 * M_PI, N, &(Polycone_Z[0]), &(Polycone_r[0]), &(Polycone_R[0]));
221 G4VSolid* geo_polyconex =
new G4SubtractionSolid(geo_polyconex_name, geo_polyconexx, elements[subtract].geo,
222 polycone.
transform.inverse()*elements[subtract].transform);
223 geo_polycone =
new G4IntersectionSolid(geo_polycone_name, geo_polyconex, elements[
intersect].geo,
225 }
else if (subtract !=
"")
226 geo_polycone =
new G4SubtractionSolid(geo_polycone_name, geo_polyconexx, elements[subtract].geo,
227 polycone.
transform.inverse()*elements[subtract].transform);
229 geo_polycone =
new G4IntersectionSolid(geo_polycone_name, geo_polyconexx, elements[
intersect].geo,
232 polycone.
geo = geo_polycone;
237 string logi_polycone_name =
"logi_" + name +
"_name";
238 G4LogicalVolume* logi_polycone =
new G4LogicalVolume(polycone.
geo, mat_polycone, logi_polycone_name);
239 setColor(*logi_polycone,
"#CC0000");
242 polycone.
logi = logi_polycone;
245 string phys_polycone_name =
"phys_" + name +
"_name";
246 new G4PVPlacement(polycone.
transform, logi_polycone, phys_polycone_name, &topVolume,
false, 0);
248 elements[name] = polycone;
251 for (
int i = 0; i < N; ++i)
252 sum += Polycone_r[i];
253 if (type ==
"pipe" && sum != 0) {
257 string nameVac = name +
"Vac";
260 string geo_vacuumxx_name =
"geo_" + nameVac +
"xx_name";
261 string geo_vacuum_name =
"geo_" + nameVac +
"_name";
263 G4VSolid* geo_vacuumxx, *geo_vacuum;
265 geo_vacuumxx =
new G4Polycone(geo_vacuumxx_name, 0.0, 2 * M_PI, N, &(Polycone_Z[0]), &(zero_r[0]), &(Polycone_r[0]));
266 geo_vacuum =
new G4IntersectionSolid(geo_vacuumxx_name, geo_vacuumxx, geo_polycone);
268 vacuum.
geo = geo_vacuum;
273 string logi_vacuum_name =
"logi_" + nameVac +
"_name";
274 G4LogicalVolume* logi_vacuum =
new G4LogicalVolume(vacuum.
geo, mat_vacuum, logi_vacuum_name);
275 if (flag_limitStep) logi_vacuum->SetUserLimits(
new G4UserLimits(stepMax));
279 vacuum.
logi = logi_vacuum;
282 string phys_vacuum_name =
"phys_" + nameVac +
"_name";
284 new G4PVPlacement(0, G4ThreeVector(0, 0, 0), logi_vacuum, phys_vacuum_name, logi_polycone,
false, 0);
286 elements[nameVac] = vacuum;
291 std::vector<std::string> bendingSections;
293 for (
const auto& name : bendingSections) {
310 torus.
transform = G4Translate3D(torus_X0, 0.0, torus_Z0);
317 string geo_torusxx_name =
"geo_" + name +
"xx_name";
318 string geo_torusx_name =
"geo_" + name +
"x_name";
319 string geo_torus_name =
"geo_" + name +
"_name";
321 G4VSolid* geo_torus(NULL);
324 G4VSolid* geo_torusxx(NULL);
326 geo_torusxx =
new G4Torus(geo_torusxx_name, 0, torus_R, torus_RT, torus_SPHI, torus_DPHI);
328 geo_torusxx =
new G4Torus(geo_torusxx_name, torus_r, torus_R, torus_RT, torus_SPHI, torus_DPHI);
330 G4VSolid* geo_torusx =
new G4SubtractionSolid(geo_torusx_name, geo_torusxx, elements[subtract].geo,
331 torus.
transform.inverse()*elements[subtract].transform);
332 geo_torus =
new G4IntersectionSolid(geo_torus_name, geo_torusx, elements[
intersect].geo,
334 }
else if (subtract !=
"")
335 geo_torus =
new G4SubtractionSolid(geo_torus_name, geo_torusxx, elements[subtract].geo,
336 torus.
transform.inverse()*elements[subtract].transform);
338 geo_torus =
new G4IntersectionSolid(geo_torus_name, geo_torusxx, elements[
intersect].geo,
340 }
else if (type ==
"pipe")
341 geo_torus =
new G4Torus(geo_torus_name, 0, torus_R, torus_RT, torus_SPHI, torus_DPHI);
343 geo_torus =
new G4Torus(geo_torus_name, torus_r, torus_R, torus_RT, torus_SPHI, torus_DPHI);
345 torus.
geo = geo_torus;
350 string logi_torus_name =
"logi_" + name +
"_name";
351 G4LogicalVolume* logi_torus =
new G4LogicalVolume(torus.
geo, mat_torus, logi_torus_name);
355 torus.
logi = logi_torus;
358 string phys_torus_name =
"phys_" + name +
"_name";
359 new G4PVPlacement(torus.
transform, logi_torus, phys_torus_name, &topVolume,
false, 0);
361 elements[name] = torus;
363 if (type ==
"pipe" && torus_r != 0) {
367 string nameVac = name +
"Vac";
370 string geo_vacuumxx_name =
"geo_" + nameVac +
"xx_name";
371 string geo_vacuum_name =
"geo_" + nameVac +
"_name";
373 G4VSolid* geo_vacuumxx, *geo_vacuum;
375 geo_vacuumxx =
new G4Torus(geo_vacuumxx_name, 0.0, torus_r, torus_RT, torus_SPHI, torus_DPHI);
376 geo_vacuum =
new G4IntersectionSolid(geo_vacuum_name, geo_vacuumxx, geo_torus);
378 vacuum.
geo = geo_vacuum;
383 string logi_vacuum_name =
"logi_" + nameVac +
"_name";
384 G4LogicalVolume* logi_vacuum =
new G4LogicalVolume(vacuum.
geo, mat_vacuum, logi_vacuum_name);
385 if (flag_limitStep) logi_vacuum->SetUserLimits(
new G4UserLimits(stepMax));
389 vacuum.
logi = logi_vacuum;
392 string phys_vacuum_name =
"phys_" + nameVac +
"_name";
394 new G4PVPlacement(0, G4ThreeVector(0, 0, 0), logi_vacuum, phys_vacuum_name, logi_torus,
false, 0);
396 elements[nameVac] = vacuum;
404 std::vector<std::string> shields;
406 for (
const auto& name : shields) {
430 shield.
transform = G4Translate3D(shield_X0, shield_Y0, shield_Z0);
433 G4Transform3D transform_shield_hole = G4Translate3D(shield_hole_dX, shield_hole_dY, shield_hole_dZ);
436 string geo_shieldx_name =
"geo_" + name +
"x_name";
437 string geo_shield_hole_name =
"geo_" + name +
"_hole_name";
438 string geo_shield_name =
"geo_" + name +
"_name";
440 if (shield_hole_W == 0 || shield_hole_H == 0 || shield_hole_L == 0) {
441 G4Box* geo_shield =
new G4Box(geo_shield_name, shield_W / 2.0, shield_H / 2.0, shield_L / 2.0);
443 shield.
geo = geo_shield;
445 G4Box* geo_shieldx =
new G4Box(geo_shieldx_name, shield_W / 2.0, shield_H / 2.0, shield_L / 2.0);
446 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);
447 G4SubtractionSolid* geo_shield =
new G4SubtractionSolid(geo_shield_name, geo_shieldx, geo_shield_hole,
448 transform_shield_hole);
450 shield.
geo = geo_shield;
456 string logi_shield_name =
"logi_" + name +
"_name";
457 G4LogicalVolume* logi_shield =
new G4LogicalVolume(shield.
geo, mat_shield, logi_shield_name);
459 shield.
logi = logi_shield;
464 string phys_shield_name =
"phys_" + name +
"_name";
465 new G4PVPlacement(shield.
transform, shield.
logi, phys_shield_name, &topVolume,
false, 0);
467 elements[name] = shield;
475 G4Tubs* geo_Tube =
new G4Tubs(
"geo_Tube_name", 3995 * CLHEP::mm, 4000 * CLHEP::mm, 29 * CLHEP::m, 0. * CLHEP::deg, 360.*CLHEP::deg);
477 G4LogicalVolume* logi_Tube =
new G4LogicalVolume(geo_Tube, mat_Tube,
"logi_Tube_name");
482 bool radiation_study =
false;
484 if (radiation_study && elements.count(
"GateShieldL")) {
485 new G4PVPlacement(elements[
"GateShieldL"].transform, logi_Tube,
"phys_Tube_name", &topVolume,
false, 0);
494 if (radiation_study) {
496 if (elements.count(
"PolyShieldL"))
498 if (elements.count(
"PolyShieldR"))
502 if (elements.count(
"ConcreteShieldL"))
504 if (elements.count(
"ConcreteShieldR"))
508 if (elements.count(
"GateShieldL"))
510 if (elements.count(
"GateShieldR"))
521 std::vector<std::string> collimators;
523 for (
const auto& name : collimators) {
531 string motherVolumeVacuum = motherVolume +
"Vac";
536 G4Translate3D translation;
537 elements[motherVolumeVacuum].transform.getDecomposition(scale, rotation, translation);
538 double zz = rotation.zz();
551 B2WARNING(
"Collimator " << name <<
" displacement d1 is set to " << collimator_d1 <<
"mm (must be negative)");
552 B2WARNING(
"Collimator " << name <<
" displacement d2 is set to " << collimator_d2 <<
"mm (must be positive)");
562 double head_dz = collimator_headH / 2.0;
563 if (type ==
"vertical") {
564 head_dx1 = collimator_th / 2.0;
565 head_dx2 = collimator_th / 2.0;
566 head_dy1 = ((collimator_maxW - collimator_minW) * collimator_headH / collimator_fullH + collimator_minW) / 2.0;
567 head_dy2 = collimator_minW / 2.0;
569 head_dx1 = ((collimator_maxW - collimator_minW) * collimator_headH / collimator_fullH + collimator_minW) / 2.0;
570 head_dx2 = collimator_minW / 2.0;
571 head_dy1 = collimator_th / 2.0;
572 head_dy2 = collimator_th / 2.0;
580 G4Transform3D transform_head1 = G4Translate3D(0.0, 0.0, collimator_Z);
581 G4Transform3D transform_head2 = G4Translate3D(0.0, 0.0, collimator_Z);
584 if (type ==
"vertical") {
585 transform_head1 = transform_head1 * G4Translate3D(0.0, -head_dz + collimator_d1, 0.0);
586 transform_head1 = transform_head1 * G4RotateX3D(-M_PI / 2 /
Unit::rad);
588 transform_head2 = transform_head2 * G4Translate3D(0.0, head_dz + collimator_d2, 0.0);
589 transform_head2 = transform_head2 * G4RotateX3D(M_PI / 2 /
Unit::rad);
592 transform_head1 = transform_head1 * G4Translate3D(-head_dz + collimator_d1, 0.0, 0.0);
593 transform_head1 = transform_head1 * G4RotateY3D(M_PI / 2 /
Unit::rad);
595 transform_head2 = transform_head2 * G4Translate3D(head_dz + collimator_d2, 0.0, 0.0);
596 transform_head2 = transform_head2 * G4RotateY3D(-M_PI / 2 /
Unit::rad);
598 transform_head1 = transform_head1 * G4Translate3D(head_dz - collimator_d1, 0.0, 0.0);
599 transform_head1 = transform_head1 * G4RotateY3D(-M_PI / 2 /
Unit::rad);
601 transform_head2 = transform_head2 * G4Translate3D(-head_dz - collimator_d2, 0.0, 0.0);
602 transform_head2 = transform_head2 * G4RotateY3D(M_PI / 2 /
Unit::rad);
606 collimator_head1.
transform = transform_head1;
607 collimator_head2.
transform = transform_head2;
610 string geo_headx_name =
"geo_" + name +
"_headx_name";
612 string geo_head1_name =
"geo_" + name +
"_head1_name";
613 string geo_head2_name =
"geo_" + name +
"_head2_name";
615 G4VSolid* geo_headx =
new G4Trd(geo_headx_name, head_dx1, head_dx2, head_dy1, head_dy2, head_dz);
617 G4VSolid* geo_head1 =
new G4IntersectionSolid(geo_head1_name, geo_headx, elements[motherVolumeVacuum].geo,
619 G4VSolid* geo_head2 =
new G4IntersectionSolid(geo_head2_name, geo_headx, elements[motherVolumeVacuum].geo,
622 collimator_head1.
geo = geo_head1;
623 collimator_head2.
geo = geo_head2;
628 string logi_head1_name =
"logi_" + name +
"_head1_name";
629 string logi_head2_name =
"logi_" + name +
"_head2_name";
630 G4LogicalVolume* logi_head1 =
new G4LogicalVolume(geo_head1, mat_head, logi_head1_name);
631 G4LogicalVolume* logi_head2 =
new G4LogicalVolume(geo_head2, mat_head, logi_head2_name);
638 double volume_head1 = logi_head1->GetSolid()->GetCubicVolume();
639 double volume_head2 = logi_head2->GetSolid()->GetCubicVolume();
641 collimator_head1.
logi = logi_head1;
642 collimator_head2.
logi = logi_head2;
645 string phys_head1_name =
"phys_" + name +
"_head1" +
"_name";
646 string phys_head2_name =
"phys_" + name +
"_head2" +
"_name";
647 if (volume_head1 != 0)
648 new G4PVPlacement(collimator_head1.
transform, logi_head1, phys_head1_name, elements[motherVolumeVacuum].logi,
false, 0);
649 if (volume_head2 != 0)
650 new G4PVPlacement(collimator_head2.
transform, logi_head2, phys_head2_name, elements[motherVolumeVacuum].logi,
false, 0);
653 collimator_head1.
transform = collimator_head1.
transform * elements[motherVolumeVacuum].transform;
654 collimator_head2.
transform = collimator_head2.
transform * elements[motherVolumeVacuum].transform;
656 string name_head1 = name +
"_head1";
657 string name_head2 = name +
"_head2";
658 elements[name_head1] = collimator_head1;
659 elements[name_head2] = collimator_head2;
669 double body_dz = (collimator_fullH - collimator_headH) / 2.0;
670 if (type ==
"vertical") {
671 body_dx1 = collimator_th / 2.0;
672 body_dx2 = collimator_th / 2.0;
673 body_dy1 = collimator_maxW / 2.0;
674 body_dy2 = ((collimator_maxW - collimator_minW) * collimator_headH / collimator_fullH + collimator_minW) / 2.0;
676 body_dx1 = collimator_maxW / 2.0;
677 body_dx2 = ((collimator_maxW - collimator_minW) * collimator_headH / collimator_fullH + collimator_minW) / 2.0;
678 body_dy1 = collimator_th / 2.0;
679 body_dy2 = collimator_th / 2.0;
687 if (type ==
"vertical") {
688 collimator_body1.
transform = G4Translate3D(0.0, -head_dz - body_dz, 0.0) * transform_head1;
689 collimator_body2.
transform = G4Translate3D(0.0, head_dz + body_dz, 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;
695 collimator_body1.
transform = G4Translate3D(head_dz + body_dz, 0.0, 0.0) * transform_head1;
696 collimator_body2.
transform = G4Translate3D(-head_dz - body_dz, 0.0, 0.0) * transform_head2;
701 string geo_bodyx_name =
"geo_" + name +
"_bodyx_name";
703 string geo_body1_name =
"geo_" + name +
"_body1_name";
704 string geo_body2_name =
"geo_" + name +
"_body2_name";
706 G4VSolid* geo_bodyx =
new G4Trd(geo_bodyx_name, body_dx1, body_dx2, body_dy1, body_dy2, body_dz);
708 G4VSolid* geo_body1 =
new G4IntersectionSolid(geo_body1_name, geo_bodyx, elements[motherVolumeVacuum].geo,
710 G4VSolid* geo_body2 =
new G4IntersectionSolid(geo_body2_name, geo_bodyx, elements[motherVolumeVacuum].geo,
713 collimator_body1.
geo = geo_body1;
714 collimator_body2.
geo = geo_body2;
719 string logi_body1_name =
"logi_" + name +
"_body1_name";
720 string logi_body2_name =
"logi_" + name +
"_body2_name";
721 G4LogicalVolume* logi_body1 =
new G4LogicalVolume(geo_body1, mat_body, logi_body1_name);
722 G4LogicalVolume* logi_body2 =
new G4LogicalVolume(geo_body2, mat_body, logi_body2_name);
729 double volume_body1 = logi_body1->GetSolid()->GetCubicVolume();
730 double volume_body2 = logi_body2->GetSolid()->GetCubicVolume();
732 collimator_body1.
logi = logi_body1;
733 collimator_body2.
logi = logi_body2;
736 string phys_body1_name =
"phys_" + name +
"_body1" +
"_name";
737 string phys_body2_name =
"phys_" + name +
"_body2" +
"_name";
738 if (volume_body1 != 0)
739 new G4PVPlacement(collimator_body1.
transform, logi_body1, phys_body1_name, elements[motherVolumeVacuum].logi,
false, 0);
740 if (volume_body2 != 0)
741 new G4PVPlacement(collimator_body2.
transform, logi_body2, phys_body2_name, elements[motherVolumeVacuum].logi,
false, 0);
744 collimator_body1.
transform = collimator_body1.
transform * elements[motherVolumeVacuum].transform;
745 collimator_body2.
transform = collimator_body2.
transform * elements[motherVolumeVacuum].transform;
747 string name_body1 = name +
"_body1";
748 string name_body2 = name +
"_body2";
749 elements[name_body1] = collimator_body1;
750 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.