From 07e741c9fdb11a22ee91aa24dea35c4840c3ee62 Mon Sep 17 00:00:00 2001 From: Runze Ren Date: Tue, 25 Apr 2023 13:57:03 -0700 Subject: [PATCH 01/12] Changing to 6 TOP layers --- src/geometry/box/Box.cc | 48 +++++++++++++++++++++++++---------------- 1 file changed, 30 insertions(+), 18 deletions(-) diff --git a/src/geometry/box/Box.cc b/src/geometry/box/Box.cc index 33c483c..7c6ffb4 100644 --- a/src/geometry/box/Box.cc +++ b/src/geometry/box/Box.cc @@ -100,10 +100,11 @@ constexpr auto steel_height = 0.03*m; constexpr auto air_gap = 30*m; -constexpr auto scintillator_casing_thickness = 0.005*m; +constexpr auto scintillator_casing_thickness = 0.003*m; constexpr auto layer_spacing = 1.0*m; -constexpr auto layer_count = 7UL; +constexpr auto layer_spacing_top = 0.8*m; // Top layers have different spacing +constexpr auto layer_count = 8UL; constexpr auto module_x_edge_length = 9.0*m; constexpr auto module_y_edge_length = 9.0*m; @@ -115,13 +116,13 @@ constexpr auto x_edge_increase = 2*full_layer_height + 4*wall_gap; constexpr auto layer_w_case = full_layer_height; -constexpr auto full_module_height = (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 5.0*layer_w_case + 4.0*layer_spacing; +constexpr auto full_module_height = (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 6.0*layer_w_case + 5.0*layer_spacing_top; constexpr auto scintillator_z_position = 0.00; constexpr auto wall_height = 20*m; -constexpr int NBEAMLAYERS = 7; +constexpr int NBEAMLAYERS = 8; constexpr auto beam_x_edge_length = 0.10*m; constexpr auto beam_y_edge_length = 0.10*m; constexpr auto beam_thickness = 0.02*m; @@ -129,21 +130,24 @@ constexpr auto beam_thickness = 0.02*m; constexpr auto full_detector_height = full_module_height + steel_height + 3.0*layer_w_case + 2.0*layer_spacing; constexpr auto half_detector_height = 0.5L * full_detector_height; -constexpr double layer_z_displacement[7] = {-0.5*full_module_height + (20.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 0.5*layer_w_case, +constexpr double layer_z_displacement[8] = {-0.5*full_module_height + (20.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 0.5*layer_w_case, -0.5*full_module_height + (20.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + layer_spacing + 1.5*layer_w_case, -0.5*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 0.5*layer_w_case, - -0.5*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + layer_spacing + 1.5*layer_w_case, - -0.5*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 2*layer_spacing + 2.5*layer_w_case, - -0.5*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 3*layer_spacing + 3.5*layer_w_case, - -0.5*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 4*layer_spacing + 4.5*layer_w_case}; + -0.5*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + layer_spacing_top + 1.5*layer_w_case, + -0.5*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 2*layer_spacing_top + 2.5*layer_w_case, + -0.5*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 3*layer_spacing_top + 3.5*layer_w_case, + -0.5*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 4*layer_spacing_top + 4.5*layer_w_case, + -0.5*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 5*layer_spacing_top + 5.5*layer_w_case}; -constexpr double module_beam_heights[7] = {20.0*m - 3*layer_w_case - 2*layer_spacing, + +constexpr double module_beam_heights[8] = {20.0*m - 3*layer_w_case - 2*layer_spacing, layer_spacing, 5.0*m - 2*layer_w_case - layer_spacing, - layer_spacing, - layer_spacing, - layer_spacing, - layer_spacing}; + layer_spacing_top, + layer_spacing_top, + layer_spacing_top, + layer_spacing_top, + layer_spacing_top}; // constexpr double module_beam_z_pos[9] = {-0.50*full_module_height + 0.50*module_beam_heights[0] + layer_w_case, // -0.50*full_module_height + 2*layer_w_case + layer_spacing + 0.50*module_beam_heights[1], @@ -155,13 +159,14 @@ constexpr double module_beam_heights[7] = {20.0*m - 3*layer_w_case - 2*layer_spa // -0.50*full_module_height + 25.0*m + 3*layer_w_case + 2*layer_spacing + 0.50*module_beam_heights[7], // -0.50*full_module_height + 25.0*m + 4*layer_w_case + 3*layer_spacing + 0.50*module_beam_heights[8]}; -constexpr double module_beam_z_pos[7] = {-0.50*full_module_height + 0.50*module_beam_heights[0], +constexpr double module_beam_z_pos[8] = {-0.50*full_module_height + 0.50*module_beam_heights[0], -0.50*full_module_height + (20.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + layer_w_case + 0.50*module_beam_heights[1], -0.50*full_module_height + (20.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 2*layer_w_case + layer_spacing + 0.50*module_beam_heights[2], -0.50*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + layer_w_case + 0.50*module_beam_heights[3], - -0.50*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 2*layer_w_case + layer_spacing + 0.50*module_beam_heights[4], - -0.50*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 3*layer_w_case + 2*layer_spacing + 0.50*module_beam_heights[5], - -0.50*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 4*layer_w_case + 3*layer_spacing + 0.50*module_beam_heights[6]}; + -0.50*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 2*layer_w_case + 1*layer_spacing_top + 0.50*module_beam_heights[4], + -0.50*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 3*layer_w_case + 2*layer_spacing_top + 0.50*module_beam_heights[5], + -0.50*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 4*layer_w_case + 3*layer_spacing_top + 0.50*module_beam_heights[6], + -0.50*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 5*layer_w_case + 3*layer_spacing_top + 0.50*module_beam_heights[7]}; @@ -211,6 +216,13 @@ Detector::Detector() : G4VSensitiveDetector("MATHUSLA/MU/Box") { collectionName.insert("Box_HC"); for (auto& scintillator : _scintillators) scintillator->Register(this); + + // Print out the detector dislacement in z direction + std::cout<<"Layer z displacement:\n \n"< Date: Thu, 27 Apr 2023 13:30:26 +0200 Subject: [PATCH 02/12] Update analysis.hh --- include/analysis.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/analysis.hh b/include/analysis.hh index cd77c8a..a31d2cf 100644 --- a/include/analysis.hh +++ b/include/analysis.hh @@ -23,7 +23,7 @@ #include #include -#include +#include namespace MATHUSLA { namespace MU { From 4d8ee31cc771cfc56ba1d71b4c90647a272b3d6d Mon Sep 17 00:00:00 2001 From: Abdulrahman Mohamed <88555369+Abdulrahma-M@users.noreply.github.com> Date: Thu, 27 Apr 2023 13:34:37 +0200 Subject: [PATCH 03/12] Update Box.cc --- src/geometry/box/Box.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/geometry/box/Box.cc b/src/geometry/box/Box.cc index 33c483c..9dc2ebd 100644 --- a/src/geometry/box/Box.cc +++ b/src/geometry/box/Box.cc @@ -558,9 +558,9 @@ G4VPhysicalVolume* Detector::Construct(G4LogicalVolume* world) { // pre_data->Branch("Y_H", &Y_POS_HIT, "Y_H/D"); auto DetectorVolume = Construction::BoxVolume("Box", x_edge_length + x_edge_increase, y_edge_length, full_detector_height, - Construction::Material::Air, G4VisAttributes::Invisible); + Construction::Material::Air, G4VisAttributes::GetInvisible()); - //DetectorVolume->SetVisAttributes(G4VisAttributes::Invisible); + //DetectorVolume->SetVisAttributes(G4VisAttributes::GetInvisible()); for (int module_number = 0; module_number < NMODULES; module_number++){ auto current = Detector::ConstructModule(DetectorVolume, module_number, @@ -761,7 +761,7 @@ G4VPhysicalVolume* Detector::ConstructEarth(G4LogicalVolume* world){ auto CMS_Detector_logical = CMSRingVolume(); auto UXC_55_air_v1 = new G4SubtractionSolid("UXC_55_air_v1", UXC_55_cavern_solid, UXC55_outer_solid); auto UXC_55_air_v2 = new G4SubtractionSolid("UXC_55_air_v2", UXC_55_air_v1, CMS_Detector_logical->GetSolid()); - auto UXC55_air_logical = Volume("UXC55_air", UXC_55_air_v2, Construction::Material::Air, G4VisAttributes::Invisible); + auto UXC55_air_logical = Volume("UXC55_air", UXC_55_air_v2, Construction::Material::Air, G4VisAttributes::GetInvisible()); Construction::PlaceVolume(UXC55_outer_logical, earth, Cavern_Transform()*Construction::Rotate(0, 1, 0, 90*deg) ); Construction::PlaceVolume(CMS_Detector_logical, earth, Cavern_Transform()*Construction::Rotate(0, 1, 0, 90*deg) ); @@ -785,7 +785,7 @@ G4VPhysicalVolume* Detector::ConstructEarth(G4LogicalVolume* world){ AS_Depth - 2* AS_Thickness, AS_Height - 2* AS_Thickness, Construction::Material::Air, - G4VisAttributes::Invisible); + G4VisAttributes::GetInvisible()); Construction::PlaceVolume(Access_Shaft_outer_logical, earth, Access_Shaft_Transform() ); Construction::PlaceVolume(Access_Shaft_Air, earth, Access_Shaft_Transform()); From f50f9f262eda0f1873d7a41f2aedf4ab11aa0071 Mon Sep 17 00:00:00 2001 From: Abdulrahman Mohamed <88555369+Abdulrahma-M@users.noreply.github.com> Date: Thu, 27 Apr 2023 13:35:22 +0200 Subject: [PATCH 04/12] Update Cosmic.cc --- src/geometry/cosmic/Cosmic.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/geometry/cosmic/Cosmic.cc b/src/geometry/cosmic/Cosmic.cc index 904b34c..c3c6f4b 100644 --- a/src/geometry/cosmic/Cosmic.cc +++ b/src/geometry/cosmic/Cosmic.cc @@ -559,9 +559,9 @@ G4VPhysicalVolume* Detector::Construct(G4LogicalVolume* world) { // pre_data->Branch("Y_H", &Y_POS_HIT, "Y_H/D"); auto DetectorVolume = Construction::BoxVolume("Box", x_edge_length + x_edge_increase, y_edge_length, full_detector_height, - Construction::Material::Air, G4VisAttributes::Invisible); + Construction::Material::Air, G4VisAttributes::GetInvisible()); - //DetectorVolume->SetVisAttributes(G4VisAttributes::Invisible); + //DetectorVolume->SetVisAttributes(G4VisAttributes::GetInvisible(); for (int module_number = 0; module_number < NMODULES; module_number++){ auto current = Detector::ConstructModule(DetectorVolume, module_number, @@ -772,7 +772,7 @@ G4VPhysicalVolume* Detector::ConstructEarth(G4LogicalVolume* world){ // auto CMS_Detector_logical = CMSRingVolume(); // auto UXC_55_air_v1 = new G4SubtractionSolid("UXC_55_air_v1", UXC_55_cavern_solid, UXC55_outer_solid); // auto UXC_55_air_v2 = new G4SubtractionSolid("UXC_55_air_v2", UXC_55_air_v1, CMS_Detector_logical->GetSolid()); - // auto UXC55_air_logical = Volume("UXC55_air", UXC_55_air_v2, Construction::Material::Air, G4VisAttributes::Invisible); + // auto UXC55_air_logical = Volume("UXC55_air", UXC_55_air_v2, Construction::Material::Air, G4VisAttributes::GetInvisible()); // Construction::PlaceVolume(UXC55_outer_logical, earth, Cavern_Transform()*Construction::Rotate(0, 1, 0, 90*deg) ); // Construction::PlaceVolume(CMS_Detector_logical, earth, Cavern_Transform()*Construction::Rotate(0, 1, 0, 90*deg) ); @@ -796,7 +796,7 @@ G4VPhysicalVolume* Detector::ConstructEarth(G4LogicalVolume* world){ // AS_Depth - 2* AS_Thickness, // AS_Height - 2* AS_Thickness, // Construction::Material::Air, - // G4VisAttributes::Invisible); + // G4VisAttributes::GetInvisible()); // Construction::PlaceVolume(Access_Shaft_outer_logical, earth, Access_Shaft_Transform() ); // Construction::PlaceVolume(Access_Shaft_Air, earth, Access_Shaft_Transform()); From efac13140d6cabb28c45602948e8e7c4397b010e Mon Sep 17 00:00:00 2001 From: Runze Ren Date: Tue, 2 May 2023 06:45:19 -0700 Subject: [PATCH 05/12] Changing geometry, 2nd try --- src/geometry/box/Box.cc | 157 ++++++++++++++++++++++------------------ 1 file changed, 87 insertions(+), 70 deletions(-) diff --git a/src/geometry/box/Box.cc b/src/geometry/box/Box.cc index 7c6ffb4..2ae00e4 100644 --- a/src/geometry/box/Box.cc +++ b/src/geometry/box/Box.cc @@ -82,7 +82,7 @@ G4ThreadLocal Tracking::HitCollection* _hit_collection; constexpr int scintillators_per_layer{400}; constexpr int NMODULES{100}; -constexpr int n_top_layers{5}; +constexpr int n_top_layers{6}; constexpr auto x_edge_length = 99.0*m; constexpr auto y_edge_length = 99.0*m; constexpr auto x_displacement = 70.0*m; @@ -116,7 +116,7 @@ constexpr auto x_edge_increase = 2*full_layer_height + 4*wall_gap; constexpr auto layer_w_case = full_layer_height; -constexpr auto full_module_height = (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 6.0*layer_w_case + 5.0*layer_spacing_top; +constexpr auto full_module_height = (25.0*m) + 6.0*layer_w_case + 5.0*layer_spacing_top; constexpr auto scintillator_z_position = 0.00; @@ -127,22 +127,35 @@ constexpr auto beam_x_edge_length = 0.10*m; constexpr auto beam_y_edge_length = 0.10*m; constexpr auto beam_thickness = 0.02*m; -constexpr auto full_detector_height = full_module_height + steel_height + 3.0*layer_w_case + 2.0*layer_spacing; +constexpr auto full_detector_height = full_module_height + steel_height + 2.0*layer_w_case + 1.0*layer_spacing; constexpr auto half_detector_height = 0.5L * full_detector_height; -constexpr double layer_z_displacement[8] = {-0.5*full_module_height + (20.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 0.5*layer_w_case, - -0.5*full_module_height + (20.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + layer_spacing + 1.5*layer_w_case, - -0.5*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 0.5*layer_w_case, - -0.5*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + layer_spacing_top + 1.5*layer_w_case, - -0.5*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 2*layer_spacing_top + 2.5*layer_w_case, - -0.5*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 3*layer_spacing_top + 3.5*layer_w_case, - -0.5*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 4*layer_spacing_top + 4.5*layer_w_case, - -0.5*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 5*layer_spacing_top + 5.5*layer_w_case}; - - -constexpr double module_beam_heights[8] = {20.0*m - 3*layer_w_case - 2*layer_spacing, +constexpr double layer_z_displacement[8] = {-0.5*full_module_height + (19.0*m ) + 0.5*layer_w_case, + -0.5*full_module_height + (19.0*m ) + layer_spacing + 1.5*layer_w_case, + -0.5*full_module_height + (25.0*m ) + 0.5*layer_w_case, + -0.5*full_module_height + (25.0*m ) + layer_spacing_top + 1.5*layer_w_case, + -0.5*full_module_height + (25.0*m ) + 2*layer_spacing_top + 2.5*layer_w_case, + -0.5*full_module_height + (25.0*m ) + 3*layer_spacing_top + 3.5*layer_w_case, + -0.5*full_module_height + (25.0*m ) + 4*layer_spacing_top + 4.5*layer_w_case, + -0.5*full_module_height + (25.0*m ) + 5*layer_spacing_top + 5.5*layer_w_case}; + +// The z coordinates of the BOTTOM of all layers in the world +constexpr double layer_z_world[10] = {19*m+layer_spacing+2.*layer_w_case, + 19*m+layer_w_case, + -(layer_z_displacement[0]+0.5*full_module_height-19*m-0.5*layer_w_case), + -(layer_z_displacement[1]+0.5*full_module_height-19*m-0.5*layer_w_case), + -(layer_z_displacement[2]+0.5*full_module_height-19*m-0.5*layer_w_case), + -(layer_z_displacement[3]+0.5*full_module_height-19*m-0.5*layer_w_case), + -(layer_z_displacement[4]+0.5*full_module_height-19*m-0.5*layer_w_case), + -(layer_z_displacement[5]+0.5*full_module_height-19*m-0.5*layer_w_case), + -(layer_z_displacement[6]+0.5*full_module_height-19*m-0.5*layer_w_case), + -(layer_z_displacement[7]+0.5*full_module_height-19*m-0.5*layer_w_case) +}; + + +constexpr double module_beam_heights[8] = {20.0*m, layer_spacing, - 5.0*m - 2*layer_w_case - layer_spacing, + 6.0*m - 2*layer_w_case - layer_spacing, layer_spacing_top, layer_spacing_top, layer_spacing_top, @@ -160,13 +173,13 @@ constexpr double module_beam_heights[8] = {20.0*m - 3*layer_w_case - 2*layer_spa // -0.50*full_module_height + 25.0*m + 4*layer_w_case + 3*layer_spacing + 0.50*module_beam_heights[8]}; constexpr double module_beam_z_pos[8] = {-0.50*full_module_height + 0.50*module_beam_heights[0], - -0.50*full_module_height + (20.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + layer_w_case + 0.50*module_beam_heights[1], - -0.50*full_module_height + (20.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 2*layer_w_case + layer_spacing + 0.50*module_beam_heights[2], - -0.50*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + layer_w_case + 0.50*module_beam_heights[3], - -0.50*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 2*layer_w_case + 1*layer_spacing_top + 0.50*module_beam_heights[4], - -0.50*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 3*layer_w_case + 2*layer_spacing_top + 0.50*module_beam_heights[5], - -0.50*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 4*layer_w_case + 3*layer_spacing_top + 0.50*module_beam_heights[6], - -0.50*full_module_height + (25.0*m - 3.0*layer_w_case - 2.0*layer_spacing) + 5*layer_w_case + 3*layer_spacing_top + 0.50*module_beam_heights[7]}; + -0.50*full_module_height + 20.0*m + layer_w_case + 0.50*module_beam_heights[1], + -0.50*full_module_height + 20.0*m + 2*layer_w_case + layer_spacing + 0.50*module_beam_heights[2], + -0.50*full_module_height + 25.0*m + layer_w_case + 0.50*module_beam_heights[3], + -0.50*full_module_height + 25.0*m + 2*layer_w_case + 1*layer_spacing_top + 0.50*module_beam_heights[4], + -0.50*full_module_height + 25.0*m + 3*layer_w_case + 2*layer_spacing_top + 0.50*module_beam_heights[5], + -0.50*full_module_height + 25.0*m + 4*layer_w_case + 3*layer_spacing_top + 0.50*module_beam_heights[6], + -0.50*full_module_height + 25.0*m + 5*layer_w_case + 3*layer_spacing_top + 0.50*module_beam_heights[7]}; @@ -218,9 +231,9 @@ Detector::Detector() : G4VSensitiveDetector("MATHUSLA/MU/Box") { scintillator->Register(this); // Print out the detector dislacement in z direction - std::cout<<"Layer z displacement:\n \n"<(std::floor(+local_position.y() / (layer_spacing + scintillator_height))); size_t y_index = 0; - if (new_position.y() < 6050.0L*cm) { +// if (new_position.y() < 6050.0L*cm) { +// y_index = static_cast(std::floor(+local_position.y() / (layer_spacing + scintillator_height))); +// } else if (new_position.y() > 6060.0L*cm && new_position.y() < 6150.0L*cm) { +// y_index = static_cast(std::floor((+local_position.y() - 1.0L*cm) / (layer_spacing + scintillator_height))); +// } else if (new_position.y() > 7900.0L*cm && new_position.y() < 8050.0L*cm) { +// y_index = static_cast(std::floor((+local_position.y() - 1796.0L*cm) / (layer_spacing + scintillator_height))); +// } else if (new_position.y() > 8050.0L*cm && new_position.y() < 8150.0L*cm) { +// y_index = static_cast(std::floor((+local_position.y() - 1797.0L*cm) / (layer_spacing + scintillator_height))); +// } else if (new_position.y() > 8400.0L*cm && new_position.y() < 8550.0L*cm) { +// y_index = static_cast(std::floor((+local_position.y() - 2092.0L*cm) / (layer_spacing + scintillator_height))); +// } else if (new_position.y() > 8560.0L*cm && new_position.y() < 8650.0L*cm) { +// y_index = static_cast(std::floor((+local_position.y() - 2093.0L*cm) / (layer_spacing + scintillator_height))); +// } else if (new_position.y() > 8660.0L*cm && new_position.y() < 8750.0L*cm) { +// y_index = static_cast(std::floor((+local_position.y() - 2094.0L*cm) / (layer_spacing + scintillator_height))); +// } else if (new_position.y() > 8760.0L*cm && new_position.y() < 8850.0L*cm) { +// y_index = static_cast(std::floor((+local_position.y() - 2095.0L*cm) / (layer_spacing + scintillator_height))); +// } else { +// y_index = static_cast(std::floor((+local_position.y() - 2096.0L*cm) / (layer_spacing + scintillator_height))); + // } + + if (new_position.y() <= 6003.5L*cm) { y_index = static_cast(std::floor(+local_position.y() / (layer_spacing + scintillator_height))); - } else if (new_position.y() > 6060.0L*cm && new_position.y() < 6150.0L*cm) { - y_index = static_cast(std::floor((+local_position.y() - 1.0L*cm) / (layer_spacing + scintillator_height))); - } else if (new_position.y() > 7900.0L*cm && new_position.y() < 8050.0L*cm) { - y_index = static_cast(std::floor((+local_position.y() - 1796.0L*cm) / (layer_spacing + scintillator_height))); - } else if (new_position.y() > 8050.0L*cm && new_position.y() < 8150.0L*cm) { - y_index = static_cast(std::floor((+local_position.y() - 1797.0L*cm) / (layer_spacing + scintillator_height))); - } else if (new_position.y() > 8400.0L*cm && new_position.y() < 8550.0L*cm) { - y_index = static_cast(std::floor((+local_position.y() - 2092.0L*cm) / (layer_spacing + scintillator_height))); - } else if (new_position.y() > 8560.0L*cm && new_position.y() < 8650.0L*cm) { - y_index = static_cast(std::floor((+local_position.y() - 2093.0L*cm) / (layer_spacing + scintillator_height))); - } else if (new_position.y() > 8660.0L*cm && new_position.y() < 8750.0L*cm) { - y_index = static_cast(std::floor((+local_position.y() - 2094.0L*cm) / (layer_spacing + scintillator_height))); - } else if (new_position.y() > 8760.0L*cm && new_position.y() < 8850.0L*cm) { - y_index = static_cast(std::floor((+local_position.y() - 2095.0L*cm) / (layer_spacing + scintillator_height))); - } else { + } else if (new_position.y() >= 6103.5L*cm && new_position.y() <= 6105.5L*cm) { + y_index = static_cast(std::floor((+local_position.y() - 100.0L*cm) / (layer_spacing + scintillator_height))); + } else if (new_position.y() >= 8101.5L*cm && new_position.y() <= 8103.5L*cm) { + y_index = static_cast(std::floor((+local_position.y() - 1996.0L*cm) / (layer_spacing + scintillator_height))); + } else if (new_position.y() >= 8203.5L*cm && new_position.y() <= 8205.5L*cm) { y_index = static_cast(std::floor((+local_position.y() - 2096.0L*cm) / (layer_spacing + scintillator_height))); + } else if (new_position.y() >= 8601.5L*cm && new_position.y() <= 8603.5L*cm) { + y_index = static_cast(std::floor((+local_position.y() - 2492.0L*cm) / (layer_spacing + scintillator_height))); + } else if (new_position.y() >= 8683.5L*cm && new_position.y() <= 8685.5L*cm) { + y_index = static_cast(std::floor((+local_position.y() - 2492.0L*cm +80.0L*cm) / (layer_spacing + scintillator_height))); + } else if (new_position.y() >= 8765.5L*cm && new_position.y() <= 8767.5L*cm) { + y_index = static_cast(std::floor((+local_position.y() - 2492.0L*cm +2*80.0L*cm) / (layer_spacing + scintillator_height))); + } else if (new_position.y() >= 8847.5L*cm && new_position.y() <= 8849.5L*cm) { + y_index = static_cast(std::floor((+local_position.y() - 2492.0L*cm +3*80.0L*cm) / (layer_spacing + scintillator_height))); + } else if (new_position.y() >= 8929.5L*cm && new_position.y() <= 8931.5L*cm) { + y_index = static_cast(std::floor((+local_position.y() - 2492.0L*cm +4*80.0L*cm) / (layer_spacing + scintillator_height))); + } else if (new_position.y() >= 9011.5L*cm && new_position.y() <= 9013.5L*cm) { + y_index = static_cast(std::floor((+local_position.y() - 2492.0L*cm +5*80.0L*cm) / (layer_spacing + scintillator_height))); + } else { + y_index = static_cast(std::floor((+local_position.y() - 2492.0L*cm +6*80.0L*cm) / (layer_spacing + scintillator_height))); } - // if (new_position.y() <= 6003.5L*cm) { - // y_index = static_cast(std::floor(+local_position.y() / (layer_spacing + scintillator_height))); - // } else if (new_position.y() >= 6103.5L*cm && new_position.y() <= 6105.5L*cm) { - // y_index = static_cast(std::floor((+local_position.y() - 100.0L*cm) / (layer_spacing + scintillator_height))); - // } else if (new_position.y() >= 8001.5L*cm && new_position.y() <= 8003.5L*cm) { - // y_index = static_cast(std::floor((+local_position.y() - 1996.0L*cm) / (layer_spacing + scintillator_height))); - // } else if (new_position.y() >= 8103.5L*cm && new_position.y() <= 8105.5L*cm) { - // y_index = static_cast(std::floor((+local_position.y() - 2096.0L*cm) / (layer_spacing + scintillator_height))); - // } else if (new_position.y() >= 8501.5L*cm && new_position.y() <= 8503.5L*cm) { - // y_index = static_cast(std::floor((+local_position.y() - 2492.0L*cm) / (layer_spacing + scintillator_height))); - // } else if (new_position.y() >= 8603.5L*cm && new_position.y() <= 8605.5L*cm) { - // y_index = static_cast(std::floor((+local_position.y() - 2592.0L*cm) / (layer_spacing + scintillator_height))); - // } else if (new_position.y() >= 8705.5L*cm && new_position.y() <= 8707.5L*cm) { - // y_index = static_cast(std::floor((+local_position.y() - 2692.0L*cm) / (layer_spacing + scintillator_height))); - // } else if (new_position.y() >= 8807.5L*cm && new_position.y() <= 8809.5L*cm) { - // y_index = static_cast(std::floor((+local_position.y() - 2792.0L*cm) / (layer_spacing + scintillator_height))); - // } else { - // y_index = static_cast(std::floor((+local_position.y() - 2892.0L*cm) / (layer_spacing + scintillator_height))); - // } - int _rotation = (1UL + y_index) % 2; size_t x_index; size_t z_index; @@ -554,7 +571,7 @@ G4VPhysicalVolume* Detector::ConstructModule(G4LogicalVolume* DetectorVolume, in return Construction::PlaceVolume(ModuleVolume, DetectorVolume, Construction::Transform(get_module_x_displacement(tag_number), get_module_y_displacement(tag_number), - half_detector_height - steel_height - 3.0*layer_w_case - 2.0*layer_spacing - 0.5*full_module_height, + half_detector_height - steel_height - 2.0*layer_w_case - 1*layer_spacing - 0.5*full_module_height, 0.0, 0.0, 1.0, 0.0)); @@ -597,13 +614,13 @@ G4VPhysicalVolume* Detector::Construct(G4LogicalVolume* world) { _scintillators.push_back(second_hermetic_floor); second_hermetic_floor->PlaceIn(DetectorVolume, G4Translate3D(0.0, 0.0, half_detector_height - 1.5*layer_w_case - layer_spacing - steel_height)); - auto third_hermetic_floor = new Scintillator("HF3", - x_edge_length, - y_edge_length, - full_layer_height, - scintillator_casing_thickness); - _scintillators.push_back(third_hermetic_floor); - third_hermetic_floor->PlaceIn(DetectorVolume, G4Translate3D(0.0, 0.0, half_detector_height - 2.5*layer_w_case - 2*layer_spacing - steel_height)); + // auto third_hermetic_floor = new Scintillator("HF3", + // x_edge_length, + // y_edge_length, + // full_layer_height, + // scintillator_casing_thickness); + // _scintillators.push_back(third_hermetic_floor); + // third_hermetic_floor->PlaceIn(DetectorVolume, G4Translate3D(0.0, 0.0, half_detector_height - 2.5*layer_w_case - 2*layer_spacing - steel_height)); auto hermetic_wall = new Scintillator("HW1", full_layer_height, @@ -622,7 +639,7 @@ G4VPhysicalVolume* Detector::Construct(G4LogicalVolume* world) { // Construction::Export(DetectorVolume, folder, file, arg4 ); return Construction::PlaceVolume(DetectorVolume, world, - Construction::Transform(0.5L*x_edge_length + x_displacement, 0.5L*y_edge_length + y_displacement, -0.50*full_detector_height + 20*m)); + Construction::Transform(0.5L*x_edge_length + x_displacement, 0.5L*y_edge_length + y_displacement, -0.50*full_detector_height + 19*m + layer_spacing + 2.*layer_w_case + steel_height)); } From 9d404ed052c5539c1f6dd17b93c85ca8c14b7be5 Mon Sep 17 00:00:00 2001 From: Runze Ren Date: Tue, 2 May 2023 15:08:28 -0700 Subject: [PATCH 06/12] change all layer gap to 80cm; save layer position --- include/geometry/Box.hh | 1 + include/geometry/Construction.hh | 2 + src/action/RunAction.cc | 2 + src/geometry/Construction.cc | 9 +++ src/geometry/box/Box.cc | 99 ++++++++++++++++++++++---------- 5 files changed, 83 insertions(+), 30 deletions(-) diff --git a/include/geometry/Box.hh b/include/geometry/Box.hh index 1d9b2de..8e2459a 100644 --- a/include/geometry/Box.hh +++ b/include/geometry/Box.hh @@ -87,6 +87,7 @@ public: static bool SaveAll; static void WritePreData(); + static void SaveInfo(std::string & prefix); }; } /* namespace Box */ ////////////////////////////////////////////////////////////////////////// diff --git a/include/geometry/Construction.hh b/include/geometry/Construction.hh index 995a2ea..d9ed6d1 100644 --- a/include/geometry/Construction.hh +++ b/include/geometry/Construction.hh @@ -62,6 +62,8 @@ public: static void SetDetector(const std::string& detector); static void SetSaveOption(const bool save_option, const bool cut_save_option); + + static void SaveInfo(std::string& prefix); static const std::string& GetDetectorName(); static bool IsDetectorDataPerEvent(); diff --git a/src/action/RunAction.cc b/src/action/RunAction.cc index 43ca57b..84a936e 100644 --- a/src/action/RunAction.cc +++ b/src/action/RunAction.cc @@ -125,6 +125,8 @@ void RunAction::BeginOfRunAction(const G4Run* run) { Construction::Builder::GetDetectorDataKeys(), Construction::Builder::GetDetectorDataKeyTypes()); +Construction::Builder::SaveInfo(_prefix); + if (!G4Threading::IsWorkerThread()) std::cout << "\n\n"; } diff --git a/src/geometry/Construction.cc b/src/geometry/Construction.cc index 4a5470a..38d7594 100644 --- a/src/geometry/Construction.cc +++ b/src/geometry/Construction.cc @@ -284,6 +284,15 @@ void Builder::SetSaveOption(bool save_option, bool cut_save_option) { } //---------------------------------------------------------------------------------------------- +//__Save some information in a txt file____________________________________________________________ +void Builder::SaveInfo(std::string & prefix){ + if (_detector == "Flat") {;} + else if (_detector == "Box") { + Box::Detector::SaveInfo(prefix); + } +} +//---------------------------------------------------------------------------------------------- + //__Get Current Detector Name___________________________________________________________________ const std::string& Builder::GetDetectorName() { return _detector; diff --git a/src/geometry/box/Box.cc b/src/geometry/box/Box.cc index 2ae00e4..6784ede 100644 --- a/src/geometry/box/Box.cc +++ b/src/geometry/box/Box.cc @@ -16,7 +16,11 @@ #include "TTree.h" #include "MuonDataController.hh" +#include +#include +#include +#define UNUSED(x) (void)(x) using dimension = double; namespace MATHUSLA { namespace MU { @@ -102,7 +106,7 @@ constexpr auto air_gap = 30*m; constexpr auto scintillator_casing_thickness = 0.003*m; -constexpr auto layer_spacing = 1.0*m; +constexpr auto layer_spacing = 0.8*m; constexpr auto layer_spacing_top = 0.8*m; // Top layers have different spacing constexpr auto layer_count = 8UL; @@ -130,8 +134,8 @@ constexpr auto beam_thickness = 0.02*m; constexpr auto full_detector_height = full_module_height + steel_height + 2.0*layer_w_case + 1.0*layer_spacing; constexpr auto half_detector_height = 0.5L * full_detector_height; -constexpr double layer_z_displacement[8] = {-0.5*full_module_height + (19.0*m ) + 0.5*layer_w_case, - -0.5*full_module_height + (19.0*m ) + layer_spacing + 1.5*layer_w_case, +constexpr double layer_z_displacement[8] = {-0.5*full_module_height + (25.0*m -5.0*m - layer_spacing - 2*layer_w_case) + 0.5*layer_w_case, + -0.5*full_module_height + (25.0*m -5.0*m - layer_spacing - 2*layer_w_case) + layer_spacing + 1.5*layer_w_case, -0.5*full_module_height + (25.0*m ) + 0.5*layer_w_case, -0.5*full_module_height + (25.0*m ) + layer_spacing_top + 1.5*layer_w_case, -0.5*full_module_height + (25.0*m ) + 2*layer_spacing_top + 2.5*layer_w_case, @@ -139,23 +143,25 @@ constexpr double layer_z_displacement[8] = {-0.5*full_module_height + (19.0*m ) -0.5*full_module_height + (25.0*m ) + 4*layer_spacing_top + 4.5*layer_w_case, -0.5*full_module_height + (25.0*m ) + 5*layer_spacing_top + 5.5*layer_w_case}; +// The z distance from the TOP of the floor layer to the BOTTOM of the middle layer +constexpr double floor_z_top = 25.0*m -5.0*m - layer_spacing - 2*layer_w_case; // The z coordinates of the BOTTOM of all layers in the world -constexpr double layer_z_world[10] = {19*m+layer_spacing+2.*layer_w_case, - 19*m+layer_w_case, - -(layer_z_displacement[0]+0.5*full_module_height-19*m-0.5*layer_w_case), - -(layer_z_displacement[1]+0.5*full_module_height-19*m-0.5*layer_w_case), - -(layer_z_displacement[2]+0.5*full_module_height-19*m-0.5*layer_w_case), - -(layer_z_displacement[3]+0.5*full_module_height-19*m-0.5*layer_w_case), - -(layer_z_displacement[4]+0.5*full_module_height-19*m-0.5*layer_w_case), - -(layer_z_displacement[5]+0.5*full_module_height-19*m-0.5*layer_w_case), - -(layer_z_displacement[6]+0.5*full_module_height-19*m-0.5*layer_w_case), - -(layer_z_displacement[7]+0.5*full_module_height-19*m-0.5*layer_w_case) +constexpr double layer_z_world[10] = {floor_z_top +layer_spacing+2.*layer_w_case, + floor_z_top +layer_w_case, + -(layer_z_displacement[0]+0.5*full_module_height-floor_z_top-0.5*layer_w_case), + -(layer_z_displacement[1]+0.5*full_module_height-floor_z_top-0.5*layer_w_case), + -(layer_z_displacement[2]+0.5*full_module_height-floor_z_top-0.5*layer_w_case), + -(layer_z_displacement[3]+0.5*full_module_height-floor_z_top-0.5*layer_w_case), + -(layer_z_displacement[4]+0.5*full_module_height-floor_z_top-0.5*layer_w_case), + -(layer_z_displacement[5]+0.5*full_module_height-floor_z_top-0.5*layer_w_case), + -(layer_z_displacement[6]+0.5*full_module_height-floor_z_top-0.5*layer_w_case), + -(layer_z_displacement[7]+0.5*full_module_height-floor_z_top-0.5*layer_w_case) }; -constexpr double module_beam_heights[8] = {20.0*m, +constexpr double module_beam_heights[8] = {floor_z_top, layer_spacing, - 6.0*m - 2*layer_w_case - layer_spacing, + 5.0*m, layer_spacing_top, layer_spacing_top, layer_spacing_top, @@ -173,13 +179,13 @@ constexpr double module_beam_heights[8] = {20.0*m, // -0.50*full_module_height + 25.0*m + 4*layer_w_case + 3*layer_spacing + 0.50*module_beam_heights[8]}; constexpr double module_beam_z_pos[8] = {-0.50*full_module_height + 0.50*module_beam_heights[0], - -0.50*full_module_height + 20.0*m + layer_w_case + 0.50*module_beam_heights[1], - -0.50*full_module_height + 20.0*m + 2*layer_w_case + layer_spacing + 0.50*module_beam_heights[2], + -0.50*full_module_height + floor_z_top + layer_w_case + 0.50*module_beam_heights[1], + -0.50*full_module_height + floor_z_top + 2*layer_w_case + layer_spacing + 0.50*module_beam_heights[2], -0.50*full_module_height + 25.0*m + layer_w_case + 0.50*module_beam_heights[3], -0.50*full_module_height + 25.0*m + 2*layer_w_case + 1*layer_spacing_top + 0.50*module_beam_heights[4], -0.50*full_module_height + 25.0*m + 3*layer_w_case + 2*layer_spacing_top + 0.50*module_beam_heights[5], -0.50*full_module_height + 25.0*m + 4*layer_w_case + 3*layer_spacing_top + 0.50*module_beam_heights[6], - -0.50*full_module_height + 25.0*m + 5*layer_w_case + 3*layer_spacing_top + 0.50*module_beam_heights[7]}; + -0.50*full_module_height + 25.0*m + 5*layer_w_case + 4*layer_spacing_top + 0.50*module_beam_heights[7]}; @@ -231,14 +237,39 @@ Detector::Detector() : G4VSensitiveDetector("MATHUSLA/MU/Box") { scintillator->Register(this); // Print out the detector dislacement in z direction - std::cout<<"Layer z displacement [mm]:\n world coord.---- CMS coord.\n"<PlaceIn(ModuleVolume, G4Translate3D(0.0, 0.0, layer_z_displacement) ); //G4Translate3D(0.5*layer_x_edge_length, 0.5*layer_y_edge_length, layer_z_displacement) } @@ -533,8 +567,13 @@ G4VPhysicalVolume* Detector::ConstructModule(G4LogicalVolume* DetectorVolume, in 0*m, 0*m, get_layer_z_displacement(layer)); + UNUSED(current); } + UNUSED(detector_x); + UNUSED(detector_y); + UNUSED(detector_z); + //CONSTRUCTING AND INSERTING STEEL BEAMS for (int beam_layer = 0; beam_layer < NBEAMLAYERS; beam_layer++){ @@ -587,15 +626,15 @@ G4VPhysicalVolume* Detector::Construct(G4LogicalVolume* world) { // pre_data->Branch("Y_H", &Y_POS_HIT, "Y_H/D"); auto DetectorVolume = Construction::BoxVolume("Box", x_edge_length + x_edge_increase, y_edge_length, full_detector_height, - Construction::Material::Air, G4VisAttributes::Invisible); + Construction::Material::Air, G4VisAttributes::GetInvisible()); //DetectorVolume->SetVisAttributes(G4VisAttributes::Invisible); for (int module_number = 0; module_number < NMODULES; module_number++){ - auto current = Detector::ConstructModule(DetectorVolume, module_number, - 0.5L*x_edge_length + x_displacement, //add extra terms for displacement from center here - 0.5L*y_edge_length + y_displacement, - -half_detector_height + steel_height + 3.0*layer_w_case + 2.0*layer_spacing); + Detector::ConstructModule(DetectorVolume, module_number, + 0.5L, //These arguments are not actually used + 0.5L, + 0.0L); } auto first_hermetic_floor = new Scintillator("HF1", @@ -639,7 +678,7 @@ G4VPhysicalVolume* Detector::Construct(G4LogicalVolume* world) { // Construction::Export(DetectorVolume, folder, file, arg4 ); return Construction::PlaceVolume(DetectorVolume, world, - Construction::Transform(0.5L*x_edge_length + x_displacement, 0.5L*y_edge_length + y_displacement, -0.50*full_detector_height + 19*m + layer_spacing + 2.*layer_w_case + steel_height)); + Construction::Transform(0.5L*x_edge_length + x_displacement, 0.5L*y_edge_length + y_displacement, -0.50*full_detector_height + floor_z_top + layer_spacing + 2.*layer_w_case + steel_height)); } @@ -790,7 +829,7 @@ G4VPhysicalVolume* Detector::ConstructEarth(G4LogicalVolume* world){ auto CMS_Detector_logical = CMSRingVolume(); auto UXC_55_air_v1 = new G4SubtractionSolid("UXC_55_air_v1", UXC_55_cavern_solid, UXC55_outer_solid); auto UXC_55_air_v2 = new G4SubtractionSolid("UXC_55_air_v2", UXC_55_air_v1, CMS_Detector_logical->GetSolid()); - auto UXC55_air_logical = Volume("UXC55_air", UXC_55_air_v2, Construction::Material::Air, G4VisAttributes::Invisible); + auto UXC55_air_logical = Volume("UXC55_air", UXC_55_air_v2, Construction::Material::Air, G4VisAttributes::GetInvisible()); Construction::PlaceVolume(UXC55_outer_logical, earth, Cavern_Transform()*Construction::Rotate(0, 1, 0, 90*deg) ); Construction::PlaceVolume(CMS_Detector_logical, earth, Cavern_Transform()*Construction::Rotate(0, 1, 0, 90*deg) ); @@ -814,7 +853,7 @@ G4VPhysicalVolume* Detector::ConstructEarth(G4LogicalVolume* world){ AS_Depth - 2* AS_Thickness, AS_Height - 2* AS_Thickness, Construction::Material::Air, - G4VisAttributes::Invisible); + G4VisAttributes::GetInvisible()); Construction::PlaceVolume(Access_Shaft_outer_logical, earth, Access_Shaft_Transform() ); Construction::PlaceVolume(Access_Shaft_Air, earth, Access_Shaft_Transform()); From 14012d4e24fc55e0aaad87abde1684edf9277636 Mon Sep 17 00:00:00 2001 From: Abdulrahman Mohamed <88555369+Abdulrahma-M@users.noreply.github.com> Date: Mon, 8 May 2023 20:35:52 +0300 Subject: [PATCH 07/12] 0.8m seperation --- Box.cc | 906 ++++++++++++++++++++++++++++++++++++++++++++++++ Box.hh | 97 ++++++ Construction.cc | 667 +++++++++++++++++++++++++++++++++++ Construction.hh | 294 ++++++++++++++++ 4 files changed, 1964 insertions(+) create mode 100644 Box.cc create mode 100644 Box.hh create mode 100644 Construction.cc create mode 100644 Construction.hh diff --git a/Box.cc b/Box.cc new file mode 100644 index 0000000..a1f2cf0 --- /dev/null +++ b/Box.cc @@ -0,0 +1,906 @@ +#include "geometry/Box.hh" + +#include +#include + +#include "action.hh" +#include "analysis.hh" +#include "geometry/Earth.hh" +#include "physics/Units.hh" +#include "tracking.hh" +#include "geometry/Cavern.hh" +#include +#include +#include +#include "TROOT.h" +#include "TTree.h" +#include "MuonDataController.hh" + +#include +#include +#include + +#define UNUSED(x) (void)(x) +using dimension = double; + +namespace MATHUSLA { namespace MU { + +namespace Box { ////////////////////////////////////////////////////////////////////////////////////////////////////// + +// double X_POS_STEP; +// double Y_POS_STEP; +// double X_POS_HIT; +// double Y_POS_HIT; + +// TTree* Detector::pre_data = new TTree("pre_data", "pre_data"); + +// void Detector::WritePreData(){ + +// TFile* f = new TFile("box_pre_data.root", "RECREATE"); +// f->cd(); + +// pre_data->Write(); + +// f->Write(); + +//} + +namespace Box::CavernConstruction{ //////////////////////////////////////////////////////////////////////////////////// + +//__Check Between-Ness__________________________________________________________________________ +bool _between(const double min_layer, + const double max_layer, + const double target) { + return min_layer < target && target < max_layer; +} +//---------------------------------------------------------------------------------------------- + +//__Calculate Subtraction of Volumes____________________________________________________________ +G4LogicalVolume* _calculate_modification(const std::string& name, + G4LogicalVolume* earth_component, + const double base_depth, + const double top_depth) { + return Construction::Volume(new G4SubtractionSolid(name, + earth_component->GetSolid(), + Cavern::Volume()->GetSolid(), + Construction::Transform(0, 1.7 * m, -0.5 * (base_depth - top_depth) + Cavern::CenterDepth() - top_depth)), + earth_component->GetMaterial()); +} + +} + + + +namespace { //////////////////////////////////////////////////////////////////////////////////// + +//__Box Sensitive Material______________________________________________________________________ +std::vector _scintillators; +G4LogicalVolume* _steel; +//---------------------------------------------------------------------------------------------- + +//__Box Hit Collection__________________________________________________________________________ +G4ThreadLocal Tracking::HitCollection* _hit_collection; +//---------------------------------------------------------------------------------------------- + +//__Box Specification Variables_________________________________________________________________ + +constexpr int scintillators_per_layer{400}; +constexpr int NMODULES{100}; +constexpr int n_top_layers{6}; +constexpr auto x_edge_length = 99.0*m; +constexpr auto y_edge_length = 99.0*m; +constexpr auto x_displacement = 70.0*m; +constexpr auto y_displacement = -49.5*m; +constexpr auto z_displacement = 6001.5*cm; + +constexpr auto layer_x_edge_length = 9.0*m; +constexpr auto layer_y_edge_length = 9.0*m; + +constexpr auto scint_x_edge_length = 4.5*m; +constexpr auto scint_y_edge_length = 0.045*m; +constexpr auto scintillator_height = 0.02*m; + +constexpr auto steel_height = 0.03*m; + +constexpr auto air_gap = 30*m; + +constexpr auto scintillator_casing_thickness = 0.003*m; + +constexpr auto layer_spacing = 0.8*m; +constexpr auto layer_spacing_top = 0.8*m; // Top layers have different spacing +constexpr auto layer_count = 8UL; + +constexpr auto module_x_edge_length = 9.0*m; +constexpr auto module_y_edge_length = 9.0*m; +constexpr auto module_case_thickness = 0.02*m; + +constexpr auto full_layer_height = scintillator_height + 2*scintillator_casing_thickness; +constexpr auto wall_gap = 0.01*m; +constexpr auto x_edge_increase = 2*full_layer_height + 4*wall_gap; + +constexpr auto layer_w_case = full_layer_height; + +constexpr auto full_module_height = (25.0*m) + 6.0*layer_w_case + 5.0*layer_spacing_top; + +constexpr auto scintillator_z_position = 0.00; + +constexpr auto wall_height = 20*m; + +constexpr int NBEAMLAYERS = 8; +constexpr auto beam_x_edge_length = 0.10*m; +constexpr auto beam_y_edge_length = 0.10*m; +constexpr auto beam_thickness = 0.02*m; + +constexpr auto full_detector_height = full_module_height + steel_height + 2.0*layer_w_case + 1.0*layer_spacing; +constexpr auto half_detector_height = 0.5L * full_detector_height; + +constexpr double layer_z_displacement[8] = {-0.5*full_module_height + (25.0*m -5.0*m - layer_spacing - 2*layer_w_case) + 0.5*layer_w_case, + -0.5*full_module_height + (25.0*m -5.0*m - layer_spacing - 2*layer_w_case) + layer_spacing + 1.5*layer_w_case, + -0.5*full_module_height + (25.0*m ) + 0.5*layer_w_case, + -0.5*full_module_height + (25.0*m ) + layer_spacing_top + 1.5*layer_w_case, + -0.5*full_module_height + (25.0*m ) + 2*layer_spacing_top + 2.5*layer_w_case, + -0.5*full_module_height + (25.0*m ) + 3*layer_spacing_top + 3.5*layer_w_case, + -0.5*full_module_height + (25.0*m ) + 4*layer_spacing_top + 4.5*layer_w_case, + -0.5*full_module_height + (25.0*m ) + 5*layer_spacing_top + 5.5*layer_w_case}; + +// The z distance from the TOP of the floor layer to the BOTTOM of the middle layer +constexpr double floor_z_top = 25.0*m -5.0*m - layer_spacing - 2*layer_w_case; +// The z coordinates of the BOTTOM of all layers in the world +constexpr double layer_z_world[10] = {floor_z_top +layer_spacing+2.*layer_w_case, + floor_z_top +layer_w_case, + -(layer_z_displacement[0]+0.5*full_module_height-floor_z_top-0.5*layer_w_case), + -(layer_z_displacement[1]+0.5*full_module_height-floor_z_top-0.5*layer_w_case), + -(layer_z_displacement[2]+0.5*full_module_height-floor_z_top-0.5*layer_w_case), + -(layer_z_displacement[3]+0.5*full_module_height-floor_z_top-0.5*layer_w_case), + -(layer_z_displacement[4]+0.5*full_module_height-floor_z_top-0.5*layer_w_case), + -(layer_z_displacement[5]+0.5*full_module_height-floor_z_top-0.5*layer_w_case), + -(layer_z_displacement[6]+0.5*full_module_height-floor_z_top-0.5*layer_w_case), + -(layer_z_displacement[7]+0.5*full_module_height-floor_z_top-0.5*layer_w_case) +}; + + +constexpr double module_beam_heights[8] = {floor_z_top, + layer_spacing, + 5.0*m, + layer_spacing_top, + layer_spacing_top, + layer_spacing_top, + layer_spacing_top, + layer_spacing_top}; + + // constexpr double module_beam_z_pos[9] = {-0.50*full_module_height + 0.50*module_beam_heights[0] + layer_w_case, + // -0.50*full_module_height + 2*layer_w_case + layer_spacing + 0.50*module_beam_heights[1], + // -0.50*full_module_height + 3*layer_w_case + 2*layer_spacing + 0.50*module_beam_heights[2], + // -0.50*full_module_height + 20.0*m + layer_w_case + 0.50*module_beam_heights[3], + // -0.50*full_module_height + 20.0*m + 2*layer_w_case + layer_spacing + 0.50*module_beam_heights[4], + // -0.50*full_module_height + 25.0*m + layer_w_case + 0.50*module_beam_heights[5], + // -0.50*full_module_height + 25.0*m + 2*layer_w_case + layer_spacing + 0.50*module_beam_heights[6], + // -0.50*full_module_height + 25.0*m + 3*layer_w_case + 2*layer_spacing + 0.50*module_beam_heights[7], + // -0.50*full_module_height + 25.0*m + 4*layer_w_case + 3*layer_spacing + 0.50*module_beam_heights[8]}; + +constexpr double module_beam_z_pos[8] = {-0.50*full_module_height + 0.50*module_beam_heights[0], + -0.50*full_module_height + floor_z_top + layer_w_case + 0.50*module_beam_heights[1], + -0.50*full_module_height + floor_z_top + 2*layer_w_case + layer_spacing + 0.50*module_beam_heights[2], + -0.50*full_module_height + 25.0*m + layer_w_case + 0.50*module_beam_heights[3], + -0.50*full_module_height + 25.0*m + 2*layer_w_case + 1*layer_spacing_top + 0.50*module_beam_heights[4], + -0.50*full_module_height + 25.0*m + 3*layer_w_case + 2*layer_spacing_top + 0.50*module_beam_heights[5], + -0.50*full_module_height + 25.0*m + 4*layer_w_case + 3*layer_spacing_top + 0.50*module_beam_heights[6], + -0.50*full_module_height + 25.0*m + 5*layer_w_case + 4*layer_spacing_top + 0.50*module_beam_heights[7]}; + + + +// const std::string folder = "detector_geo"; +// const std::string file = "box.gdml"; +// const std::string file2 ="mod.gdml"; +// const std::string file3 ="layer.gdml"; +// const std::string file4 ="earth.gdml"; +// const std::string file5 ="modified.gdml"; +// const std::string arg4 = "http://service-spi.web.cern.ch/service-spi/app/releases/GDML/Schema/gdml.xsd"; + + +auto get_module_x_displacement(int tag_number){ + if (tag_number < 10) return -0.5 * x_edge_length + 0.5*module_x_edge_length; + else if (tag_number < 20) return -0.5*x_edge_length + 0.5*module_x_edge_length + 1.00*(module_x_edge_length + 1.0*m ); + else if (tag_number < 30) return -0.5*x_edge_length + 0.5*module_x_edge_length + 2.00*(module_x_edge_length + 1.0*m ); + else if (tag_number < 40) return -0.5*x_edge_length + 0.5*module_x_edge_length + 3.00*(module_x_edge_length + 1.0*m ); + else if (tag_number < 50) return -0.5*x_edge_length + 0.5*module_x_edge_length + 4.00*(module_x_edge_length + 1.0*m ); + else if (tag_number < 60) return -0.5*x_edge_length + 0.5*module_x_edge_length + 5.00*(module_x_edge_length + 1.0*m ); + else if (tag_number < 70) return -0.5*x_edge_length + 0.5*module_x_edge_length + 6.00*(module_x_edge_length + 1.0*m ); + else if (tag_number < 80) return -0.5*x_edge_length + 0.5*module_x_edge_length + 7.00*(module_x_edge_length + 1.0*m ); + else if (tag_number < 90) return -0.5*x_edge_length + 0.5*module_x_edge_length + 8.00*(module_x_edge_length + 1.0*m ); + else return -0.5*x_edge_length + 0.5*module_x_edge_length + 9.00*(module_x_edge_length + 1.0*m ); +} + +auto get_module_y_displacement(int tag_number){ + return (((double) (tag_number % 10))*(module_y_edge_length + 1.0*m) -0.5 * y_edge_length + 0.5*module_y_edge_length ); +} + +auto get_layer_z_displacement(int layer_number){ + return -1.0*layer_z_displacement[layer_number]; +} + +//---------------------------------------------------------------------------------------------- + +} /* anonymous namespace */ //////////////////////////////////////////////////////////////////// + +//__Box Data Variables__________________________________________________________________________ +const std::string& Detector::DataName = "box_run"; +const Analysis::ROOT::DataKeyList Detector::DataKeys = Analysis::ROOT::DefaultDataKeyList; +const Analysis::ROOT::DataKeyTypeList Detector::DataKeyTypes = Analysis::ROOT::DefaultDataKeyTypeList; +bool Detector::SaveAll = false; +//---------------------------------------------------------------------------------------------- + +//__Detector Constructor________________________________________________________________________ +Detector::Detector() : G4VSensitiveDetector("MATHUSLA/MU/Box") { + collectionName.insert("Box_HC"); + for (auto& scintillator : _scintillators) + scintillator->Register(this); + + // Print out the detector dislacement in z direction + std::cout<<"Layer z displacement [cm]:\n world coord.---- CMS coord.\n"<GetTotalEnergyDeposit(); + + //const auto step_point = step->GetPreStepPoint(); + //const auto position = G4LorentzVector(step_point->GetGlobalTime(), step_point->GetPosition()); + // X_POS_STEP = position.x(); + // Y_POS_STEP = position.y(); + // X_POS_HIT = 0.0; + // Y_POS_HIT = 0.0; + + if (deposit == 0.0L){ + // pre_data->Fill(); + return false; + } + + // X_POS_HIT = position.x(); + // Y_POS_HIT = position.y(); + // pre_data->Fill(); + // const auto track = step->GetTrack(); + // const auto particle = track->GetParticleDefinition(); + // const auto trackID = track->GetTrackID(); + // const auto parentID = track->GetParentID(); + // const auto momentum = G4LorentzVector(step_point->GetTotalEnergy(), step_point->GetMomentum()); + +//////////////////////// + const auto track = step->GetTrack(); + const auto step_point = step->GetPostStepPoint(); + const auto particle = track->GetParticleDefinition(); + const auto trackID = track->GetTrackID(); + const auto parentID = track->GetParentID(); + const auto position = G4LorentzVector(step_point->GetGlobalTime(), step_point->GetPosition()); + const auto momentum = G4LorentzVector(step_point->GetTotalEnergy(), step_point->GetMomentum()); + + //______Tranfomation to CMS Coordinates_____________________________________________________ + const auto transformed_z = -(position.z() - Box_IP_Depth); + const auto position_transformed = G4ThreeVector(position.y(), transformed_z, position.x()); + const auto new_position = G4LorentzVector(step_point->GetGlobalTime(), position_transformed); + + const auto transformed_pz = -(momentum.pz()); + const auto momentum_transformed = G4ThreeVector(momentum.py(), transformed_pz, momentum.px()); + const auto new_momentum = G4LorentzVector(step_point->GetTotalEnergy(), momentum_transformed); + //__________________________________________________________________________________________ + + const auto local_position = new_position.vect() - G4ThreeVector(y_displacement, z_displacement, x_displacement); + + // auto z_index = static_cast(std::floor(+local_position.y() / (layer_spacing + scintillator_height))); + size_t y_index = 0; + +// if (new_position.y() < 6050.0L*cm) { +// y_index = static_cast(std::floor(+local_position.y() / (layer_spacing + scintillator_height))); +// } else if (new_position.y() > 6060.0L*cm && new_position.y() < 6150.0L*cm) { +// y_index = static_cast(std::floor((+local_position.y() - 1.0L*cm) / (layer_spacing + scintillator_height))); +// } else if (new_position.y() > 7900.0L*cm && new_position.y() < 8050.0L*cm) { +// y_index = static_cast(std::floor((+local_position.y() - 1796.0L*cm) / (layer_spacing + scintillator_height))); +// } else if (new_position.y() > 8050.0L*cm && new_position.y() < 8150.0L*cm) { +// y_index = static_cast(std::floor((+local_position.y() - 1797.0L*cm) / (layer_spacing + scintillator_height))); +// } else if (new_position.y() > 8400.0L*cm && new_position.y() < 8550.0L*cm) { +// y_index = static_cast(std::floor((+local_position.y() - 2092.0L*cm) / (layer_spacing + scintillator_height))); +// } else if (new_position.y() > 8560.0L*cm && new_position.y() < 8650.0L*cm) { +// y_index = static_cast(std::floor((+local_position.y() - 2093.0L*cm) / (layer_spacing + scintillator_height))); +// } else if (new_position.y() > 8660.0L*cm && new_position.y() < 8750.0L*cm) { +// y_index = static_cast(std::floor((+local_position.y() - 2094.0L*cm) / (layer_spacing + scintillator_height))); +// } else if (new_position.y() > 8760.0L*cm && new_position.y() < 8850.0L*cm) { +// y_index = static_cast(std::floor((+local_position.y() - 2095.0L*cm) / (layer_spacing + scintillator_height))); +// } else { +// y_index = static_cast(std::floor((+local_position.y() - 2096.0L*cm) / (layer_spacing + scintillator_height))); + // } + + if (new_position.y() <= 6003.5L*cm) { + y_index = static_cast(std::floor(+local_position.y() / (layer_spacing + scintillator_height))); + } else if (new_position.y() >= 6103.5L*cm && new_position.y() <= 6105.5L*cm) { + y_index = static_cast(std::floor((+local_position.y() - 100.0L*cm) / (layer_spacing + scintillator_height))); + } else if (new_position.y() >= 8101.5L*cm && new_position.y() <= 8103.5L*cm) { + y_index = static_cast(std::floor((+local_position.y() - 1996.0L*cm) / (layer_spacing + scintillator_height))); + } else if (new_position.y() >= 8203.5L*cm && new_position.y() <= 8205.5L*cm) { + y_index = static_cast(std::floor((+local_position.y() - 2096.0L*cm) / (layer_spacing + scintillator_height))); + } else if (new_position.y() >= 8601.5L*cm && new_position.y() <= 8603.5L*cm) { + y_index = static_cast(std::floor((+local_position.y() - 2492.0L*cm) / (layer_spacing + scintillator_height))); + } else if (new_position.y() >= 8683.5L*cm && new_position.y() <= 8685.5L*cm) { + y_index = static_cast(std::floor((+local_position.y() - 2492.0L*cm +80.0L*cm) / (layer_spacing + scintillator_height))); + } else if (new_position.y() >= 8765.5L*cm && new_position.y() <= 8767.5L*cm) { + y_index = static_cast(std::floor((+local_position.y() - 2492.0L*cm +2*80.0L*cm) / (layer_spacing + scintillator_height))); + } else if (new_position.y() >= 8847.5L*cm && new_position.y() <= 8849.5L*cm) { + y_index = static_cast(std::floor((+local_position.y() - 2492.0L*cm +3*80.0L*cm) / (layer_spacing + scintillator_height))); + } else if (new_position.y() >= 8929.5L*cm && new_position.y() <= 8931.5L*cm) { + y_index = static_cast(std::floor((+local_position.y() - 2492.0L*cm +4*80.0L*cm) / (layer_spacing + scintillator_height))); + } else if (new_position.y() >= 9011.5L*cm && new_position.y() <= 9013.5L*cm) { + y_index = static_cast(std::floor((+local_position.y() - 2492.0L*cm +5*80.0L*cm) / (layer_spacing + scintillator_height))); + } else { + y_index = static_cast(std::floor((+local_position.y() - 2492.0L*cm +6*80.0L*cm) / (layer_spacing + scintillator_height))); + } + + int _rotation = (1UL + y_index) % 2; + size_t x_index = 0; + size_t z_index = 0; + + if (_rotation == 0){ + + if (new_position.x() <= -4050.0L*cm) { + x_index = static_cast(std::floor(+local_position.x() / (scint_y_edge_length) )); + } else if (new_position.x() >= -3950.0L*cm && new_position.x() <= -3050.0L*cm) { + x_index = static_cast(std::floor((+local_position.x() - 100.0L*cm) / (scint_y_edge_length) )); + } else if (new_position.x() >= -2950.0L*cm && new_position.x() <= -2050.0L*cm) { + x_index = static_cast(std::floor((+local_position.x() - 200.0L*cm) / (scint_y_edge_length) )); + } else if (new_position.x() >= -1950.0L*cm && new_position.x() <= -1050.0L*cm) { + x_index = static_cast(std::floor((+local_position.x() - 300.0L*cm) / (scint_y_edge_length) )); + } else if (new_position.x() >= -950.0L*cm && new_position.x() <= -50.0L*cm) { + x_index = static_cast(std::floor((+local_position.x() - 400.0L*cm) / (scint_y_edge_length) )); + } else if (new_position.x() >= 50.0L*cm && new_position.x() <= 950.0L*cm) { + x_index = static_cast(std::floor((+local_position.x() - 500.0L*cm) / (scint_y_edge_length) )); + } else if (new_position.x() >= 1050.0L*cm && new_position.x() <= 1950.0L*cm) { + x_index = static_cast(std::floor((+local_position.x() - 600.0L*cm) / (scint_y_edge_length) )); + } else if (new_position.x() >= 2050.0L*cm && new_position.x() <= 2950.0L*cm) { + x_index = static_cast(std::floor((+local_position.x() - 700.0L*cm) / (scint_y_edge_length) )); + } else if (new_position.x() >= 3050.0L*cm && new_position.x() <= 3950.0L*cm) { + x_index = static_cast(std::floor((+local_position.x() - 800.0L*cm) / (scint_y_edge_length) )); + } else if (new_position.x() >= 4050.0L*cm && new_position.x() <= 4950.0L*cm) { + x_index = static_cast(std::floor((+local_position.x() - 900.0L*cm) / (scint_y_edge_length) )); + } + if (new_position.z() <= 7900.0L*cm) { + z_index = static_cast(std::floor(+local_position.z() / (scint_x_edge_length) )); + } else if (new_position.z() >= 8000.0L*cm && new_position.z() <= 8900.0L*cm) { + z_index = static_cast(std::floor((+local_position.z() - 100.0L*cm) / (scint_x_edge_length) )); + } else if (new_position.z() >= 9000.0L*cm && new_position.z() <= 9900.0L*cm) { + z_index = static_cast(std::floor((+local_position.z() - 200.0L*cm) / (scint_x_edge_length) )); + } else if (new_position.z() >= 10000.0L*cm && new_position.z() <= 10900.0L*cm) { + z_index = static_cast(std::floor((+local_position.z() - 300.0L*cm) / (scint_x_edge_length) )); + } else if (new_position.z() >= 11000.0L*cm && new_position.z() <= 11900.0L*cm) { + z_index = static_cast(std::floor((+local_position.z() - 400.0L*cm) / (scint_x_edge_length) )); + } else if (new_position.z() >= 12000.0L*cm && new_position.z() <= 12900.0L*cm) { + z_index = static_cast(std::floor((+local_position.z() - 500.0L*cm) / (scint_x_edge_length) )); + } else if (new_position.z() >= 13000.0L*cm && new_position.z() <= 13900.0L*cm) { + z_index = static_cast(std::floor((+local_position.z() - 600.0L*cm) / (scint_x_edge_length) )); + } else if (new_position.z() >= 14000.0L*cm && new_position.z() <= 14900.0L*cm) { + z_index = static_cast(std::floor((+local_position.z() - 700.0L*cm) / (scint_x_edge_length) )); + } else if (new_position.z() >= 15000.0L*cm && new_position.z() <= 15900.0L*cm) { + z_index = static_cast(std::floor((+local_position.z() - 800.0L*cm) / (scint_x_edge_length) )); + } else if (new_position.z() >= 16000.0L*cm && new_position.z() <= 16900.0L*cm) { + z_index = static_cast(std::floor((+local_position.z() - 900.0L*cm) / (scint_x_edge_length) )); + } + + } else if (_rotation == 1){ + + if (new_position.x() <= -4050.0L*cm) { + x_index = static_cast(std::floor(+local_position.x() / (scint_x_edge_length) )); + } else if (new_position.x() >= -3950.0L*cm && new_position.x() <= -3050.0L*cm) { + x_index = static_cast(std::floor((+local_position.x() - 100.0L*cm) / (scint_x_edge_length) )); + } else if (new_position.x() >= -2950.0L*cm && new_position.x() <= -2050.0L*cm) { + x_index = static_cast(std::floor((+local_position.x() - 200.0L*cm) / (scint_x_edge_length) )); + } else if (new_position.x() >= -1950.0L*cm && new_position.x() <= -1050.0L*cm) { + x_index = static_cast(std::floor((+local_position.x() - 300.0L*cm) / (scint_x_edge_length) )); + } else if (new_position.x() >= -950.0L*cm && new_position.x() <= -50.0L*cm) { + x_index = static_cast(std::floor((+local_position.x() - 400.0L*cm) / (scint_x_edge_length) )); + } else if (new_position.x() >= 50.0L*cm && new_position.x() <= 950.0L*cm) { + x_index = static_cast(std::floor((+local_position.x() - 500.0L*cm) / (scint_x_edge_length) )); + } else if (new_position.x() >= 1050.0L*cm && new_position.x() <= 1950.0L*cm) { + x_index = static_cast(std::floor((+local_position.x() - 600.0L*cm) / (scint_x_edge_length) )); + } else if (new_position.x() >= 2050.0L*cm && new_position.x() <= 2950.0L*cm) { + x_index = static_cast(std::floor((+local_position.x() - 700.0L*cm) / (scint_x_edge_length) )); + } else if (new_position.x() >= 3050.0L*cm && new_position.x() <= 3950.0L*cm) { + x_index = static_cast(std::floor((+local_position.x() - 800.0L*cm) / (scint_x_edge_length) )); + } else if (new_position.x() >= 4050.0L*cm && new_position.x() <= 4950.0L*cm) { + x_index = static_cast(std::floor((+local_position.x() - 900.0L*cm) / (scint_x_edge_length) )); + } + if (new_position.z() <= 7900.0L*cm) { + z_index = static_cast(std::floor(+local_position.z() / (scint_y_edge_length) )); + } else if (new_position.z() >= 8000.0L*cm && new_position.z() <= 8900.0L*cm) { + z_index = static_cast(std::floor((+local_position.z() - 100.0L*cm) / (scint_y_edge_length) )); + } else if (new_position.z() >= 9000.0L*cm && new_position.z() <= 9900.0L*cm) { + z_index = static_cast(std::floor((+local_position.z() - 200.0L*cm) / (scint_y_edge_length) )); + } else if (new_position.z() >= 10000.0L*cm && new_position.z() <= 10900.0L*cm) { + z_index = static_cast(std::floor((+local_position.z() - 300.0L*cm) / (scint_y_edge_length) )); + } else if (new_position.z() >= 11000.0L*cm && new_position.z() <= 11900.0L*cm) { + z_index = static_cast(std::floor((+local_position.z() - 400.0L*cm) / (scint_y_edge_length) )); + } else if (new_position.z() >= 12000.0L*cm && new_position.z() <= 12900.0L*cm) { + z_index = static_cast(std::floor((+local_position.z() - 500.0L*cm) / (scint_y_edge_length) )); + } else if (new_position.z() >= 13000.0L*cm && new_position.z() <= 13900.0L*cm) { + z_index = static_cast(std::floor((+local_position.z() - 600.0L*cm) / (scint_y_edge_length) )); + } else if (new_position.z() >= 14000.0L*cm && new_position.z() <= 14900.0L*cm) { + z_index = static_cast(std::floor((+local_position.z() - 700.0L*cm) / (scint_y_edge_length) )); + } else if (new_position.z() >= 15000.0L*cm && new_position.z() <= 15900.0L*cm) { + z_index = static_cast(std::floor((+local_position.z() - 800.0L*cm) / (scint_y_edge_length) )); + } else if (new_position.z() >= 16000.0L*cm && new_position.z() <= 16900.0L*cm) { + z_index = static_cast(std::floor((+local_position.z() - 900.0L*cm) / (scint_y_edge_length) )); + } + } + + const auto x_name = std::to_string(x_index); + const auto z_name = std::to_string(z_index); + + _hit_collection->insert(new Tracking::Hit( + particle, + trackID, + parentID, + std::to_string(1UL + y_index) + + (z_index < 10UL ? "000" + z_name : (z_index < 100UL ? "00" + z_name : (z_index < 1000UL ? "0" + z_name : z_name))) + + (x_index < 10UL ? "000" + x_name : (x_index < 100UL ? "00" + x_name : (x_index < 1000UL ? "0" + x_name : x_name))), + deposit / Units::Energy, + G4LorentzVector(new_position.t() / Units::Time, new_position.vect() / Units::Length), + G4LorentzVector(new_momentum.e() / Units::Energy, new_momentum.vect() / Units::Momentum))); + + return true; +} +//---------------------------------------------------------------------------------------------- + +//__Post-Event Processing_______________________________________________________________________ +void Detector::EndOfEvent(G4HCofThisEvent*) { + if (_hit_collection->GetSize() == 0) + return; + + MuonDataController* controller = MuonDataController::getMuonDataController(); + if(controller->getOn() ==true){ + if(controller->getDecayInEvent() == false){ + return; + } + if(controller->getDecayInZone() == false){ + G4cout<<"Decay In Zone is false"<ExtraDetails()); + root_data.insert(root_data.cend(), gen_particle_data.cbegin(), gen_particle_data.cend()); + root_data.insert(root_data.cend(), extra_gen_data.cbegin(), extra_gen_data.cend()); + + Analysis::ROOT::DataEntry metadata; + metadata.reserve(2UL); + metadata.push_back(collection_data[0UL].size()); + metadata.push_back(gen_particle_data[0UL].size()); + + Analysis::ROOT::FillNTuple(DataName, Detector::DataKeyTypes, metadata, root_data); + if (verboseLevel >= 2 && _hit_collection) + std::cout << *_hit_collection; +} +//---------------------------------------------------------------------------------------------- + +//Build 1 Module for detector +G4VPhysicalVolume* Detector::ConstructScintillatorLayer(G4LogicalVolume* ModuleVolume, int module_number, int layer_number, dimension module_x_displacement, dimension module_y_displacement, dimension layer_z_displacement){ + + auto current = new Scintillator("M" + std::to_string(module_number) + "L" + std::to_string(layer_number), + layer_x_edge_length, + layer_y_edge_length, + full_layer_height, + scintillator_casing_thickness); + + _scintillators.push_back(current); + + UNUSED(module_x_displacement); + UNUSED(module_y_displacement); + + return current->PlaceIn(ModuleVolume, G4Translate3D(0.0, 0.0, layer_z_displacement) ); + //G4Translate3D(0.5*layer_x_edge_length, 0.5*layer_y_edge_length, layer_z_displacement) +} + +G4VPhysicalVolume* Detector::ConstructModule(G4LogicalVolume* DetectorVolume, int tag_number, dimension detector_x, dimension detector_y, dimension detector_z){ + + auto ModuleVolume = Construction::BoxVolume("Module" + std::to_string(tag_number), module_x_edge_length + module_case_thickness, module_y_edge_length + module_case_thickness, full_module_height); + + // auto ModuleCaseVolume = Construction::OpenBoxVolume("Module" + std::to_string(tag_number) + "_Case", module_x_edge_length, module_y_edge_length, full_module_height, + // module_case_thickness, Construction::Material::Iron, *ModuleVisAttr()); + + // Construction::PlaceVolume(ModuleCaseVolume, ModuleVolume, Construction::Transform(0.0, 0.0, 0.0)); + + + for (std::size_t layer{}; layer < layer_count; ++layer) { + auto current = Detector::ConstructScintillatorLayer(ModuleVolume, tag_number, layer, + 0*m, + 0*m, + get_layer_z_displacement(layer)); + UNUSED(current); + } + + UNUSED(detector_x); + UNUSED(detector_y); + UNUSED(detector_z); + + //CONSTRUCTING AND INSERTING STEEL BEAMS + + for (int beam_layer = 0; beam_layer < NBEAMLAYERS; beam_layer++){ + auto BeamL1 = Construction::OpenBoxVolume("Module" + std::to_string(tag_number) + "BL" + std::to_string(beam_layer) + "PL1", beam_x_edge_length, beam_y_edge_length, module_beam_heights[beam_layer], + beam_thickness, Construction::Material::Iron, Construction::CasingAttributes()); + auto BeamL2 = Construction::OpenBoxVolume("Module" + std::to_string(tag_number) + "BL" + std::to_string(beam_layer) + "PL2", beam_x_edge_length, beam_y_edge_length, module_beam_heights[beam_layer], + beam_thickness, Construction::Material::Iron, Construction::CasingAttributes()); + auto BeamR1 = Construction::OpenBoxVolume("Module" + std::to_string(tag_number) + "BL" + std::to_string(beam_layer) + "PR1", beam_x_edge_length, beam_y_edge_length, module_beam_heights[beam_layer], + beam_thickness, Construction::Material::Iron, Construction::CasingAttributes()); + auto BeamR2 = Construction::OpenBoxVolume("Module" + std::to_string(tag_number) + "BL" + std::to_string(beam_layer) + "PR2", beam_x_edge_length, beam_y_edge_length, module_beam_heights[beam_layer], + beam_thickness, Construction::Material::Iron, Construction::CasingAttributes()); + + Construction::PlaceVolume(BeamL1, ModuleVolume, Construction::Transform(-0.50*module_x_edge_length + 0.50*beam_x_edge_length, + -0.50*module_y_edge_length + 0.50*beam_y_edge_length, + -1.0*module_beam_z_pos[beam_layer])); + Construction::PlaceVolume(BeamL2, ModuleVolume, Construction::Transform(-0.50*module_x_edge_length + 0.50*beam_x_edge_length, + 0.50*module_y_edge_length - 0.50*beam_y_edge_length, + -1.0*module_beam_z_pos[beam_layer])); + Construction::PlaceVolume(BeamR1, ModuleVolume, Construction::Transform(0.50*module_x_edge_length - 0.50*beam_x_edge_length, + -0.50*module_y_edge_length + 0.50*beam_y_edge_length, + -1.0* module_beam_z_pos[beam_layer])); + Construction::PlaceVolume(BeamR2, ModuleVolume, Construction::Transform(0.50*module_x_edge_length - 0.50*beam_x_edge_length, + 0.50*module_y_edge_length - 0.50*beam_y_edge_length, + -1.0* module_beam_z_pos[beam_layer])); + } + + + // if (tag_number == 0) { + // std::cout << "ABOUT TO WRITE GDML FOR MODULE" << std::endl; + // Construction::Export(ModuleVolume, folder, file2, arg4 ); + // } + + + return Construction::PlaceVolume(ModuleVolume, DetectorVolume, + Construction::Transform(get_module_x_displacement(tag_number), + get_module_y_displacement(tag_number), + half_detector_height - steel_height - 2.0*layer_w_case - 1*layer_spacing - 0.5*full_module_height, + 0.0, 0.0, 1.0, 0.0)); + + +} + +//__Build Detector______________________________________________________________________________ +G4VPhysicalVolume* Detector::Construct(G4LogicalVolume* world) { + Scintillator::Material::Define(); + _scintillators.clear(); + // pre_data->Branch("X_S", &X_POS_STEP, "X_S/D"); + // pre_data->Branch("Y_S", &Y_POS_STEP, "Y_S/D"); + // pre_data->Branch("X_H", &X_POS_HIT, "X_H/D"); + // pre_data->Branch("Y_H", &Y_POS_HIT, "Y_H/D"); + + auto DetectorVolume = Construction::BoxVolume("Box", x_edge_length + x_edge_increase, y_edge_length, full_detector_height, + Construction::Material::Air, G4VisAttributes::GetInvisible()); + + //DetectorVolume->SetVisAttributes(G4VisAttributes::Invisible); + + for (int module_number = 0; module_number < NMODULES; module_number++){ + Detector::ConstructModule(DetectorVolume, module_number, + 0.5L, //These arguments are not actually used + 0.5L, + 0.0L); + } + + auto first_hermetic_floor = new Scintillator("HF1", + x_edge_length, + y_edge_length, + full_layer_height, + scintillator_casing_thickness); + _scintillators.push_back(first_hermetic_floor); + first_hermetic_floor->PlaceIn(DetectorVolume, G4Translate3D(0.0, 0.0, half_detector_height - 0.5*layer_w_case - steel_height)); + + auto second_hermetic_floor = new Scintillator("HF2", + x_edge_length, + y_edge_length, + full_layer_height, + scintillator_casing_thickness); + _scintillators.push_back(second_hermetic_floor); + second_hermetic_floor->PlaceIn(DetectorVolume, G4Translate3D(0.0, 0.0, half_detector_height - 1.5*layer_w_case - layer_spacing - steel_height)); + + // auto third_hermetic_floor = new Scintillator("HF3", + // x_edge_length, + // y_edge_length, + // full_layer_height, + // scintillator_casing_thickness); + // _scintillators.push_back(third_hermetic_floor); + // third_hermetic_floor->PlaceIn(DetectorVolume, G4Translate3D(0.0, 0.0, half_detector_height - 2.5*layer_w_case - 2*layer_spacing - steel_height)); + + auto hermetic_wall = new Scintillator("HW1", + full_layer_height, + y_edge_length, + wall_height, + scintillator_casing_thickness); + _scintillators.push_back(hermetic_wall); + hermetic_wall->PlaceIn(DetectorVolume, G4Translate3D(-0.5L*x_edge_length - 0.5L*full_layer_height - wall_gap, 0.0, half_detector_height - 0.5L*wall_height)); + + _steel = Construction::BoxVolume("SteelPlate", + x_edge_length, y_edge_length, steel_height, + Construction::Material::Iron, + Construction::CasingAttributes()); + Construction::PlaceVolume(_steel, DetectorVolume, Construction::Transform(0.0, 0.0, half_detector_height - 0.5*steel_height)); + + // Construction::Export(DetectorVolume, folder, file, arg4 ); + + return Construction::PlaceVolume(DetectorVolume, world, + Construction::Transform(0.5L*x_edge_length + x_displacement, 0.5L*y_edge_length + y_displacement, -0.50*full_detector_height + floor_z_top + layer_spacing + 2.*layer_w_case + steel_height)); + +} + + +//---------------------------------------------------------------------------------------------- + +namespace CMS{ + + constexpr auto ____DEFINE_ME____ = 0.0*m; + + constexpr auto earth_total_depth = 4530.0L*cm + 1825.0L*cm + 3645.0L*cm; + + constexpr auto uxc55_cavern_length = 53.0*m; + constexpr auto uxc55_inner_radius = 13.250*m; + constexpr auto uxc55_outer_radius = 14.530*m; + constexpr auto IPDepth = earth_total_depth - uxc55_outer_radius; + constexpr auto _base_depth = earth_total_depth; + constexpr auto access_shaft_x = -14.00*m; + constexpr auto access_shaft_y = 0.00*m; + constexpr auto access_shaft_z = -1.0* uxc55_outer_radius; + + constexpr auto CMSSteelThickness = 1.48L*m; + constexpr auto CMSDetectorLength = 20.00L*m; + constexpr auto CMSDetectorRadius = 8.00L*m; + + constexpr auto AS_Depth = 20.50*m; + constexpr auto AS_Width = 20.50*m; + constexpr auto AS_Height = earth_total_depth - 2 * uxc55_outer_radius; + constexpr auto AS_Thickness = 1.5*m; + + + G4LogicalVolume* EarthVolume() { + using namespace Earth; + + auto earth_box = Construction::Box("", LayerWidthX(), LayerWidthY(), TotalDepth()); + + auto modified = new G4SubtractionSolid("", + earth_box, + Construction::Box("AirBox", x_edge_length + x_edge_increase, y_edge_length, air_gap), + Construction::Transform(0.5L*x_edge_length + x_displacement, + 0.5L*y_edge_length + y_displacement, + 0.5L*(air_gap-Earth::TotalDepth()) -9.50*m )); + + return Construction::Volume(modified); + } + + G4LogicalVolume* SandstoneVolume() { + using namespace Earth; + auto sandstone_box = Construction::Box("", LayerWidthX(), LayerWidthY(), SandstoneDepth()); + + return Construction::Volume(sandstone_box, Material::SiO2, Construction::BorderAttributes()); + } + + long double BaseDepth() { + return _base_depth - Earth::TotalShift(); + } + + long double TopDepth() { + return BaseDepth() - uxc55_outer_radius; + } + long double CenterDepth() { + return BaseDepth() - uxc55_outer_radius; + } + G4Translate3D IPTransform() { + return G4Translate3D(0.0, 0.0, IPDepth); + } + + G4Translate3D Access_Shaft_Transform(){ + return G4Translate3D(access_shaft_x, access_shaft_y, static_cast(access_shaft_z)); + } + + G4Translate3D Cavern_Transform(){ + return G4Translate3D(0, 0, -0.5 * Earth::TotalDepth() + IPDepth); + //* Construction::Rotate(0, 1, 0, 90*deg)) + } + + + G4LogicalVolume* CMSRingVolume() { + return Construction::Volume(Construction::Cylinder("DetectorRing", + CMSDetectorLength, CMSDetectorRadius - CMSSteelThickness, CMSDetectorRadius), + Construction::Material::Iron, + Construction::CasingAttributes()); + } + + G4LogicalVolume* CMSVolume(){ + using namespace Construction; + auto cavern = Cylinder("cavern", uxc55_cavern_length, 0.0*m, uxc55_outer_radius); + + auto access_shaft = Construction::Box("shaft", AS_Width, AS_Depth, AS_Height ); + + return Construction::Volume(new G4UnionSolid("fake_cms", + access_shaft, + cavern, + Construction::Rotate(1, 0, 0, 90*deg) + *G4Translate3D(0.0, -1.0*static_cast(uxc55_outer_radius + 0.5*AS_Height), -0.5*uxc55_cavern_length + 0.5*AS_Width) )); + + } + + //__Calculate Subtraction of Volumes____________________________________________________________ + G4LogicalVolume* _calculate_modification(const std::string& name, + G4LogicalVolume* earth_component, + const double base_depth, + const double top_depth) { + return Construction::Volume(new G4SubtractionSolid(name, + earth_component->GetSolid(), + CMSVolume()->GetSolid(), + Construction::Transform(-0.5*uxc55_cavern_length + 0.5*AS_Width, 0.0, 0.0) + *Construction::Rotate(0, 0, 1, 90*deg) + *Construction::Transform(0.0, 0.0, -0.5 * (base_depth - top_depth) + CenterDepth() - top_depth - uxc55_outer_radius - 0.5*AS_Height)), + earth_component->GetMaterial()); + } + +} // NAMESPACE CMS + + +//__Build Earth for Detector____________________________________________________________________ +G4VPhysicalVolume* Detector::ConstructEarth(G4LogicalVolume* world){ + + //EDIT EARTH TO ACCOMODATE VOLUME + using namespace Box::CavernConstruction; + using namespace Construction; + using namespace CMS; + + Earth::Material::Define(); + + auto earth = CMS::EarthVolume(); + + const auto mix_top = Earth::TotalDepth() - Earth::MixDepth(); + const auto marl_top = mix_top - Earth::MarlDepth(); + const auto sandstone_top = marl_top - Earth::SandstoneDepth(); + + Construction::PlaceVolume(CMS::_calculate_modification("modified_mix", Earth::MixVolume(), + mix_top + Earth::MixDepth(), mix_top), + earth, Earth::MixTransform()); + + Construction::PlaceVolume(CMS::_calculate_modification("modified_marl", Earth::MarlVolume(), + marl_top + Earth::MarlDepth(), marl_top), + earth, Earth::MarlTransform()); + + auto sandstone = CMS::_calculate_modification("modified_sandstone", CMS::SandstoneVolume(), + sandstone_top + Earth::SandstoneDepth(), sandstone_top); + + /////////////// UXC55 AND CMS DETECTOR CONSTRUCTION //////////////////////////////////////// + + auto UXC_55_cavern_solid = Cylinder("UXC55_outer", uxc55_cavern_length, 0.0*m, uxc55_outer_radius); + auto UXC55_outer_solid = Cylinder("UXC55_outer", uxc55_cavern_length, uxc55_inner_radius, uxc55_outer_radius); + auto UXC55_outer_logical = Volume(UXC55_outer_solid, Construction::Material::Concrete, Construction::CasingAttributes()); + auto CMS_Detector_logical = CMSRingVolume(); + auto UXC_55_air_v1 = new G4SubtractionSolid("UXC_55_air_v1", UXC_55_cavern_solid, UXC55_outer_solid); + auto UXC_55_air_v2 = new G4SubtractionSolid("UXC_55_air_v2", UXC_55_air_v1, CMS_Detector_logical->GetSolid()); + auto UXC55_air_logical = Volume("UXC55_air", UXC_55_air_v2, Construction::Material::Air, G4VisAttributes::GetInvisible()); + + Construction::PlaceVolume(UXC55_outer_logical, earth, Cavern_Transform()*Construction::Rotate(0, 1, 0, 90*deg) ); + Construction::PlaceVolume(CMS_Detector_logical, earth, Cavern_Transform()*Construction::Rotate(0, 1, 0, 90*deg) ); + Construction::PlaceVolume(UXC55_air_logical, earth, Cavern_Transform()*Construction::Rotate(0, 1, 0, 90*deg) ); + + + /////////////// ACCESS SHAFT CONSTRUCTION //////////////////////////////////////// + + //MAKE WHOLE SHAFT HERE + + auto Access_Shaft_outer_logical = OpenBoxVolume("Access_Shaft_outer", + AS_Width, + AS_Depth, + AS_Height, + AS_Thickness, + Construction::Material::Concrete, + Construction::CasingAttributes()); + + auto Access_Shaft_Air = BoxVolume("Acess_Shaft_Air", + AS_Width - 2* AS_Thickness, + AS_Depth - 2* AS_Thickness, + AS_Height - 2* AS_Thickness, + Construction::Material::Air, + G4VisAttributes::GetInvisible()); + + Construction::PlaceVolume(Access_Shaft_outer_logical, earth, Access_Shaft_Transform() ); + Construction::PlaceVolume(Access_Shaft_Air, earth, Access_Shaft_Transform()); + + auto modified = Construction::Volume(new G4SubtractionSolid("ModifiedSandstone", + sandstone->GetSolid(), + Construction::Box("AirBox", x_edge_length + x_edge_increase, y_edge_length, air_gap), + Construction::Transform(0.5L*x_edge_length + x_displacement, + 0.5L*y_edge_length + y_displacement, + 0.5L*(air_gap-Earth::SandstoneDepth()) - 9.50*m)), + Earth::Material::SiO2); + + Construction::PlaceVolume(modified, earth, Earth::SandstoneTransform()); + + //PLACE WHOLE THING IN WORLD + + //auto mod = CMS::_calculate_modification("modified_marl", Earth::MarlVolume(), + // marl_top + Earth::MarlDepth(), marl_top); + + + ////export geometry to gdml files + // Construction::Export(CMSVolume(), folder, file5, arg4 ); + // Construction::Export(earth, folder, file4, arg4 ); + + + //// Put Range Cuts on earth volume + // G4Region* cut_region = new G4Region("Earth_Cut_Region"); + // cut_region->AddRootLogicalVolume(earth); + // G4ProductionCuts* cuts = new G4ProductionCuts; + // cuts->SetProductionCut(1.0*m,G4ProductionCuts::GetIndex("gamma")); + // cuts->SetProductionCut(1.0*m,G4ProductionCuts::GetIndex("e-")); + // cuts->SetProductionCut(1.0*m,G4ProductionCuts::GetIndex("e+")); + // cuts->SetProductionCut(1.0*m,G4ProductionCuts::GetIndex("proton")); + // cut_region->SetProductionCuts(cuts); + + + + return Construction::PlaceVolume(earth, world, Earth::Transform()); + +} + + +//---------------------------------------------------------------------------------------------- + +} /* namespace Box */ ////////////////////////////////////////////////////////////////////////// + + + +} } /* namespace MATHUSLA::MU */ + diff --git a/Box.hh b/Box.hh new file mode 100644 index 0000000..8e2459a --- /dev/null +++ b/Box.hh @@ -0,0 +1,97 @@ +/* include/geometry/Box.hh + * + * Copyright 2018 Brandon Gomes + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#ifndef MU__GEOMETRY_BOX_HH +#define MU__GEOMETRY_BOX_HH +#pragma once +#include "TROOT.h" +#include "TH1.h" +#include "TTree.h" +#include "TFile.h" +#include "analysis.hh" +#include "geometry/Construction.hh" + +namespace MATHUSLA { namespace MU { + +namespace Box { //////////////////////////////////////////////////////////////////////////////// + +constexpr auto Box_IP_Depth = 85.47*m; + +class Scintillator { +public: + Scintillator(const std::string& name, + const double length, + const double height, + const double width, + const double thickness); + + struct Material { + static G4Material* Casing; + static G4Material* Scintillator; + static void Define(); + private: + Material(); + }; + + double GetLength() const { return _length; } + double GetHeight() const { return _height; } + double GetWidth() const { return _width; } + double GetCasingThickness() const { return _thickness; } + G4LogicalVolume* GetVolume() const { return _lvolume; } + G4VPhysicalVolume* GetSensitiveVolume() const { return _sensitive; } + + void Register(G4VSensitiveDetector* detector); + + G4VPhysicalVolume* PlaceIn(G4LogicalVolume* parent, + const G4Transform3D& transform); + +private: + G4LogicalVolume* _lvolume; + G4VPhysicalVolume* _sensitive; + std::string _name; + double _length, _height, _width, _thickness; +}; + +class Detector : public G4VSensitiveDetector { +public: + Detector(); + + void Initialize(G4HCofThisEvent* event); + G4bool ProcessHits(G4Step* step, G4TouchableHistory*); + void EndOfEvent(G4HCofThisEvent*); + + static const bool DataPerEvent = true; + static const std::string& DataName; + static const Analysis::ROOT::DataKeyList DataKeys; + static const Analysis::ROOT::DataKeyTypeList DataKeyTypes; + static TTree* pre_data; + + static G4VPhysicalVolume* Construct(G4LogicalVolume* world); + static G4VPhysicalVolume* ConstructEarth(G4LogicalVolume* world); + static G4VPhysicalVolume* ConstructModule(G4LogicalVolume* detector, int tag_number, double detector_x, double detector_y, double detector_z); + static G4VPhysicalVolume* ConstructScintillatorLayer(G4LogicalVolume* Module_volume, int module_number, int layer_number, double module_x_displacement, double module_y_displacement, double layer_z_displacement); + static bool SaveAll; + + static void WritePreData(); + static void SaveInfo(std::string & prefix); +}; + +} /* namespace Box */ ////////////////////////////////////////////////////////////////////////// + +} } /* namespace MATHUSLA::MU */ + +#endif /* MU__GEOMETRY_BOX_HH */ diff --git a/Construction.cc b/Construction.cc new file mode 100644 index 0000000..f8fd8e7 --- /dev/null +++ b/Construction.cc @@ -0,0 +1,667 @@ + +#include "geometry/Construction.hh" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "geometry/Box.hh" +#include "geometry/Cosmic.hh" +#include "geometry/Prototype.hh" +#include "geometry/Flat.hh" +#include "geometry/MuonMapper.hh" + +#include "MultiParticleChangeCrossSection.hh" + +#include "util/io.hh" + +namespace MATHUSLA { namespace MU { + +namespace { //////////////////////////////////////////////////////////////////////////////////// + +//__NIST Material Manager_______________________________________________________________________ +const auto _nist = G4NistManager::Instance(); +//---------------------------------------------------------------------------------------------- + +//__Detector Details for Builder________________________________________________________________ +std::string _detector; +std::string _export_dir; +bool _data_per_event; +std::string _data_name; +const Analysis::ROOT::DataKeyList* _data_keys; +const Analysis::ROOT::DataKeyTypeList* _data_key_types; +bool _save_option; +bool _cut_save_option; +//---------------------------------------------------------------------------------------------- + +//__Detector List_______________________________________________________________________________ +const std::string& _detectors = "Prototype Flat Cosmic Box MuonMapper"; +//---------------------------------------------------------------------------------------------- + +} /* anonymous namespace */ //////////////////////////////////////////////////////////////////// + +namespace Construction { /////////////////////////////////////////////////////////////////////// + +//__Construction Materials______________________________________________________________________ +G4Element* Material::H = _nist->FindOrBuildElement("H"); +G4Element* Material::C = _nist->FindOrBuildElement("C"); +G4Element* Material::N = _nist->FindOrBuildElement("N"); +G4Element* Material::O = _nist->FindOrBuildElement("O"); +G4Element* Material::F = _nist->FindOrBuildElement("F"); +G4Element* Material::Al = _nist->FindOrBuildElement("Al"); +G4Element* Material::Si = _nist->FindOrBuildElement("Si"); +G4Element* Material::S = _nist->FindOrBuildElement("S"); +G4Element* Material::Ar = _nist->FindOrBuildElement("Ar"); +G4Element* Material::Ca = _nist->FindOrBuildElement("Ca"); +G4Material* Material::Air = _nist->FindOrBuildMaterial("G4_AIR"); +G4Material* Material::Aluminum = _nist->FindOrBuildMaterial("G4_Al"); +G4Material* Material::Bakelite = _nist->FindOrBuildMaterial("G4_BAKELITE"); +G4Material* Material::Copper = _nist->FindOrBuildMaterial("G4_Cu"); +G4Material* Material::Concrete = _nist->FindOrBuildMaterial("G4_CONCRETE"); +G4Material* Material::Iron = _nist->FindOrBuildMaterial("G4_Fe"); +G4Material* Material::PolystyreneFoam = _nist->BuildMaterialWithNewDensity("PolystyreneFoam", "G4_POLYSTYRENE", 32.0*kg/m3); +G4Material* Material::Polyvinyltoluene = _nist->FindOrBuildMaterial("G4_PLASTIC_SC_VINYLTOLUENE"); +//---------------------------------------------------------------------------------------------- + +//__Detector Messenger Directory Path___________________________________________________________ +const std::string Builder::MessengerDirectory = "/det/"; +//---------------------------------------------------------------------------------------------- + +//__Builder Constructor_________________________________________________________________________ +Builder::Builder(const std::string& detector, + const std::string& export_dir, + const bool save_option, + const bool cut_save_option) + : G4VUserDetectorConstruction(), G4UImessenger(MessengerDirectory, "Particle Detectors.") { + _detector = detector; + _export_dir = export_dir; + _save_option = save_option; + _cut_save_option = cut_save_option; + + _select = CreateCommand("select", "Select Detector."); + _select->SetParameterName("detector", false); + _select->SetDefaultValue("Prototype"); + _select->SetCandidates(_detectors.c_str()); + _select->AvailableForStates(G4State_PreInit, G4State_Idle); + + _list = CreateCommand("list", "List Avaliable Detector."); + _list->AvailableForStates(G4State_PreInit, G4State_Idle); + + _current = CreateCommand("current", "Current Detector."); + _current->AvailableForStates(G4State_PreInit, G4State_Idle); +} +//---------------------------------------------------------------------------------------------- + +//__Build World and Detector Geometry___________________________________________________________ +G4VPhysicalVolume* Builder::Construct() { + G4GeometryManager::GetInstance()->OpenGeometry(); + G4PhysicalVolumeStore::GetInstance()->Clean(); + G4LogicalVolumeStore::GetInstance()->Clean(); + G4SolidStore::GetInstance()->Clean(); + + G4GeometryManager::GetInstance()->SetWorldMaximumExtent(WorldLength); + + std::cout << "Computed tolerance = " + << G4GeometryTolerance::GetInstance()->GetSurfaceTolerance() / m << " m\n"; + + + auto CosmicWorldLV = BoxVolume("CosmicWorld", WorldLength, WorldLength, WorldLength - 1552*m); + auto worldLV = BoxVolume("World", WorldLength, WorldLength, WorldLength - 700*m); + + + if (!_export_dir.empty()) { + if (_detector == "Flat") { + Export(Flat::Detector::Construct(worldLV), _export_dir, "flat.gdml"); + Export(Flat::Detector::ConstructEarth(worldLV), _export_dir, "flat.earth.gdml"); + } else if (_detector == "Cosmic") { + Export(Cosmic::Detector::Construct(CosmicWorldLV), _export_dir, "cosmic.gdml"); + Export(Cosmic::Detector::ConstructEarth(CosmicWorldLV), _export_dir, "cosmic.earth.gdml"); + } else if (_detector == "Box") { + Export(Box::Detector::Construct(worldLV), _export_dir, "box.gdml"); + Export(Box::Detector::ConstructEarth(worldLV), _export_dir, "box.earth.gdml"); + } else if (_detector == "MuonMapper") { + Export(MuonMapper::Detector::Construct(worldLV), _export_dir, "muon_mapper.gdml"); + Export(MuonMapper::Detector::ConstructEarth(worldLV), _export_dir, "muon_mapper.earth.gdml"); + } else { + Export(Prototype::Detector::Construct(worldLV), _export_dir, "prototype.gdml"); + Export(Prototype::Detector::ConstructEarth(worldLV), _export_dir, "prototype.earth.gdml"); + } + } else { + if (_detector == "Flat") { + Flat::Detector::Construct(worldLV); + Flat::Detector::ConstructEarth(worldLV); + } else if (_detector == "Cosmic") { + Cosmic::Detector::Construct(CosmicWorldLV); + Cosmic::Detector::ConstructEarth(CosmicWorldLV); + } else if (_detector == "Box") { + Box::Detector::Construct(worldLV); + Box::Detector::ConstructEarth(worldLV); + } else if (_detector == "MuonMapper") { + MuonMapper::Detector::Construct(worldLV); + MuonMapper::Detector::ConstructEarth(worldLV); + } else { + Prototype::Detector::Construct(worldLV); + Prototype::Detector::ConstructEarth(worldLV); + } + } + + Builder::SetSaveOption(_save_option, _cut_save_option); + + auto world = PlaceVolume(worldLV, nullptr); + auto cosmic_world = PlaceVolume(CosmicWorldLV, nullptr); + + if (!_export_dir.empty()) { + if (_detector == "Flat") { + Export(world, _export_dir, "world.flat.gdml"); + } else if (_detector == "Cosmic") { + Export(cosmic_world, _export_dir, "world.cosmic.gdml"); + } else if (_detector == "Box") { + Export(world, _export_dir, "world.box.gdml"); + } else if (_detector == "MuonMapper") { + Export(world, _export_dir, "world.muon_mapper.gdml"); + } else { + Export(world, _export_dir, "world.prototype.gdml"); + } + } + + std::cout << "Materials: " + << *G4Material::GetMaterialTable() << '\n'; + + // const std::string folder = "detector_geo"; + // const std::string file = "world.gdml"; + // const std::string arg4 = "http://service-spi.web.cern.ch/service-spi/app/releases/GDML/Schema/gdml.xsd"; + + if (_detector == "Cosmic") { + // Construction::Export(cosmic_world, folder, file, arg4); + return cosmic_world; + } else { + // Construction::Export(world, folder, file, arg4); + return world; + } + +} +//---------------------------------------------------------------------------------------------- + +//__Build Detector______________________________________________________________________________ +void Builder::ConstructSDandField() { + if (_detector == "Flat") { + _data_per_event = Flat::Detector::DataPerEvent; + _data_name = Flat::Detector::DataName; + _data_keys = &Flat::Detector::DataKeys; + _data_key_types = &Flat::Detector::DataKeyTypes; + G4SDManager::GetSDMpointer()->AddNewDetector(new Flat::Detector); + } else if (_detector == "Cosmic") { + _data_per_event = Cosmic::Detector::DataPerEvent; + _data_name = Cosmic::Detector::DataName; + _data_keys = &Cosmic::Detector::DataKeys; + _data_key_types = &Cosmic::Detector::DataKeyTypes; + G4SDManager::GetSDMpointer()->AddNewDetector(new Cosmic::Detector); + + G4LogicalVolume* earth_volume = G4LogicalVolumeStore::GetInstance()->GetVolume("ModifiedSandstoneCosmic"); + MultiParticleChangeCrossSection* biasingOperator = new MultiParticleChangeCrossSection; + biasingOperator->AddParticle( "mu+" ); + biasingOperator->AddParticle( "mu-" ); + biasingOperator->AttachTo(earth_volume); + G4cout << " Attaching biasing operator " << biasingOperator->GetName() + << " to logical volume " << earth_volume->GetName() + << G4endl; + } else if (_detector == "Box") { + _data_per_event = Box::Detector::DataPerEvent; + _data_name = Box::Detector::DataName; + _data_keys = &Box::Detector::DataKeys; + _data_key_types = &Box::Detector::DataKeyTypes; + G4SDManager::GetSDMpointer()->AddNewDetector(new Box::Detector); + + G4LogicalVolume* earth_volume = G4LogicalVolumeStore::GetInstance()->GetVolume("ModifiedSandstone"); + MultiParticleChangeCrossSection* biasingOperator = new MultiParticleChangeCrossSection; + biasingOperator->AddParticle( "mu+" ); + biasingOperator->AddParticle( "mu-" ); + biasingOperator->AttachTo(earth_volume); + G4cout << " Attaching biasing operator " << biasingOperator->GetName() + << " to logical volume " << earth_volume->GetName() + << G4endl; + } else if (_detector == "MuonMapper") { + _data_per_event = MuonMapper::Detector::DataPerEvent; + _data_name = MuonMapper::Detector::DataName; + _data_keys = &MuonMapper::Detector::DataKeys; + _data_key_types = &MuonMapper::Detector::DataKeyTypes; + G4SDManager::GetSDMpointer()->AddNewDetector(new MuonMapper::Detector); + } else { + _data_per_event = Prototype::Detector::DataPerEvent; + _data_name = Prototype::Detector::DataName; + _data_keys = &Prototype::Detector::DataKeys; + _data_key_types = &Prototype::Detector::DataKeyTypes; + G4SDManager::GetSDMpointer()->AddNewDetector(new Prototype::Detector); + } +} +//---------------------------------------------------------------------------------------------- + +//__Builder Messenger Set New Value_____________________________________________________________ +void Builder::SetNewValue(G4UIcommand* command, G4String value) { + if (command == _select && value != _detector.c_str()) { + SetDetector(value); + } else if (command == _list) { + std::cout << "Detectors: " << _detectors << "\n"; + } else if (command == _current) { + std::cout << "Current Detector: " << _detector << "\n"; + } +} +//---------------------------------------------------------------------------------------------- + +//__Set Current Detector________________________________________________________________________ +void Builder::SetDetector(const std::string& detector) { + _detector = detector; + Command::Execute("/run/reinitializeGeometry", + "/run/geometryModified", + "/run/initialize", + "/vis/viewer/clearTransients"); +} +//---------------------------------------------------------------------------------------------- + +//__Set Current Detector Save Option____________________________________________________________ +void Builder::SetSaveOption(bool save_option, bool cut_save_option) { + if (_detector == "Flat") { + Flat::Detector::SaveAll = save_option; + } else if (_detector == "Box") { + Box::Detector::SaveAll = save_option; + } else if (_detector == "Cosmic") { + Cosmic::Detector::SaveAll = save_option; + Cosmic::Detector::SaveCut = cut_save_option; + } else if (_detector == "MuonMapper") { + MuonMapper::Detector::SaveAll = save_option; + } else { + Prototype::Detector::SaveAll = save_option; + } +} +//---------------------------------------------------------------------------------------------- + +//__Save some information in a txt file____________________________________________________________ +void Builder::SaveInfo(std::string & prefix){ + if (_detector == "Flat") {;} + else if (_detector == "Box") { + Box::Detector::SaveInfo(prefix); + } +} +//---------------------------------------------------------------------------------------------- + + +//__Get Current Detector Name___________________________________________________________________ +const std::string& Builder::GetDetectorName() { + return _detector; + +} +//---------------------------------------------------------------------------------------------- + +//__Get Current Detector Data Preference________________________________________________________ +bool Builder::IsDetectorDataPerEvent() { + return _data_per_event; +} +//---------------------------------------------------------------------------------------------- + +//__Get Current Detector Data Prefix____________________________________________________________ +const std::string& Builder::GetDetectorDataName() { + return _data_name; +} +//---------------------------------------------------------------------------------------------- + +//__Get Current Detector Data Keys______________________________________________________________ +const Analysis::ROOT::DataKeyList& Builder::GetDetectorDataKeys() { + return *_data_keys; +} +//---------------------------------------------------------------------------------------------- + +//__Get Current Detector Data Key Types_________________________________________________________ +const Analysis::ROOT::DataKeyTypeList& Builder::GetDetectorDataKeyTypes() { + return *_data_key_types; +} +//---------------------------------------------------------------------------------------------- + +//////////////////////////////////////////////////////////////////////////////////////////////// + +//__Sensitive Material Attribute Definition_____________________________________________________ +const G4VisAttributes SensitiveAttributes() { + auto attr = G4VisAttributes(G4Colour(0., 1., 0., 1.0)); + attr.SetForceSolid(true); + return attr; +} + +const G4VisAttributes SensitiveAttributes2() { + auto attr = G4VisAttributes(G4Colour(1., 0., 0., 1.0)); + attr.SetForceSolid(true); + return attr; +} + +const G4VisAttributes HighlightRed() { + auto attr = G4VisAttributes(G4Colour(1., 0., 0.)); + attr.SetForceSolid(true); + return attr; +} + + +//---------------------------------------------------------------------------------------------- + +//__Iron Plating Attributes +// +const G4VisAttributes IronAttributes() { + auto attr = G4VisAttributes(G4Colour(0.5, 0.5, 0.5)); + attr.SetForceSolid(true); + return attr; +} + +const G4VisAttributes AlAttributes() { + auto attr = G4VisAttributes(G4Colour(0.4, 0.6, 0.7)); + attr.SetForceSolid(true); + return attr; +} + + +//__Casing Material Attribute Definition________________________________________________________ +const G4VisAttributes CasingAttributes() { + auto attr = G4VisAttributes(G4Colour(0., 0., 1., 0.2)); + attr.SetForceWireframe(true); + return attr; +} +//---------------------------------------------------------------------------------------------- + +//__Border Attribute Definition_________________________________________________________________ +const G4VisAttributes BorderAttributes() { + return G4VisAttributes(false); +} +//---------------------------------------------------------------------------------------------- + +//__Box Builder_________________________________________________________________________________ +G4Box* Box(const std::string& name, + const double width, + const double height, + const double depth) { + return new G4Box(name, 0.5 * width, 0.5 * height, 0.5 * depth); +} +//---------------------------------------------------------------------------------------------- + +//__Trapezoid Builder___________________________________________________________________________ +G4Trap* Trap(const std::string& name, + const double height, + const double minwidth, + const double maxwidth, + const double depth) { + return new G4Trap(name, + 0.5 * height, 0, 0, + 0.5 * depth, 0.5 * minwidth, 0.5 * minwidth, 0, + 0.5 * depth, 0.5 * maxwidth, 0.5 * maxwidth, 0); +} +//---------------------------------------------------------------------------------------------- + +//__Cylinder Builder____________________________________________________________________________ +G4Tubs* Cylinder(const std::string& name, + const double height, + const double inner_radius, + const double outer_radius) { + return new G4Tubs(name, inner_radius, outer_radius, 0.5 * height, 0, 360*deg); +} +//---------------------------------------------------------------------------------------------- + +//__Volume Builder______________________________________________________________________________ +G4LogicalVolume* Volume(const std::string& name, + G4VSolid* solid, + G4Material* material, + const G4VisAttributes& attr) { + auto out = new G4LogicalVolume(solid, material, name); + out->SetVisAttributes(attr); + return out; +} +//---------------------------------------------------------------------------------------------- + +//__Volume Builder______________________________________________________________________________ +G4LogicalVolume* Volume(G4VSolid* solid, + G4Material* material, + const G4VisAttributes& attr) { + return Volume(solid->GetName(), solid, material, attr); +} +//---------------------------------------------------------------------------------------------- + +//__Volume Builder______________________________________________________________________________ +G4LogicalVolume* Volume(const std::string& name, + G4VSolid* solid, + const G4VisAttributes& attr) { + return Volume(name, solid, Material::Air, attr); +} +//---------------------------------------------------------------------------------------------- + +//__Volume Builder______________________________________________________________________________ +G4LogicalVolume* Volume(G4VSolid* solid, + const G4VisAttributes& attr) { + return Volume(solid, Material::Air, attr); +} +//---------------------------------------------------------------------------------------------- + +//__Box Volume Builder__________________________________________________________________________ +G4LogicalVolume* BoxVolume(const std::string& name, + const double width, + const double height, + const double depth, + G4Material* material, + const G4VisAttributes& attr) { + return Volume( + Box(name, width, height, depth), + material, attr); +} +//---------------------------------------------------------------------------------------------- + +//__Open Box Volume Builder_____________________________________________________________________ +G4LogicalVolume* OpenBoxVolume(const std::string& name, + const double width, + const double height, + const double depth, + const double thickness, + G4Material* material, + const G4VisAttributes& attr) { + auto outer = Box(name, + width, + height, + depth); + auto inner = Box(name, + width - 2 * thickness, + height - 2 * thickness, + depth - 2 * thickness); + return Volume(new G4SubtractionSolid(name, outer, inner), material, attr); +} +//---------------------------------------------------------------------------------------------- + +//__Physical Volume Placer______________________________________________________________________ +G4VPhysicalVolume* PlaceVolume(const std::string& name, + G4LogicalVolume* current, + G4LogicalVolume* parent, + const G4Transform3D& transform) { + return new G4PVPlacement(transform, current, name, parent, false, 0); +} +//---------------------------------------------------------------------------------------------- + +//__Physical Volume Placer______________________________________________________________________ +G4VPhysicalVolume* PlaceVolume(G4LogicalVolume* current, + G4LogicalVolume* parent, + const G4Transform3D& transform) { + return PlaceVolume(current->GetName(), current, parent, transform); +} +//---------------------------------------------------------------------------------------------- + +//__Physical Volume Placer______________________________________________________________________ +G4VPhysicalVolume* PlaceVolume(const std::string& name, + G4VSolid* solid, + G4Material* material, + G4LogicalVolume* parent, + const G4Transform3D& transform) { + return PlaceVolume(name, Volume(name, solid, material), parent, transform); +} +//---------------------------------------------------------------------------------------------- + +//__Physical Volume Placer______________________________________________________________________ +G4VPhysicalVolume* PlaceVolume(G4VSolid* solid, + G4Material* material, + G4LogicalVolume* parent, + const G4Transform3D& transform) { + return PlaceVolume(solid->GetName(), solid, material, parent, transform); +} +//---------------------------------------------------------------------------------------------- + +//__Physical Volume Placer______________________________________________________________________ +G4VPhysicalVolume* PlaceVolume(const std::string& name, + G4LogicalVolume* current, + const G4VisAttributes& attr, + G4LogicalVolume* parent, + const G4Transform3D& transform) { + current->SetVisAttributes(attr); + return PlaceVolume(name, current, parent, transform); +} +//---------------------------------------------------------------------------------------------- + +//__Physical Volume Placer______________________________________________________________________ +G4VPhysicalVolume* PlaceVolume(G4LogicalVolume* current, + const G4VisAttributes& attr, + G4LogicalVolume* parent, + const G4Transform3D& transform) { + current->SetVisAttributes(attr); + return PlaceVolume(current, parent, transform); +} +//---------------------------------------------------------------------------------------------- + +//__Physical Volume Placer______________________________________________________________________ +G4VPhysicalVolume* PlaceVolume(const std::string& name, + G4VSolid* solid, + G4Material* material, + const G4VisAttributes& attr, + G4LogicalVolume* parent, + const G4Transform3D& transform) { + return PlaceVolume(name, Volume(solid, material, attr), parent, transform); +} +//---------------------------------------------------------------------------------------------- + +//__Physical Volume Placer______________________________________________________________________ +G4VPhysicalVolume* PlaceVolume(G4VSolid* solid, + G4Material* material, + const G4VisAttributes& attr, + G4LogicalVolume* parent, + const G4Transform3D& transform) { + return PlaceVolume(Volume(solid, material, attr), parent, transform); +} +//---------------------------------------------------------------------------------------------- + +//__Translation Transformation Generator_______________________________________________________ +G4Transform3D Transform(const double x, + const double y, + const double z) { + return G4Transform3D(G4RotationMatrix(), G4ThreeVector(x, y, z)); +} +//---------------------------------------------------------------------------------------------- + +//__Rotation/Translation Transformation Generator_______________________________________________ +G4Transform3D Transform(const G4ThreeVector& translate, + const G4ThreeVector& axis, + const double angle) { + return G4Transform3D(G4RotationMatrix(axis, angle), translate); +} +//---------------------------------------------------------------------------------------------- + +//__Rotation/Translation Transformation Generator_______________________________________________ +G4Transform3D Transform(const double x, + const double y, + const double z, + const double axisx, + const double axisy, + const double axisz, + const double angle) { + return Transform( + G4ThreeVector(x, y, z), G4ThreeVector(axisx, axisy, axisz), angle); +} +//---------------------------------------------------------------------------------------------- + +//__Rotation Transformation Generator___________________________________________________________ +G4Transform3D Rotate(const double axisx, + const double axisy, + const double axisz, + const double angle) { + return Transform(0, 0, 0, axisx, axisy, axisz, angle); +} +//---------------------------------------------------------------------------------------------- + +//__Matrix Transformation Generator_____________________________________________________________ +G4RotationMatrix Matrix(const double th1, + const double phi1, + const double th2, + const double phi2, + const double th3, + const double phi3) { + const double sinth1 = std::sin(th1); + const double sinth2 = std::sin(th2); + const double sinth3 = std::sin(th3); + auto matrix = G4RotationMatrix(); + matrix.rotateAxes( + G4ThreeVector(sinth1*std::cos(phi1), sinth1*std::sin(phi1), std::cos(th1)), + G4ThreeVector(sinth2*std::cos(phi2), sinth2*std::sin(phi2), std::cos(th2)), + G4ThreeVector(sinth3*std::cos(phi3), sinth3*std::sin(phi3), std::cos(th3))); + + if (matrix != G4RotationMatrix()) matrix.invert(); + return matrix; +} +//---------------------------------------------------------------------------------------------- + +//__Matrix Transformation Generator_____________________________________________________________ +G4RotationMatrix Matrix(const double mxx, + const double mxy, + const double mxz, + const double myx, + const double myy, + const double myz, + const double mzx, + const double mzy, + const double mzz) { + return G4RotationMatrix( + G4ThreeVector(mxx, myx, mzx), + G4ThreeVector(mxy, myy, mzy), + G4ThreeVector(mxz, myz, mzz)); +} +//---------------------------------------------------------------------------------------------- + +//__GDML File Export____________________________________________________________________________ +void Export(const G4LogicalVolume* volume, + const std::string& dir, + const std::string& file, + const std::string& schema) { + util::io::create_directory(dir); + auto path = dir + "/" + file; + if (util::io::path_exists(path)) + util::io::remove_file(path); + + static G4ThreadLocal G4GDMLParser _parser; + _parser.Write(path, volume, true, + !schema.empty() ? schema : G4GDML_DEFAULT_SCHEMALOCATION); +} +//---------------------------------------------------------------------------------------------- + +//__GDML File Export____________________________________________________________________________ +void Export(const G4VPhysicalVolume* volume, + const std::string& dir, + const std::string& file, + const std::string& schema) { + util::io::create_directory(dir); + auto path = dir + "/" + file; + if (util::io::path_exists(path)) + util::io::remove_file(path); + + static G4ThreadLocal G4GDMLParser _parser; + _parser.Write(path, volume, true, + !schema.empty() ? schema : G4GDML_DEFAULT_SCHEMALOCATION); +} +//---------------------------------------------------------------------------------------------- + +} /* namespace Construction */ ///////////////////////////////////////////////////////////////// + +} } /* namespace MATHUSLA::MU */ diff --git a/Construction.hh b/Construction.hh new file mode 100644 index 0000000..58a9d3b --- /dev/null +++ b/Construction.hh @@ -0,0 +1,294 @@ +#ifndef MU__GEOMETRY_CONSTRUCTION_HH +#define MU__GEOMETRY_CONSTRUCTION_HH +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "analysis.hh" +#include "ui.hh" + +namespace MATHUSLA { namespace MU { + +namespace Construction { /////////////////////////////////////////////////////////////////////// + +namespace Material { /////////////////////////////////////////////////////////////////////////// +extern G4Element* H; +extern G4Element* C; +extern G4Element* N; +extern G4Element* O; +extern G4Element* F; +extern G4Element* Al; +extern G4Element* Si; +extern G4Element* S; +extern G4Element* Ar; +extern G4Element* Ca; +extern G4Material* Air; +extern G4Material* Aluminum; +extern G4Material* Bakelite; +extern G4Material* Copper; +extern G4Material* Concrete; +extern G4Material* Iron; +extern G4Material* PolystyreneFoam; +extern G4Material* Polyvinyltoluene; +} /* namespace Material */ ///////////////////////////////////////////////////////////////////// + +//__Size of The World___________________________________________________________________________ +static constexpr const auto WorldLength = 1600*m; +//---------------------------------------------------------------------------------------------- + +//__Geometry Builder Class______________________________________________________________________ +class Builder : public G4VUserDetectorConstruction, public G4UImessenger { +public: + Builder(const std::string& detector, + const std::string& export_dir, + const bool save_option, + const bool cut_save_option); + G4VPhysicalVolume* Construct(); + void ConstructSDandField(); + + void SetNewValue(G4UIcommand* command, G4String value); + + static const std::string MessengerDirectory; + + static void SetDetector(const std::string& detector); + static void SetSaveOption(const bool save_option, const bool cut_save_option); + static void SaveInfo(std::string& prefix); + + static const std::string& GetDetectorName(); + static bool IsDetectorDataPerEvent(); + static const std::string& GetDetectorDataName(); + static const Analysis::ROOT::DataKeyList& GetDetectorDataKeys(); + static const Analysis::ROOT::DataKeyTypeList& GetDetectorDataKeyTypes(); + +private: + Command::NoArg* _list; + Command::NoArg* _current; + Command::StringArg* _select; +}; +//---------------------------------------------------------------------------------------------- + +//__Sensitive Material Attribute Definition_____________________________________________________ +const G4VisAttributes SensitiveAttributes(); +const G4VisAttributes SensitiveAttributes2(); +const G4VisAttributes HighlightRed(); + +//---------------------------------------------------------------------------------------------- + +//__Casing Material Attribute Definition________________________________________________________ +const G4VisAttributes CasingAttributes(); +//---------------------------------------------------------------------------------------------- + +//__Border Attribute Definition_________________________________________________________________ +const G4VisAttributes BorderAttributes(); +//---------------------------------------------------------------------------------------------- +//__Iron and Aluminum Material Attribute Defintions + +const G4VisAttributes AlAttributes(); +//---------------------------------------------------------------------------------------------- +const G4VisAttributes IronAttributes(); + +//---------------------------------------------------------------------------------------------- +//__Box Builder_________________________________________________________________________________ +G4Box* Box(const std::string& name, + const double width, + const double height, + const double depth); +//---------------------------------------------------------------------------------------------- + +//__Trapezoid Builder___________________________________________________________________________ +G4Trap* Trap(const std::string& name, + const double height, + const double minwidth, + const double maxwidth, + const double depth); +//---------------------------------------------------------------------------------------------- + +//__Cylinder Builder____________________________________________________________________________ +G4Tubs* Cylinder(const std::string& name, + const double height, + const double inner_radius, + const double outer_radius); +//---------------------------------------------------------------------------------------------- + +//__Volume Builder______________________________________________________________________________ +G4LogicalVolume* Volume(const std::string& name, + G4VSolid* solid, + G4Material* material=Material::Air, + const G4VisAttributes& attr=G4VisAttributes()); +//---------------------------------------------------------------------------------------------- + +//__Volume Builder______________________________________________________________________________ +G4LogicalVolume* Volume(G4VSolid* solid, + G4Material* material=Material::Air, + const G4VisAttributes& attr=G4VisAttributes()); +//---------------------------------------------------------------------------------------------- + +//__Volume Builder______________________________________________________________________________ +G4LogicalVolume* Volume(const std::string& name, + G4VSolid* solid, + const G4VisAttributes& attr); +//---------------------------------------------------------------------------------------------- + +//__Volume Builder______________________________________________________________________________ +G4LogicalVolume* Volume(G4VSolid* solid, + const G4VisAttributes& attr); +//---------------------------------------------------------------------------------------------- + +//__Box Volume Builder__________________________________________________________________________ +G4LogicalVolume* BoxVolume(const std::string& name, + const double width, + const double height, + const double depth, + G4Material* material=Material::Air, + const G4VisAttributes& attr=G4VisAttributes()); +//---------------------------------------------------------------------------------------------- + +//__Open Box Volume Builder_____________________________________________________________________ +G4LogicalVolume* OpenBoxVolume(const std::string& name, + const double width, + const double height, + const double depth, + const double thickness, + G4Material* material=Material::Air, + const G4VisAttributes& attr=G4VisAttributes()); +//---------------------------------------------------------------------------------------------- + +//__Physical Volume Placer______________________________________________________________________ +G4VPhysicalVolume* PlaceVolume(const std::string& name, + G4LogicalVolume* current, + G4LogicalVolume* parent, + const G4Transform3D& transform=G4Transform3D()); +//---------------------------------------------------------------------------------------------- + +//__Physical Volume Placer______________________________________________________________________ +G4VPhysicalVolume* PlaceVolume(G4LogicalVolume* current, + G4LogicalVolume* parent, + const G4Transform3D& transform=G4Transform3D()); +//---------------------------------------------------------------------------------------------- + +//__Physical Volume Placer______________________________________________________________________ +G4VPhysicalVolume* PlaceVolume(const std::string& name, + G4VSolid* solid, + G4Material* material, + G4LogicalVolume* parent, + const G4Transform3D& transform=G4Transform3D()); +//---------------------------------------------------------------------------------------------- + +//__Physical Volume Placer______________________________________________________________________ +G4VPhysicalVolume* PlaceVolume(G4VSolid* solid, + G4Material* material, + G4LogicalVolume* parent, + const G4Transform3D& transform=G4Transform3D()); +//---------------------------------------------------------------------------------------------- + +//__Physical Volume Placer______________________________________________________________________ +G4VPhysicalVolume* PlaceVolume(const std::string& name, + G4LogicalVolume* current, + const G4VisAttributes& attr, + G4LogicalVolume* parent, + const G4Transform3D& transform=G4Transform3D()); +//---------------------------------------------------------------------------------------------- + +//__Physical Volume Placer______________________________________________________________________ +G4VPhysicalVolume* PlaceVolume(G4LogicalVolume* current, + const G4VisAttributes& attr, + G4LogicalVolume* parent, + const G4Transform3D& transform=G4Transform3D()); +//---------------------------------------------------------------------------------------------- + +//__Physical Volume Placer______________________________________________________________________ +G4VPhysicalVolume* PlaceVolume(const std::string& name, + G4VSolid* solid, + G4Material* material, + const G4VisAttributes& attr, + G4LogicalVolume* parent, + const G4Transform3D& transform=G4Transform3D()); +//---------------------------------------------------------------------------------------------- + +//__Physical Volume Placer______________________________________________________________________ +G4VPhysicalVolume* PlaceVolume(G4VSolid* solid, + G4Material* material, + const G4VisAttributes& attr, + G4LogicalVolume* parent, + const G4Transform3D& transform=G4Transform3D()); +//---------------------------------------------------------------------------------------------- + +//__Translation Transformation Generator_______________________________________________________ +G4Transform3D Transform(const double x, + const double y, + const double z); +//---------------------------------------------------------------------------------------------- + +//__Rotation/Translation Transformation Generator_______________________________________________ +G4Transform3D Transform(const G4ThreeVector& translate, + const G4ThreeVector& axis, + const double angle); +//---------------------------------------------------------------------------------------------- + +//__Rotation/Translation Transformation Generator_______________________________________________ +G4Transform3D Transform(const double x, + const double y, + const double z, + const double axisx, + const double axisy, + const double axisz, + const double angle); +//---------------------------------------------------------------------------------------------- + +//__Rotation Transformation Generator___________________________________________________________ +G4Transform3D Rotate(const double axisx, + const double axisy, + const double axisz, + const double angle); +//---------------------------------------------------------------------------------------------- + +//__Matrix Transformation Generator_____________________________________________________________ +G4RotationMatrix Matrix(const double th1, + const double phi1, + const double th2, + const double phi2, + const double th3, + const double phi3); +//---------------------------------------------------------------------------------------------- + +//__Matrix Transformation Generator_____________________________________________________________ +G4RotationMatrix Matrix(const double mxx, + const double mxy, + const double mxz, + const double myx, + const double myy, + const double myz, + const double mzx, + const double mzy, + const double mzz); +//---------------------------------------------------------------------------------------------- + +//__GDML File Export____________________________________________________________________________ +void Export(const G4LogicalVolume* volume, + const std::string& dir, + const std::string& file, + const std::string& schema=""); +//---------------------------------------------------------------------------------------------- + +//__GDML File Export____________________________________________________________________________ +void Export(const G4VPhysicalVolume* volume, + const std::string& dir, + const std::string& file, + const std::string& schema=""); +//---------------------------------------------------------------------------------------------- + +} /* namespace Construction */ ///////////////////////////////////////////////////////////////// + +} } /* namespace MATHUSLA::MU */ + +#endif /* MU__GEOMETRY_CONSTRUCTION_HH */ From 942ec27d4631f63aa258274fc7b688609008fafa Mon Sep 17 00:00:00 2001 From: Abdulrahman Mohamed <88555369+Abdulrahma-M@users.noreply.github.com> Date: Mon, 8 May 2023 20:37:12 +0300 Subject: [PATCH 08/12] Delete Construction.cc --- Construction.cc | 667 ------------------------------------------------ 1 file changed, 667 deletions(-) delete mode 100644 Construction.cc diff --git a/Construction.cc b/Construction.cc deleted file mode 100644 index f8fd8e7..0000000 --- a/Construction.cc +++ /dev/null @@ -1,667 +0,0 @@ - -#include "geometry/Construction.hh" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "geometry/Box.hh" -#include "geometry/Cosmic.hh" -#include "geometry/Prototype.hh" -#include "geometry/Flat.hh" -#include "geometry/MuonMapper.hh" - -#include "MultiParticleChangeCrossSection.hh" - -#include "util/io.hh" - -namespace MATHUSLA { namespace MU { - -namespace { //////////////////////////////////////////////////////////////////////////////////// - -//__NIST Material Manager_______________________________________________________________________ -const auto _nist = G4NistManager::Instance(); -//---------------------------------------------------------------------------------------------- - -//__Detector Details for Builder________________________________________________________________ -std::string _detector; -std::string _export_dir; -bool _data_per_event; -std::string _data_name; -const Analysis::ROOT::DataKeyList* _data_keys; -const Analysis::ROOT::DataKeyTypeList* _data_key_types; -bool _save_option; -bool _cut_save_option; -//---------------------------------------------------------------------------------------------- - -//__Detector List_______________________________________________________________________________ -const std::string& _detectors = "Prototype Flat Cosmic Box MuonMapper"; -//---------------------------------------------------------------------------------------------- - -} /* anonymous namespace */ //////////////////////////////////////////////////////////////////// - -namespace Construction { /////////////////////////////////////////////////////////////////////// - -//__Construction Materials______________________________________________________________________ -G4Element* Material::H = _nist->FindOrBuildElement("H"); -G4Element* Material::C = _nist->FindOrBuildElement("C"); -G4Element* Material::N = _nist->FindOrBuildElement("N"); -G4Element* Material::O = _nist->FindOrBuildElement("O"); -G4Element* Material::F = _nist->FindOrBuildElement("F"); -G4Element* Material::Al = _nist->FindOrBuildElement("Al"); -G4Element* Material::Si = _nist->FindOrBuildElement("Si"); -G4Element* Material::S = _nist->FindOrBuildElement("S"); -G4Element* Material::Ar = _nist->FindOrBuildElement("Ar"); -G4Element* Material::Ca = _nist->FindOrBuildElement("Ca"); -G4Material* Material::Air = _nist->FindOrBuildMaterial("G4_AIR"); -G4Material* Material::Aluminum = _nist->FindOrBuildMaterial("G4_Al"); -G4Material* Material::Bakelite = _nist->FindOrBuildMaterial("G4_BAKELITE"); -G4Material* Material::Copper = _nist->FindOrBuildMaterial("G4_Cu"); -G4Material* Material::Concrete = _nist->FindOrBuildMaterial("G4_CONCRETE"); -G4Material* Material::Iron = _nist->FindOrBuildMaterial("G4_Fe"); -G4Material* Material::PolystyreneFoam = _nist->BuildMaterialWithNewDensity("PolystyreneFoam", "G4_POLYSTYRENE", 32.0*kg/m3); -G4Material* Material::Polyvinyltoluene = _nist->FindOrBuildMaterial("G4_PLASTIC_SC_VINYLTOLUENE"); -//---------------------------------------------------------------------------------------------- - -//__Detector Messenger Directory Path___________________________________________________________ -const std::string Builder::MessengerDirectory = "/det/"; -//---------------------------------------------------------------------------------------------- - -//__Builder Constructor_________________________________________________________________________ -Builder::Builder(const std::string& detector, - const std::string& export_dir, - const bool save_option, - const bool cut_save_option) - : G4VUserDetectorConstruction(), G4UImessenger(MessengerDirectory, "Particle Detectors.") { - _detector = detector; - _export_dir = export_dir; - _save_option = save_option; - _cut_save_option = cut_save_option; - - _select = CreateCommand("select", "Select Detector."); - _select->SetParameterName("detector", false); - _select->SetDefaultValue("Prototype"); - _select->SetCandidates(_detectors.c_str()); - _select->AvailableForStates(G4State_PreInit, G4State_Idle); - - _list = CreateCommand("list", "List Avaliable Detector."); - _list->AvailableForStates(G4State_PreInit, G4State_Idle); - - _current = CreateCommand("current", "Current Detector."); - _current->AvailableForStates(G4State_PreInit, G4State_Idle); -} -//---------------------------------------------------------------------------------------------- - -//__Build World and Detector Geometry___________________________________________________________ -G4VPhysicalVolume* Builder::Construct() { - G4GeometryManager::GetInstance()->OpenGeometry(); - G4PhysicalVolumeStore::GetInstance()->Clean(); - G4LogicalVolumeStore::GetInstance()->Clean(); - G4SolidStore::GetInstance()->Clean(); - - G4GeometryManager::GetInstance()->SetWorldMaximumExtent(WorldLength); - - std::cout << "Computed tolerance = " - << G4GeometryTolerance::GetInstance()->GetSurfaceTolerance() / m << " m\n"; - - - auto CosmicWorldLV = BoxVolume("CosmicWorld", WorldLength, WorldLength, WorldLength - 1552*m); - auto worldLV = BoxVolume("World", WorldLength, WorldLength, WorldLength - 700*m); - - - if (!_export_dir.empty()) { - if (_detector == "Flat") { - Export(Flat::Detector::Construct(worldLV), _export_dir, "flat.gdml"); - Export(Flat::Detector::ConstructEarth(worldLV), _export_dir, "flat.earth.gdml"); - } else if (_detector == "Cosmic") { - Export(Cosmic::Detector::Construct(CosmicWorldLV), _export_dir, "cosmic.gdml"); - Export(Cosmic::Detector::ConstructEarth(CosmicWorldLV), _export_dir, "cosmic.earth.gdml"); - } else if (_detector == "Box") { - Export(Box::Detector::Construct(worldLV), _export_dir, "box.gdml"); - Export(Box::Detector::ConstructEarth(worldLV), _export_dir, "box.earth.gdml"); - } else if (_detector == "MuonMapper") { - Export(MuonMapper::Detector::Construct(worldLV), _export_dir, "muon_mapper.gdml"); - Export(MuonMapper::Detector::ConstructEarth(worldLV), _export_dir, "muon_mapper.earth.gdml"); - } else { - Export(Prototype::Detector::Construct(worldLV), _export_dir, "prototype.gdml"); - Export(Prototype::Detector::ConstructEarth(worldLV), _export_dir, "prototype.earth.gdml"); - } - } else { - if (_detector == "Flat") { - Flat::Detector::Construct(worldLV); - Flat::Detector::ConstructEarth(worldLV); - } else if (_detector == "Cosmic") { - Cosmic::Detector::Construct(CosmicWorldLV); - Cosmic::Detector::ConstructEarth(CosmicWorldLV); - } else if (_detector == "Box") { - Box::Detector::Construct(worldLV); - Box::Detector::ConstructEarth(worldLV); - } else if (_detector == "MuonMapper") { - MuonMapper::Detector::Construct(worldLV); - MuonMapper::Detector::ConstructEarth(worldLV); - } else { - Prototype::Detector::Construct(worldLV); - Prototype::Detector::ConstructEarth(worldLV); - } - } - - Builder::SetSaveOption(_save_option, _cut_save_option); - - auto world = PlaceVolume(worldLV, nullptr); - auto cosmic_world = PlaceVolume(CosmicWorldLV, nullptr); - - if (!_export_dir.empty()) { - if (_detector == "Flat") { - Export(world, _export_dir, "world.flat.gdml"); - } else if (_detector == "Cosmic") { - Export(cosmic_world, _export_dir, "world.cosmic.gdml"); - } else if (_detector == "Box") { - Export(world, _export_dir, "world.box.gdml"); - } else if (_detector == "MuonMapper") { - Export(world, _export_dir, "world.muon_mapper.gdml"); - } else { - Export(world, _export_dir, "world.prototype.gdml"); - } - } - - std::cout << "Materials: " - << *G4Material::GetMaterialTable() << '\n'; - - // const std::string folder = "detector_geo"; - // const std::string file = "world.gdml"; - // const std::string arg4 = "http://service-spi.web.cern.ch/service-spi/app/releases/GDML/Schema/gdml.xsd"; - - if (_detector == "Cosmic") { - // Construction::Export(cosmic_world, folder, file, arg4); - return cosmic_world; - } else { - // Construction::Export(world, folder, file, arg4); - return world; - } - -} -//---------------------------------------------------------------------------------------------- - -//__Build Detector______________________________________________________________________________ -void Builder::ConstructSDandField() { - if (_detector == "Flat") { - _data_per_event = Flat::Detector::DataPerEvent; - _data_name = Flat::Detector::DataName; - _data_keys = &Flat::Detector::DataKeys; - _data_key_types = &Flat::Detector::DataKeyTypes; - G4SDManager::GetSDMpointer()->AddNewDetector(new Flat::Detector); - } else if (_detector == "Cosmic") { - _data_per_event = Cosmic::Detector::DataPerEvent; - _data_name = Cosmic::Detector::DataName; - _data_keys = &Cosmic::Detector::DataKeys; - _data_key_types = &Cosmic::Detector::DataKeyTypes; - G4SDManager::GetSDMpointer()->AddNewDetector(new Cosmic::Detector); - - G4LogicalVolume* earth_volume = G4LogicalVolumeStore::GetInstance()->GetVolume("ModifiedSandstoneCosmic"); - MultiParticleChangeCrossSection* biasingOperator = new MultiParticleChangeCrossSection; - biasingOperator->AddParticle( "mu+" ); - biasingOperator->AddParticle( "mu-" ); - biasingOperator->AttachTo(earth_volume); - G4cout << " Attaching biasing operator " << biasingOperator->GetName() - << " to logical volume " << earth_volume->GetName() - << G4endl; - } else if (_detector == "Box") { - _data_per_event = Box::Detector::DataPerEvent; - _data_name = Box::Detector::DataName; - _data_keys = &Box::Detector::DataKeys; - _data_key_types = &Box::Detector::DataKeyTypes; - G4SDManager::GetSDMpointer()->AddNewDetector(new Box::Detector); - - G4LogicalVolume* earth_volume = G4LogicalVolumeStore::GetInstance()->GetVolume("ModifiedSandstone"); - MultiParticleChangeCrossSection* biasingOperator = new MultiParticleChangeCrossSection; - biasingOperator->AddParticle( "mu+" ); - biasingOperator->AddParticle( "mu-" ); - biasingOperator->AttachTo(earth_volume); - G4cout << " Attaching biasing operator " << biasingOperator->GetName() - << " to logical volume " << earth_volume->GetName() - << G4endl; - } else if (_detector == "MuonMapper") { - _data_per_event = MuonMapper::Detector::DataPerEvent; - _data_name = MuonMapper::Detector::DataName; - _data_keys = &MuonMapper::Detector::DataKeys; - _data_key_types = &MuonMapper::Detector::DataKeyTypes; - G4SDManager::GetSDMpointer()->AddNewDetector(new MuonMapper::Detector); - } else { - _data_per_event = Prototype::Detector::DataPerEvent; - _data_name = Prototype::Detector::DataName; - _data_keys = &Prototype::Detector::DataKeys; - _data_key_types = &Prototype::Detector::DataKeyTypes; - G4SDManager::GetSDMpointer()->AddNewDetector(new Prototype::Detector); - } -} -//---------------------------------------------------------------------------------------------- - -//__Builder Messenger Set New Value_____________________________________________________________ -void Builder::SetNewValue(G4UIcommand* command, G4String value) { - if (command == _select && value != _detector.c_str()) { - SetDetector(value); - } else if (command == _list) { - std::cout << "Detectors: " << _detectors << "\n"; - } else if (command == _current) { - std::cout << "Current Detector: " << _detector << "\n"; - } -} -//---------------------------------------------------------------------------------------------- - -//__Set Current Detector________________________________________________________________________ -void Builder::SetDetector(const std::string& detector) { - _detector = detector; - Command::Execute("/run/reinitializeGeometry", - "/run/geometryModified", - "/run/initialize", - "/vis/viewer/clearTransients"); -} -//---------------------------------------------------------------------------------------------- - -//__Set Current Detector Save Option____________________________________________________________ -void Builder::SetSaveOption(bool save_option, bool cut_save_option) { - if (_detector == "Flat") { - Flat::Detector::SaveAll = save_option; - } else if (_detector == "Box") { - Box::Detector::SaveAll = save_option; - } else if (_detector == "Cosmic") { - Cosmic::Detector::SaveAll = save_option; - Cosmic::Detector::SaveCut = cut_save_option; - } else if (_detector == "MuonMapper") { - MuonMapper::Detector::SaveAll = save_option; - } else { - Prototype::Detector::SaveAll = save_option; - } -} -//---------------------------------------------------------------------------------------------- - -//__Save some information in a txt file____________________________________________________________ -void Builder::SaveInfo(std::string & prefix){ - if (_detector == "Flat") {;} - else if (_detector == "Box") { - Box::Detector::SaveInfo(prefix); - } -} -//---------------------------------------------------------------------------------------------- - - -//__Get Current Detector Name___________________________________________________________________ -const std::string& Builder::GetDetectorName() { - return _detector; - -} -//---------------------------------------------------------------------------------------------- - -//__Get Current Detector Data Preference________________________________________________________ -bool Builder::IsDetectorDataPerEvent() { - return _data_per_event; -} -//---------------------------------------------------------------------------------------------- - -//__Get Current Detector Data Prefix____________________________________________________________ -const std::string& Builder::GetDetectorDataName() { - return _data_name; -} -//---------------------------------------------------------------------------------------------- - -//__Get Current Detector Data Keys______________________________________________________________ -const Analysis::ROOT::DataKeyList& Builder::GetDetectorDataKeys() { - return *_data_keys; -} -//---------------------------------------------------------------------------------------------- - -//__Get Current Detector Data Key Types_________________________________________________________ -const Analysis::ROOT::DataKeyTypeList& Builder::GetDetectorDataKeyTypes() { - return *_data_key_types; -} -//---------------------------------------------------------------------------------------------- - -//////////////////////////////////////////////////////////////////////////////////////////////// - -//__Sensitive Material Attribute Definition_____________________________________________________ -const G4VisAttributes SensitiveAttributes() { - auto attr = G4VisAttributes(G4Colour(0., 1., 0., 1.0)); - attr.SetForceSolid(true); - return attr; -} - -const G4VisAttributes SensitiveAttributes2() { - auto attr = G4VisAttributes(G4Colour(1., 0., 0., 1.0)); - attr.SetForceSolid(true); - return attr; -} - -const G4VisAttributes HighlightRed() { - auto attr = G4VisAttributes(G4Colour(1., 0., 0.)); - attr.SetForceSolid(true); - return attr; -} - - -//---------------------------------------------------------------------------------------------- - -//__Iron Plating Attributes -// -const G4VisAttributes IronAttributes() { - auto attr = G4VisAttributes(G4Colour(0.5, 0.5, 0.5)); - attr.SetForceSolid(true); - return attr; -} - -const G4VisAttributes AlAttributes() { - auto attr = G4VisAttributes(G4Colour(0.4, 0.6, 0.7)); - attr.SetForceSolid(true); - return attr; -} - - -//__Casing Material Attribute Definition________________________________________________________ -const G4VisAttributes CasingAttributes() { - auto attr = G4VisAttributes(G4Colour(0., 0., 1., 0.2)); - attr.SetForceWireframe(true); - return attr; -} -//---------------------------------------------------------------------------------------------- - -//__Border Attribute Definition_________________________________________________________________ -const G4VisAttributes BorderAttributes() { - return G4VisAttributes(false); -} -//---------------------------------------------------------------------------------------------- - -//__Box Builder_________________________________________________________________________________ -G4Box* Box(const std::string& name, - const double width, - const double height, - const double depth) { - return new G4Box(name, 0.5 * width, 0.5 * height, 0.5 * depth); -} -//---------------------------------------------------------------------------------------------- - -//__Trapezoid Builder___________________________________________________________________________ -G4Trap* Trap(const std::string& name, - const double height, - const double minwidth, - const double maxwidth, - const double depth) { - return new G4Trap(name, - 0.5 * height, 0, 0, - 0.5 * depth, 0.5 * minwidth, 0.5 * minwidth, 0, - 0.5 * depth, 0.5 * maxwidth, 0.5 * maxwidth, 0); -} -//---------------------------------------------------------------------------------------------- - -//__Cylinder Builder____________________________________________________________________________ -G4Tubs* Cylinder(const std::string& name, - const double height, - const double inner_radius, - const double outer_radius) { - return new G4Tubs(name, inner_radius, outer_radius, 0.5 * height, 0, 360*deg); -} -//---------------------------------------------------------------------------------------------- - -//__Volume Builder______________________________________________________________________________ -G4LogicalVolume* Volume(const std::string& name, - G4VSolid* solid, - G4Material* material, - const G4VisAttributes& attr) { - auto out = new G4LogicalVolume(solid, material, name); - out->SetVisAttributes(attr); - return out; -} -//---------------------------------------------------------------------------------------------- - -//__Volume Builder______________________________________________________________________________ -G4LogicalVolume* Volume(G4VSolid* solid, - G4Material* material, - const G4VisAttributes& attr) { - return Volume(solid->GetName(), solid, material, attr); -} -//---------------------------------------------------------------------------------------------- - -//__Volume Builder______________________________________________________________________________ -G4LogicalVolume* Volume(const std::string& name, - G4VSolid* solid, - const G4VisAttributes& attr) { - return Volume(name, solid, Material::Air, attr); -} -//---------------------------------------------------------------------------------------------- - -//__Volume Builder______________________________________________________________________________ -G4LogicalVolume* Volume(G4VSolid* solid, - const G4VisAttributes& attr) { - return Volume(solid, Material::Air, attr); -} -//---------------------------------------------------------------------------------------------- - -//__Box Volume Builder__________________________________________________________________________ -G4LogicalVolume* BoxVolume(const std::string& name, - const double width, - const double height, - const double depth, - G4Material* material, - const G4VisAttributes& attr) { - return Volume( - Box(name, width, height, depth), - material, attr); -} -//---------------------------------------------------------------------------------------------- - -//__Open Box Volume Builder_____________________________________________________________________ -G4LogicalVolume* OpenBoxVolume(const std::string& name, - const double width, - const double height, - const double depth, - const double thickness, - G4Material* material, - const G4VisAttributes& attr) { - auto outer = Box(name, - width, - height, - depth); - auto inner = Box(name, - width - 2 * thickness, - height - 2 * thickness, - depth - 2 * thickness); - return Volume(new G4SubtractionSolid(name, outer, inner), material, attr); -} -//---------------------------------------------------------------------------------------------- - -//__Physical Volume Placer______________________________________________________________________ -G4VPhysicalVolume* PlaceVolume(const std::string& name, - G4LogicalVolume* current, - G4LogicalVolume* parent, - const G4Transform3D& transform) { - return new G4PVPlacement(transform, current, name, parent, false, 0); -} -//---------------------------------------------------------------------------------------------- - -//__Physical Volume Placer______________________________________________________________________ -G4VPhysicalVolume* PlaceVolume(G4LogicalVolume* current, - G4LogicalVolume* parent, - const G4Transform3D& transform) { - return PlaceVolume(current->GetName(), current, parent, transform); -} -//---------------------------------------------------------------------------------------------- - -//__Physical Volume Placer______________________________________________________________________ -G4VPhysicalVolume* PlaceVolume(const std::string& name, - G4VSolid* solid, - G4Material* material, - G4LogicalVolume* parent, - const G4Transform3D& transform) { - return PlaceVolume(name, Volume(name, solid, material), parent, transform); -} -//---------------------------------------------------------------------------------------------- - -//__Physical Volume Placer______________________________________________________________________ -G4VPhysicalVolume* PlaceVolume(G4VSolid* solid, - G4Material* material, - G4LogicalVolume* parent, - const G4Transform3D& transform) { - return PlaceVolume(solid->GetName(), solid, material, parent, transform); -} -//---------------------------------------------------------------------------------------------- - -//__Physical Volume Placer______________________________________________________________________ -G4VPhysicalVolume* PlaceVolume(const std::string& name, - G4LogicalVolume* current, - const G4VisAttributes& attr, - G4LogicalVolume* parent, - const G4Transform3D& transform) { - current->SetVisAttributes(attr); - return PlaceVolume(name, current, parent, transform); -} -//---------------------------------------------------------------------------------------------- - -//__Physical Volume Placer______________________________________________________________________ -G4VPhysicalVolume* PlaceVolume(G4LogicalVolume* current, - const G4VisAttributes& attr, - G4LogicalVolume* parent, - const G4Transform3D& transform) { - current->SetVisAttributes(attr); - return PlaceVolume(current, parent, transform); -} -//---------------------------------------------------------------------------------------------- - -//__Physical Volume Placer______________________________________________________________________ -G4VPhysicalVolume* PlaceVolume(const std::string& name, - G4VSolid* solid, - G4Material* material, - const G4VisAttributes& attr, - G4LogicalVolume* parent, - const G4Transform3D& transform) { - return PlaceVolume(name, Volume(solid, material, attr), parent, transform); -} -//---------------------------------------------------------------------------------------------- - -//__Physical Volume Placer______________________________________________________________________ -G4VPhysicalVolume* PlaceVolume(G4VSolid* solid, - G4Material* material, - const G4VisAttributes& attr, - G4LogicalVolume* parent, - const G4Transform3D& transform) { - return PlaceVolume(Volume(solid, material, attr), parent, transform); -} -//---------------------------------------------------------------------------------------------- - -//__Translation Transformation Generator_______________________________________________________ -G4Transform3D Transform(const double x, - const double y, - const double z) { - return G4Transform3D(G4RotationMatrix(), G4ThreeVector(x, y, z)); -} -//---------------------------------------------------------------------------------------------- - -//__Rotation/Translation Transformation Generator_______________________________________________ -G4Transform3D Transform(const G4ThreeVector& translate, - const G4ThreeVector& axis, - const double angle) { - return G4Transform3D(G4RotationMatrix(axis, angle), translate); -} -//---------------------------------------------------------------------------------------------- - -//__Rotation/Translation Transformation Generator_______________________________________________ -G4Transform3D Transform(const double x, - const double y, - const double z, - const double axisx, - const double axisy, - const double axisz, - const double angle) { - return Transform( - G4ThreeVector(x, y, z), G4ThreeVector(axisx, axisy, axisz), angle); -} -//---------------------------------------------------------------------------------------------- - -//__Rotation Transformation Generator___________________________________________________________ -G4Transform3D Rotate(const double axisx, - const double axisy, - const double axisz, - const double angle) { - return Transform(0, 0, 0, axisx, axisy, axisz, angle); -} -//---------------------------------------------------------------------------------------------- - -//__Matrix Transformation Generator_____________________________________________________________ -G4RotationMatrix Matrix(const double th1, - const double phi1, - const double th2, - const double phi2, - const double th3, - const double phi3) { - const double sinth1 = std::sin(th1); - const double sinth2 = std::sin(th2); - const double sinth3 = std::sin(th3); - auto matrix = G4RotationMatrix(); - matrix.rotateAxes( - G4ThreeVector(sinth1*std::cos(phi1), sinth1*std::sin(phi1), std::cos(th1)), - G4ThreeVector(sinth2*std::cos(phi2), sinth2*std::sin(phi2), std::cos(th2)), - G4ThreeVector(sinth3*std::cos(phi3), sinth3*std::sin(phi3), std::cos(th3))); - - if (matrix != G4RotationMatrix()) matrix.invert(); - return matrix; -} -//---------------------------------------------------------------------------------------------- - -//__Matrix Transformation Generator_____________________________________________________________ -G4RotationMatrix Matrix(const double mxx, - const double mxy, - const double mxz, - const double myx, - const double myy, - const double myz, - const double mzx, - const double mzy, - const double mzz) { - return G4RotationMatrix( - G4ThreeVector(mxx, myx, mzx), - G4ThreeVector(mxy, myy, mzy), - G4ThreeVector(mxz, myz, mzz)); -} -//---------------------------------------------------------------------------------------------- - -//__GDML File Export____________________________________________________________________________ -void Export(const G4LogicalVolume* volume, - const std::string& dir, - const std::string& file, - const std::string& schema) { - util::io::create_directory(dir); - auto path = dir + "/" + file; - if (util::io::path_exists(path)) - util::io::remove_file(path); - - static G4ThreadLocal G4GDMLParser _parser; - _parser.Write(path, volume, true, - !schema.empty() ? schema : G4GDML_DEFAULT_SCHEMALOCATION); -} -//---------------------------------------------------------------------------------------------- - -//__GDML File Export____________________________________________________________________________ -void Export(const G4VPhysicalVolume* volume, - const std::string& dir, - const std::string& file, - const std::string& schema) { - util::io::create_directory(dir); - auto path = dir + "/" + file; - if (util::io::path_exists(path)) - util::io::remove_file(path); - - static G4ThreadLocal G4GDMLParser _parser; - _parser.Write(path, volume, true, - !schema.empty() ? schema : G4GDML_DEFAULT_SCHEMALOCATION); -} -//---------------------------------------------------------------------------------------------- - -} /* namespace Construction */ ///////////////////////////////////////////////////////////////// - -} } /* namespace MATHUSLA::MU */ From 540b0c4d8b6801fd4e6b293d9ce697b1f7e6d368 Mon Sep 17 00:00:00 2001 From: Abdulrahman Mohamed <88555369+Abdulrahma-M@users.noreply.github.com> Date: Mon, 8 May 2023 20:37:22 +0300 Subject: [PATCH 09/12] Delete Box.hh --- Box.hh | 97 ---------------------------------------------------------- 1 file changed, 97 deletions(-) delete mode 100644 Box.hh diff --git a/Box.hh b/Box.hh deleted file mode 100644 index 8e2459a..0000000 --- a/Box.hh +++ /dev/null @@ -1,97 +0,0 @@ -/* include/geometry/Box.hh - * - * Copyright 2018 Brandon Gomes - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -#ifndef MU__GEOMETRY_BOX_HH -#define MU__GEOMETRY_BOX_HH -#pragma once -#include "TROOT.h" -#include "TH1.h" -#include "TTree.h" -#include "TFile.h" -#include "analysis.hh" -#include "geometry/Construction.hh" - -namespace MATHUSLA { namespace MU { - -namespace Box { //////////////////////////////////////////////////////////////////////////////// - -constexpr auto Box_IP_Depth = 85.47*m; - -class Scintillator { -public: - Scintillator(const std::string& name, - const double length, - const double height, - const double width, - const double thickness); - - struct Material { - static G4Material* Casing; - static G4Material* Scintillator; - static void Define(); - private: - Material(); - }; - - double GetLength() const { return _length; } - double GetHeight() const { return _height; } - double GetWidth() const { return _width; } - double GetCasingThickness() const { return _thickness; } - G4LogicalVolume* GetVolume() const { return _lvolume; } - G4VPhysicalVolume* GetSensitiveVolume() const { return _sensitive; } - - void Register(G4VSensitiveDetector* detector); - - G4VPhysicalVolume* PlaceIn(G4LogicalVolume* parent, - const G4Transform3D& transform); - -private: - G4LogicalVolume* _lvolume; - G4VPhysicalVolume* _sensitive; - std::string _name; - double _length, _height, _width, _thickness; -}; - -class Detector : public G4VSensitiveDetector { -public: - Detector(); - - void Initialize(G4HCofThisEvent* event); - G4bool ProcessHits(G4Step* step, G4TouchableHistory*); - void EndOfEvent(G4HCofThisEvent*); - - static const bool DataPerEvent = true; - static const std::string& DataName; - static const Analysis::ROOT::DataKeyList DataKeys; - static const Analysis::ROOT::DataKeyTypeList DataKeyTypes; - static TTree* pre_data; - - static G4VPhysicalVolume* Construct(G4LogicalVolume* world); - static G4VPhysicalVolume* ConstructEarth(G4LogicalVolume* world); - static G4VPhysicalVolume* ConstructModule(G4LogicalVolume* detector, int tag_number, double detector_x, double detector_y, double detector_z); - static G4VPhysicalVolume* ConstructScintillatorLayer(G4LogicalVolume* Module_volume, int module_number, int layer_number, double module_x_displacement, double module_y_displacement, double layer_z_displacement); - static bool SaveAll; - - static void WritePreData(); - static void SaveInfo(std::string & prefix); -}; - -} /* namespace Box */ ////////////////////////////////////////////////////////////////////////// - -} } /* namespace MATHUSLA::MU */ - -#endif /* MU__GEOMETRY_BOX_HH */ From 05b6c9e449f1462e1b17fdb8980bf42e85298625 Mon Sep 17 00:00:00 2001 From: Abdulrahman Mohamed <88555369+Abdulrahma-M@users.noreply.github.com> Date: Mon, 8 May 2023 20:37:34 +0300 Subject: [PATCH 10/12] Delete Box.cc --- Box.cc | 906 --------------------------------------------------------- 1 file changed, 906 deletions(-) delete mode 100644 Box.cc diff --git a/Box.cc b/Box.cc deleted file mode 100644 index a1f2cf0..0000000 --- a/Box.cc +++ /dev/null @@ -1,906 +0,0 @@ -#include "geometry/Box.hh" - -#include -#include - -#include "action.hh" -#include "analysis.hh" -#include "geometry/Earth.hh" -#include "physics/Units.hh" -#include "tracking.hh" -#include "geometry/Cavern.hh" -#include -#include -#include -#include "TROOT.h" -#include "TTree.h" -#include "MuonDataController.hh" - -#include -#include -#include - -#define UNUSED(x) (void)(x) -using dimension = double; - -namespace MATHUSLA { namespace MU { - -namespace Box { ////////////////////////////////////////////////////////////////////////////////////////////////////// - -// double X_POS_STEP; -// double Y_POS_STEP; -// double X_POS_HIT; -// double Y_POS_HIT; - -// TTree* Detector::pre_data = new TTree("pre_data", "pre_data"); - -// void Detector::WritePreData(){ - -// TFile* f = new TFile("box_pre_data.root", "RECREATE"); -// f->cd(); - -// pre_data->Write(); - -// f->Write(); - -//} - -namespace Box::CavernConstruction{ //////////////////////////////////////////////////////////////////////////////////// - -//__Check Between-Ness__________________________________________________________________________ -bool _between(const double min_layer, - const double max_layer, - const double target) { - return min_layer < target && target < max_layer; -} -//---------------------------------------------------------------------------------------------- - -//__Calculate Subtraction of Volumes____________________________________________________________ -G4LogicalVolume* _calculate_modification(const std::string& name, - G4LogicalVolume* earth_component, - const double base_depth, - const double top_depth) { - return Construction::Volume(new G4SubtractionSolid(name, - earth_component->GetSolid(), - Cavern::Volume()->GetSolid(), - Construction::Transform(0, 1.7 * m, -0.5 * (base_depth - top_depth) + Cavern::CenterDepth() - top_depth)), - earth_component->GetMaterial()); -} - -} - - - -namespace { //////////////////////////////////////////////////////////////////////////////////// - -//__Box Sensitive Material______________________________________________________________________ -std::vector _scintillators; -G4LogicalVolume* _steel; -//---------------------------------------------------------------------------------------------- - -//__Box Hit Collection__________________________________________________________________________ -G4ThreadLocal Tracking::HitCollection* _hit_collection; -//---------------------------------------------------------------------------------------------- - -//__Box Specification Variables_________________________________________________________________ - -constexpr int scintillators_per_layer{400}; -constexpr int NMODULES{100}; -constexpr int n_top_layers{6}; -constexpr auto x_edge_length = 99.0*m; -constexpr auto y_edge_length = 99.0*m; -constexpr auto x_displacement = 70.0*m; -constexpr auto y_displacement = -49.5*m; -constexpr auto z_displacement = 6001.5*cm; - -constexpr auto layer_x_edge_length = 9.0*m; -constexpr auto layer_y_edge_length = 9.0*m; - -constexpr auto scint_x_edge_length = 4.5*m; -constexpr auto scint_y_edge_length = 0.045*m; -constexpr auto scintillator_height = 0.02*m; - -constexpr auto steel_height = 0.03*m; - -constexpr auto air_gap = 30*m; - -constexpr auto scintillator_casing_thickness = 0.003*m; - -constexpr auto layer_spacing = 0.8*m; -constexpr auto layer_spacing_top = 0.8*m; // Top layers have different spacing -constexpr auto layer_count = 8UL; - -constexpr auto module_x_edge_length = 9.0*m; -constexpr auto module_y_edge_length = 9.0*m; -constexpr auto module_case_thickness = 0.02*m; - -constexpr auto full_layer_height = scintillator_height + 2*scintillator_casing_thickness; -constexpr auto wall_gap = 0.01*m; -constexpr auto x_edge_increase = 2*full_layer_height + 4*wall_gap; - -constexpr auto layer_w_case = full_layer_height; - -constexpr auto full_module_height = (25.0*m) + 6.0*layer_w_case + 5.0*layer_spacing_top; - -constexpr auto scintillator_z_position = 0.00; - -constexpr auto wall_height = 20*m; - -constexpr int NBEAMLAYERS = 8; -constexpr auto beam_x_edge_length = 0.10*m; -constexpr auto beam_y_edge_length = 0.10*m; -constexpr auto beam_thickness = 0.02*m; - -constexpr auto full_detector_height = full_module_height + steel_height + 2.0*layer_w_case + 1.0*layer_spacing; -constexpr auto half_detector_height = 0.5L * full_detector_height; - -constexpr double layer_z_displacement[8] = {-0.5*full_module_height + (25.0*m -5.0*m - layer_spacing - 2*layer_w_case) + 0.5*layer_w_case, - -0.5*full_module_height + (25.0*m -5.0*m - layer_spacing - 2*layer_w_case) + layer_spacing + 1.5*layer_w_case, - -0.5*full_module_height + (25.0*m ) + 0.5*layer_w_case, - -0.5*full_module_height + (25.0*m ) + layer_spacing_top + 1.5*layer_w_case, - -0.5*full_module_height + (25.0*m ) + 2*layer_spacing_top + 2.5*layer_w_case, - -0.5*full_module_height + (25.0*m ) + 3*layer_spacing_top + 3.5*layer_w_case, - -0.5*full_module_height + (25.0*m ) + 4*layer_spacing_top + 4.5*layer_w_case, - -0.5*full_module_height + (25.0*m ) + 5*layer_spacing_top + 5.5*layer_w_case}; - -// The z distance from the TOP of the floor layer to the BOTTOM of the middle layer -constexpr double floor_z_top = 25.0*m -5.0*m - layer_spacing - 2*layer_w_case; -// The z coordinates of the BOTTOM of all layers in the world -constexpr double layer_z_world[10] = {floor_z_top +layer_spacing+2.*layer_w_case, - floor_z_top +layer_w_case, - -(layer_z_displacement[0]+0.5*full_module_height-floor_z_top-0.5*layer_w_case), - -(layer_z_displacement[1]+0.5*full_module_height-floor_z_top-0.5*layer_w_case), - -(layer_z_displacement[2]+0.5*full_module_height-floor_z_top-0.5*layer_w_case), - -(layer_z_displacement[3]+0.5*full_module_height-floor_z_top-0.5*layer_w_case), - -(layer_z_displacement[4]+0.5*full_module_height-floor_z_top-0.5*layer_w_case), - -(layer_z_displacement[5]+0.5*full_module_height-floor_z_top-0.5*layer_w_case), - -(layer_z_displacement[6]+0.5*full_module_height-floor_z_top-0.5*layer_w_case), - -(layer_z_displacement[7]+0.5*full_module_height-floor_z_top-0.5*layer_w_case) -}; - - -constexpr double module_beam_heights[8] = {floor_z_top, - layer_spacing, - 5.0*m, - layer_spacing_top, - layer_spacing_top, - layer_spacing_top, - layer_spacing_top, - layer_spacing_top}; - - // constexpr double module_beam_z_pos[9] = {-0.50*full_module_height + 0.50*module_beam_heights[0] + layer_w_case, - // -0.50*full_module_height + 2*layer_w_case + layer_spacing + 0.50*module_beam_heights[1], - // -0.50*full_module_height + 3*layer_w_case + 2*layer_spacing + 0.50*module_beam_heights[2], - // -0.50*full_module_height + 20.0*m + layer_w_case + 0.50*module_beam_heights[3], - // -0.50*full_module_height + 20.0*m + 2*layer_w_case + layer_spacing + 0.50*module_beam_heights[4], - // -0.50*full_module_height + 25.0*m + layer_w_case + 0.50*module_beam_heights[5], - // -0.50*full_module_height + 25.0*m + 2*layer_w_case + layer_spacing + 0.50*module_beam_heights[6], - // -0.50*full_module_height + 25.0*m + 3*layer_w_case + 2*layer_spacing + 0.50*module_beam_heights[7], - // -0.50*full_module_height + 25.0*m + 4*layer_w_case + 3*layer_spacing + 0.50*module_beam_heights[8]}; - -constexpr double module_beam_z_pos[8] = {-0.50*full_module_height + 0.50*module_beam_heights[0], - -0.50*full_module_height + floor_z_top + layer_w_case + 0.50*module_beam_heights[1], - -0.50*full_module_height + floor_z_top + 2*layer_w_case + layer_spacing + 0.50*module_beam_heights[2], - -0.50*full_module_height + 25.0*m + layer_w_case + 0.50*module_beam_heights[3], - -0.50*full_module_height + 25.0*m + 2*layer_w_case + 1*layer_spacing_top + 0.50*module_beam_heights[4], - -0.50*full_module_height + 25.0*m + 3*layer_w_case + 2*layer_spacing_top + 0.50*module_beam_heights[5], - -0.50*full_module_height + 25.0*m + 4*layer_w_case + 3*layer_spacing_top + 0.50*module_beam_heights[6], - -0.50*full_module_height + 25.0*m + 5*layer_w_case + 4*layer_spacing_top + 0.50*module_beam_heights[7]}; - - - -// const std::string folder = "detector_geo"; -// const std::string file = "box.gdml"; -// const std::string file2 ="mod.gdml"; -// const std::string file3 ="layer.gdml"; -// const std::string file4 ="earth.gdml"; -// const std::string file5 ="modified.gdml"; -// const std::string arg4 = "http://service-spi.web.cern.ch/service-spi/app/releases/GDML/Schema/gdml.xsd"; - - -auto get_module_x_displacement(int tag_number){ - if (tag_number < 10) return -0.5 * x_edge_length + 0.5*module_x_edge_length; - else if (tag_number < 20) return -0.5*x_edge_length + 0.5*module_x_edge_length + 1.00*(module_x_edge_length + 1.0*m ); - else if (tag_number < 30) return -0.5*x_edge_length + 0.5*module_x_edge_length + 2.00*(module_x_edge_length + 1.0*m ); - else if (tag_number < 40) return -0.5*x_edge_length + 0.5*module_x_edge_length + 3.00*(module_x_edge_length + 1.0*m ); - else if (tag_number < 50) return -0.5*x_edge_length + 0.5*module_x_edge_length + 4.00*(module_x_edge_length + 1.0*m ); - else if (tag_number < 60) return -0.5*x_edge_length + 0.5*module_x_edge_length + 5.00*(module_x_edge_length + 1.0*m ); - else if (tag_number < 70) return -0.5*x_edge_length + 0.5*module_x_edge_length + 6.00*(module_x_edge_length + 1.0*m ); - else if (tag_number < 80) return -0.5*x_edge_length + 0.5*module_x_edge_length + 7.00*(module_x_edge_length + 1.0*m ); - else if (tag_number < 90) return -0.5*x_edge_length + 0.5*module_x_edge_length + 8.00*(module_x_edge_length + 1.0*m ); - else return -0.5*x_edge_length + 0.5*module_x_edge_length + 9.00*(module_x_edge_length + 1.0*m ); -} - -auto get_module_y_displacement(int tag_number){ - return (((double) (tag_number % 10))*(module_y_edge_length + 1.0*m) -0.5 * y_edge_length + 0.5*module_y_edge_length ); -} - -auto get_layer_z_displacement(int layer_number){ - return -1.0*layer_z_displacement[layer_number]; -} - -//---------------------------------------------------------------------------------------------- - -} /* anonymous namespace */ //////////////////////////////////////////////////////////////////// - -//__Box Data Variables__________________________________________________________________________ -const std::string& Detector::DataName = "box_run"; -const Analysis::ROOT::DataKeyList Detector::DataKeys = Analysis::ROOT::DefaultDataKeyList; -const Analysis::ROOT::DataKeyTypeList Detector::DataKeyTypes = Analysis::ROOT::DefaultDataKeyTypeList; -bool Detector::SaveAll = false; -//---------------------------------------------------------------------------------------------- - -//__Detector Constructor________________________________________________________________________ -Detector::Detector() : G4VSensitiveDetector("MATHUSLA/MU/Box") { - collectionName.insert("Box_HC"); - for (auto& scintillator : _scintillators) - scintillator->Register(this); - - // Print out the detector dislacement in z direction - std::cout<<"Layer z displacement [cm]:\n world coord.---- CMS coord.\n"<GetTotalEnergyDeposit(); - - //const auto step_point = step->GetPreStepPoint(); - //const auto position = G4LorentzVector(step_point->GetGlobalTime(), step_point->GetPosition()); - // X_POS_STEP = position.x(); - // Y_POS_STEP = position.y(); - // X_POS_HIT = 0.0; - // Y_POS_HIT = 0.0; - - if (deposit == 0.0L){ - // pre_data->Fill(); - return false; - } - - // X_POS_HIT = position.x(); - // Y_POS_HIT = position.y(); - // pre_data->Fill(); - // const auto track = step->GetTrack(); - // const auto particle = track->GetParticleDefinition(); - // const auto trackID = track->GetTrackID(); - // const auto parentID = track->GetParentID(); - // const auto momentum = G4LorentzVector(step_point->GetTotalEnergy(), step_point->GetMomentum()); - -//////////////////////// - const auto track = step->GetTrack(); - const auto step_point = step->GetPostStepPoint(); - const auto particle = track->GetParticleDefinition(); - const auto trackID = track->GetTrackID(); - const auto parentID = track->GetParentID(); - const auto position = G4LorentzVector(step_point->GetGlobalTime(), step_point->GetPosition()); - const auto momentum = G4LorentzVector(step_point->GetTotalEnergy(), step_point->GetMomentum()); - - //______Tranfomation to CMS Coordinates_____________________________________________________ - const auto transformed_z = -(position.z() - Box_IP_Depth); - const auto position_transformed = G4ThreeVector(position.y(), transformed_z, position.x()); - const auto new_position = G4LorentzVector(step_point->GetGlobalTime(), position_transformed); - - const auto transformed_pz = -(momentum.pz()); - const auto momentum_transformed = G4ThreeVector(momentum.py(), transformed_pz, momentum.px()); - const auto new_momentum = G4LorentzVector(step_point->GetTotalEnergy(), momentum_transformed); - //__________________________________________________________________________________________ - - const auto local_position = new_position.vect() - G4ThreeVector(y_displacement, z_displacement, x_displacement); - - // auto z_index = static_cast(std::floor(+local_position.y() / (layer_spacing + scintillator_height))); - size_t y_index = 0; - -// if (new_position.y() < 6050.0L*cm) { -// y_index = static_cast(std::floor(+local_position.y() / (layer_spacing + scintillator_height))); -// } else if (new_position.y() > 6060.0L*cm && new_position.y() < 6150.0L*cm) { -// y_index = static_cast(std::floor((+local_position.y() - 1.0L*cm) / (layer_spacing + scintillator_height))); -// } else if (new_position.y() > 7900.0L*cm && new_position.y() < 8050.0L*cm) { -// y_index = static_cast(std::floor((+local_position.y() - 1796.0L*cm) / (layer_spacing + scintillator_height))); -// } else if (new_position.y() > 8050.0L*cm && new_position.y() < 8150.0L*cm) { -// y_index = static_cast(std::floor((+local_position.y() - 1797.0L*cm) / (layer_spacing + scintillator_height))); -// } else if (new_position.y() > 8400.0L*cm && new_position.y() < 8550.0L*cm) { -// y_index = static_cast(std::floor((+local_position.y() - 2092.0L*cm) / (layer_spacing + scintillator_height))); -// } else if (new_position.y() > 8560.0L*cm && new_position.y() < 8650.0L*cm) { -// y_index = static_cast(std::floor((+local_position.y() - 2093.0L*cm) / (layer_spacing + scintillator_height))); -// } else if (new_position.y() > 8660.0L*cm && new_position.y() < 8750.0L*cm) { -// y_index = static_cast(std::floor((+local_position.y() - 2094.0L*cm) / (layer_spacing + scintillator_height))); -// } else if (new_position.y() > 8760.0L*cm && new_position.y() < 8850.0L*cm) { -// y_index = static_cast(std::floor((+local_position.y() - 2095.0L*cm) / (layer_spacing + scintillator_height))); -// } else { -// y_index = static_cast(std::floor((+local_position.y() - 2096.0L*cm) / (layer_spacing + scintillator_height))); - // } - - if (new_position.y() <= 6003.5L*cm) { - y_index = static_cast(std::floor(+local_position.y() / (layer_spacing + scintillator_height))); - } else if (new_position.y() >= 6103.5L*cm && new_position.y() <= 6105.5L*cm) { - y_index = static_cast(std::floor((+local_position.y() - 100.0L*cm) / (layer_spacing + scintillator_height))); - } else if (new_position.y() >= 8101.5L*cm && new_position.y() <= 8103.5L*cm) { - y_index = static_cast(std::floor((+local_position.y() - 1996.0L*cm) / (layer_spacing + scintillator_height))); - } else if (new_position.y() >= 8203.5L*cm && new_position.y() <= 8205.5L*cm) { - y_index = static_cast(std::floor((+local_position.y() - 2096.0L*cm) / (layer_spacing + scintillator_height))); - } else if (new_position.y() >= 8601.5L*cm && new_position.y() <= 8603.5L*cm) { - y_index = static_cast(std::floor((+local_position.y() - 2492.0L*cm) / (layer_spacing + scintillator_height))); - } else if (new_position.y() >= 8683.5L*cm && new_position.y() <= 8685.5L*cm) { - y_index = static_cast(std::floor((+local_position.y() - 2492.0L*cm +80.0L*cm) / (layer_spacing + scintillator_height))); - } else if (new_position.y() >= 8765.5L*cm && new_position.y() <= 8767.5L*cm) { - y_index = static_cast(std::floor((+local_position.y() - 2492.0L*cm +2*80.0L*cm) / (layer_spacing + scintillator_height))); - } else if (new_position.y() >= 8847.5L*cm && new_position.y() <= 8849.5L*cm) { - y_index = static_cast(std::floor((+local_position.y() - 2492.0L*cm +3*80.0L*cm) / (layer_spacing + scintillator_height))); - } else if (new_position.y() >= 8929.5L*cm && new_position.y() <= 8931.5L*cm) { - y_index = static_cast(std::floor((+local_position.y() - 2492.0L*cm +4*80.0L*cm) / (layer_spacing + scintillator_height))); - } else if (new_position.y() >= 9011.5L*cm && new_position.y() <= 9013.5L*cm) { - y_index = static_cast(std::floor((+local_position.y() - 2492.0L*cm +5*80.0L*cm) / (layer_spacing + scintillator_height))); - } else { - y_index = static_cast(std::floor((+local_position.y() - 2492.0L*cm +6*80.0L*cm) / (layer_spacing + scintillator_height))); - } - - int _rotation = (1UL + y_index) % 2; - size_t x_index = 0; - size_t z_index = 0; - - if (_rotation == 0){ - - if (new_position.x() <= -4050.0L*cm) { - x_index = static_cast(std::floor(+local_position.x() / (scint_y_edge_length) )); - } else if (new_position.x() >= -3950.0L*cm && new_position.x() <= -3050.0L*cm) { - x_index = static_cast(std::floor((+local_position.x() - 100.0L*cm) / (scint_y_edge_length) )); - } else if (new_position.x() >= -2950.0L*cm && new_position.x() <= -2050.0L*cm) { - x_index = static_cast(std::floor((+local_position.x() - 200.0L*cm) / (scint_y_edge_length) )); - } else if (new_position.x() >= -1950.0L*cm && new_position.x() <= -1050.0L*cm) { - x_index = static_cast(std::floor((+local_position.x() - 300.0L*cm) / (scint_y_edge_length) )); - } else if (new_position.x() >= -950.0L*cm && new_position.x() <= -50.0L*cm) { - x_index = static_cast(std::floor((+local_position.x() - 400.0L*cm) / (scint_y_edge_length) )); - } else if (new_position.x() >= 50.0L*cm && new_position.x() <= 950.0L*cm) { - x_index = static_cast(std::floor((+local_position.x() - 500.0L*cm) / (scint_y_edge_length) )); - } else if (new_position.x() >= 1050.0L*cm && new_position.x() <= 1950.0L*cm) { - x_index = static_cast(std::floor((+local_position.x() - 600.0L*cm) / (scint_y_edge_length) )); - } else if (new_position.x() >= 2050.0L*cm && new_position.x() <= 2950.0L*cm) { - x_index = static_cast(std::floor((+local_position.x() - 700.0L*cm) / (scint_y_edge_length) )); - } else if (new_position.x() >= 3050.0L*cm && new_position.x() <= 3950.0L*cm) { - x_index = static_cast(std::floor((+local_position.x() - 800.0L*cm) / (scint_y_edge_length) )); - } else if (new_position.x() >= 4050.0L*cm && new_position.x() <= 4950.0L*cm) { - x_index = static_cast(std::floor((+local_position.x() - 900.0L*cm) / (scint_y_edge_length) )); - } - if (new_position.z() <= 7900.0L*cm) { - z_index = static_cast(std::floor(+local_position.z() / (scint_x_edge_length) )); - } else if (new_position.z() >= 8000.0L*cm && new_position.z() <= 8900.0L*cm) { - z_index = static_cast(std::floor((+local_position.z() - 100.0L*cm) / (scint_x_edge_length) )); - } else if (new_position.z() >= 9000.0L*cm && new_position.z() <= 9900.0L*cm) { - z_index = static_cast(std::floor((+local_position.z() - 200.0L*cm) / (scint_x_edge_length) )); - } else if (new_position.z() >= 10000.0L*cm && new_position.z() <= 10900.0L*cm) { - z_index = static_cast(std::floor((+local_position.z() - 300.0L*cm) / (scint_x_edge_length) )); - } else if (new_position.z() >= 11000.0L*cm && new_position.z() <= 11900.0L*cm) { - z_index = static_cast(std::floor((+local_position.z() - 400.0L*cm) / (scint_x_edge_length) )); - } else if (new_position.z() >= 12000.0L*cm && new_position.z() <= 12900.0L*cm) { - z_index = static_cast(std::floor((+local_position.z() - 500.0L*cm) / (scint_x_edge_length) )); - } else if (new_position.z() >= 13000.0L*cm && new_position.z() <= 13900.0L*cm) { - z_index = static_cast(std::floor((+local_position.z() - 600.0L*cm) / (scint_x_edge_length) )); - } else if (new_position.z() >= 14000.0L*cm && new_position.z() <= 14900.0L*cm) { - z_index = static_cast(std::floor((+local_position.z() - 700.0L*cm) / (scint_x_edge_length) )); - } else if (new_position.z() >= 15000.0L*cm && new_position.z() <= 15900.0L*cm) { - z_index = static_cast(std::floor((+local_position.z() - 800.0L*cm) / (scint_x_edge_length) )); - } else if (new_position.z() >= 16000.0L*cm && new_position.z() <= 16900.0L*cm) { - z_index = static_cast(std::floor((+local_position.z() - 900.0L*cm) / (scint_x_edge_length) )); - } - - } else if (_rotation == 1){ - - if (new_position.x() <= -4050.0L*cm) { - x_index = static_cast(std::floor(+local_position.x() / (scint_x_edge_length) )); - } else if (new_position.x() >= -3950.0L*cm && new_position.x() <= -3050.0L*cm) { - x_index = static_cast(std::floor((+local_position.x() - 100.0L*cm) / (scint_x_edge_length) )); - } else if (new_position.x() >= -2950.0L*cm && new_position.x() <= -2050.0L*cm) { - x_index = static_cast(std::floor((+local_position.x() - 200.0L*cm) / (scint_x_edge_length) )); - } else if (new_position.x() >= -1950.0L*cm && new_position.x() <= -1050.0L*cm) { - x_index = static_cast(std::floor((+local_position.x() - 300.0L*cm) / (scint_x_edge_length) )); - } else if (new_position.x() >= -950.0L*cm && new_position.x() <= -50.0L*cm) { - x_index = static_cast(std::floor((+local_position.x() - 400.0L*cm) / (scint_x_edge_length) )); - } else if (new_position.x() >= 50.0L*cm && new_position.x() <= 950.0L*cm) { - x_index = static_cast(std::floor((+local_position.x() - 500.0L*cm) / (scint_x_edge_length) )); - } else if (new_position.x() >= 1050.0L*cm && new_position.x() <= 1950.0L*cm) { - x_index = static_cast(std::floor((+local_position.x() - 600.0L*cm) / (scint_x_edge_length) )); - } else if (new_position.x() >= 2050.0L*cm && new_position.x() <= 2950.0L*cm) { - x_index = static_cast(std::floor((+local_position.x() - 700.0L*cm) / (scint_x_edge_length) )); - } else if (new_position.x() >= 3050.0L*cm && new_position.x() <= 3950.0L*cm) { - x_index = static_cast(std::floor((+local_position.x() - 800.0L*cm) / (scint_x_edge_length) )); - } else if (new_position.x() >= 4050.0L*cm && new_position.x() <= 4950.0L*cm) { - x_index = static_cast(std::floor((+local_position.x() - 900.0L*cm) / (scint_x_edge_length) )); - } - if (new_position.z() <= 7900.0L*cm) { - z_index = static_cast(std::floor(+local_position.z() / (scint_y_edge_length) )); - } else if (new_position.z() >= 8000.0L*cm && new_position.z() <= 8900.0L*cm) { - z_index = static_cast(std::floor((+local_position.z() - 100.0L*cm) / (scint_y_edge_length) )); - } else if (new_position.z() >= 9000.0L*cm && new_position.z() <= 9900.0L*cm) { - z_index = static_cast(std::floor((+local_position.z() - 200.0L*cm) / (scint_y_edge_length) )); - } else if (new_position.z() >= 10000.0L*cm && new_position.z() <= 10900.0L*cm) { - z_index = static_cast(std::floor((+local_position.z() - 300.0L*cm) / (scint_y_edge_length) )); - } else if (new_position.z() >= 11000.0L*cm && new_position.z() <= 11900.0L*cm) { - z_index = static_cast(std::floor((+local_position.z() - 400.0L*cm) / (scint_y_edge_length) )); - } else if (new_position.z() >= 12000.0L*cm && new_position.z() <= 12900.0L*cm) { - z_index = static_cast(std::floor((+local_position.z() - 500.0L*cm) / (scint_y_edge_length) )); - } else if (new_position.z() >= 13000.0L*cm && new_position.z() <= 13900.0L*cm) { - z_index = static_cast(std::floor((+local_position.z() - 600.0L*cm) / (scint_y_edge_length) )); - } else if (new_position.z() >= 14000.0L*cm && new_position.z() <= 14900.0L*cm) { - z_index = static_cast(std::floor((+local_position.z() - 700.0L*cm) / (scint_y_edge_length) )); - } else if (new_position.z() >= 15000.0L*cm && new_position.z() <= 15900.0L*cm) { - z_index = static_cast(std::floor((+local_position.z() - 800.0L*cm) / (scint_y_edge_length) )); - } else if (new_position.z() >= 16000.0L*cm && new_position.z() <= 16900.0L*cm) { - z_index = static_cast(std::floor((+local_position.z() - 900.0L*cm) / (scint_y_edge_length) )); - } - } - - const auto x_name = std::to_string(x_index); - const auto z_name = std::to_string(z_index); - - _hit_collection->insert(new Tracking::Hit( - particle, - trackID, - parentID, - std::to_string(1UL + y_index) - + (z_index < 10UL ? "000" + z_name : (z_index < 100UL ? "00" + z_name : (z_index < 1000UL ? "0" + z_name : z_name))) - + (x_index < 10UL ? "000" + x_name : (x_index < 100UL ? "00" + x_name : (x_index < 1000UL ? "0" + x_name : x_name))), - deposit / Units::Energy, - G4LorentzVector(new_position.t() / Units::Time, new_position.vect() / Units::Length), - G4LorentzVector(new_momentum.e() / Units::Energy, new_momentum.vect() / Units::Momentum))); - - return true; -} -//---------------------------------------------------------------------------------------------- - -//__Post-Event Processing_______________________________________________________________________ -void Detector::EndOfEvent(G4HCofThisEvent*) { - if (_hit_collection->GetSize() == 0) - return; - - MuonDataController* controller = MuonDataController::getMuonDataController(); - if(controller->getOn() ==true){ - if(controller->getDecayInEvent() == false){ - return; - } - if(controller->getDecayInZone() == false){ - G4cout<<"Decay In Zone is false"<ExtraDetails()); - root_data.insert(root_data.cend(), gen_particle_data.cbegin(), gen_particle_data.cend()); - root_data.insert(root_data.cend(), extra_gen_data.cbegin(), extra_gen_data.cend()); - - Analysis::ROOT::DataEntry metadata; - metadata.reserve(2UL); - metadata.push_back(collection_data[0UL].size()); - metadata.push_back(gen_particle_data[0UL].size()); - - Analysis::ROOT::FillNTuple(DataName, Detector::DataKeyTypes, metadata, root_data); - if (verboseLevel >= 2 && _hit_collection) - std::cout << *_hit_collection; -} -//---------------------------------------------------------------------------------------------- - -//Build 1 Module for detector -G4VPhysicalVolume* Detector::ConstructScintillatorLayer(G4LogicalVolume* ModuleVolume, int module_number, int layer_number, dimension module_x_displacement, dimension module_y_displacement, dimension layer_z_displacement){ - - auto current = new Scintillator("M" + std::to_string(module_number) + "L" + std::to_string(layer_number), - layer_x_edge_length, - layer_y_edge_length, - full_layer_height, - scintillator_casing_thickness); - - _scintillators.push_back(current); - - UNUSED(module_x_displacement); - UNUSED(module_y_displacement); - - return current->PlaceIn(ModuleVolume, G4Translate3D(0.0, 0.0, layer_z_displacement) ); - //G4Translate3D(0.5*layer_x_edge_length, 0.5*layer_y_edge_length, layer_z_displacement) -} - -G4VPhysicalVolume* Detector::ConstructModule(G4LogicalVolume* DetectorVolume, int tag_number, dimension detector_x, dimension detector_y, dimension detector_z){ - - auto ModuleVolume = Construction::BoxVolume("Module" + std::to_string(tag_number), module_x_edge_length + module_case_thickness, module_y_edge_length + module_case_thickness, full_module_height); - - // auto ModuleCaseVolume = Construction::OpenBoxVolume("Module" + std::to_string(tag_number) + "_Case", module_x_edge_length, module_y_edge_length, full_module_height, - // module_case_thickness, Construction::Material::Iron, *ModuleVisAttr()); - - // Construction::PlaceVolume(ModuleCaseVolume, ModuleVolume, Construction::Transform(0.0, 0.0, 0.0)); - - - for (std::size_t layer{}; layer < layer_count; ++layer) { - auto current = Detector::ConstructScintillatorLayer(ModuleVolume, tag_number, layer, - 0*m, - 0*m, - get_layer_z_displacement(layer)); - UNUSED(current); - } - - UNUSED(detector_x); - UNUSED(detector_y); - UNUSED(detector_z); - - //CONSTRUCTING AND INSERTING STEEL BEAMS - - for (int beam_layer = 0; beam_layer < NBEAMLAYERS; beam_layer++){ - auto BeamL1 = Construction::OpenBoxVolume("Module" + std::to_string(tag_number) + "BL" + std::to_string(beam_layer) + "PL1", beam_x_edge_length, beam_y_edge_length, module_beam_heights[beam_layer], - beam_thickness, Construction::Material::Iron, Construction::CasingAttributes()); - auto BeamL2 = Construction::OpenBoxVolume("Module" + std::to_string(tag_number) + "BL" + std::to_string(beam_layer) + "PL2", beam_x_edge_length, beam_y_edge_length, module_beam_heights[beam_layer], - beam_thickness, Construction::Material::Iron, Construction::CasingAttributes()); - auto BeamR1 = Construction::OpenBoxVolume("Module" + std::to_string(tag_number) + "BL" + std::to_string(beam_layer) + "PR1", beam_x_edge_length, beam_y_edge_length, module_beam_heights[beam_layer], - beam_thickness, Construction::Material::Iron, Construction::CasingAttributes()); - auto BeamR2 = Construction::OpenBoxVolume("Module" + std::to_string(tag_number) + "BL" + std::to_string(beam_layer) + "PR2", beam_x_edge_length, beam_y_edge_length, module_beam_heights[beam_layer], - beam_thickness, Construction::Material::Iron, Construction::CasingAttributes()); - - Construction::PlaceVolume(BeamL1, ModuleVolume, Construction::Transform(-0.50*module_x_edge_length + 0.50*beam_x_edge_length, - -0.50*module_y_edge_length + 0.50*beam_y_edge_length, - -1.0*module_beam_z_pos[beam_layer])); - Construction::PlaceVolume(BeamL2, ModuleVolume, Construction::Transform(-0.50*module_x_edge_length + 0.50*beam_x_edge_length, - 0.50*module_y_edge_length - 0.50*beam_y_edge_length, - -1.0*module_beam_z_pos[beam_layer])); - Construction::PlaceVolume(BeamR1, ModuleVolume, Construction::Transform(0.50*module_x_edge_length - 0.50*beam_x_edge_length, - -0.50*module_y_edge_length + 0.50*beam_y_edge_length, - -1.0* module_beam_z_pos[beam_layer])); - Construction::PlaceVolume(BeamR2, ModuleVolume, Construction::Transform(0.50*module_x_edge_length - 0.50*beam_x_edge_length, - 0.50*module_y_edge_length - 0.50*beam_y_edge_length, - -1.0* module_beam_z_pos[beam_layer])); - } - - - // if (tag_number == 0) { - // std::cout << "ABOUT TO WRITE GDML FOR MODULE" << std::endl; - // Construction::Export(ModuleVolume, folder, file2, arg4 ); - // } - - - return Construction::PlaceVolume(ModuleVolume, DetectorVolume, - Construction::Transform(get_module_x_displacement(tag_number), - get_module_y_displacement(tag_number), - half_detector_height - steel_height - 2.0*layer_w_case - 1*layer_spacing - 0.5*full_module_height, - 0.0, 0.0, 1.0, 0.0)); - - -} - -//__Build Detector______________________________________________________________________________ -G4VPhysicalVolume* Detector::Construct(G4LogicalVolume* world) { - Scintillator::Material::Define(); - _scintillators.clear(); - // pre_data->Branch("X_S", &X_POS_STEP, "X_S/D"); - // pre_data->Branch("Y_S", &Y_POS_STEP, "Y_S/D"); - // pre_data->Branch("X_H", &X_POS_HIT, "X_H/D"); - // pre_data->Branch("Y_H", &Y_POS_HIT, "Y_H/D"); - - auto DetectorVolume = Construction::BoxVolume("Box", x_edge_length + x_edge_increase, y_edge_length, full_detector_height, - Construction::Material::Air, G4VisAttributes::GetInvisible()); - - //DetectorVolume->SetVisAttributes(G4VisAttributes::Invisible); - - for (int module_number = 0; module_number < NMODULES; module_number++){ - Detector::ConstructModule(DetectorVolume, module_number, - 0.5L, //These arguments are not actually used - 0.5L, - 0.0L); - } - - auto first_hermetic_floor = new Scintillator("HF1", - x_edge_length, - y_edge_length, - full_layer_height, - scintillator_casing_thickness); - _scintillators.push_back(first_hermetic_floor); - first_hermetic_floor->PlaceIn(DetectorVolume, G4Translate3D(0.0, 0.0, half_detector_height - 0.5*layer_w_case - steel_height)); - - auto second_hermetic_floor = new Scintillator("HF2", - x_edge_length, - y_edge_length, - full_layer_height, - scintillator_casing_thickness); - _scintillators.push_back(second_hermetic_floor); - second_hermetic_floor->PlaceIn(DetectorVolume, G4Translate3D(0.0, 0.0, half_detector_height - 1.5*layer_w_case - layer_spacing - steel_height)); - - // auto third_hermetic_floor = new Scintillator("HF3", - // x_edge_length, - // y_edge_length, - // full_layer_height, - // scintillator_casing_thickness); - // _scintillators.push_back(third_hermetic_floor); - // third_hermetic_floor->PlaceIn(DetectorVolume, G4Translate3D(0.0, 0.0, half_detector_height - 2.5*layer_w_case - 2*layer_spacing - steel_height)); - - auto hermetic_wall = new Scintillator("HW1", - full_layer_height, - y_edge_length, - wall_height, - scintillator_casing_thickness); - _scintillators.push_back(hermetic_wall); - hermetic_wall->PlaceIn(DetectorVolume, G4Translate3D(-0.5L*x_edge_length - 0.5L*full_layer_height - wall_gap, 0.0, half_detector_height - 0.5L*wall_height)); - - _steel = Construction::BoxVolume("SteelPlate", - x_edge_length, y_edge_length, steel_height, - Construction::Material::Iron, - Construction::CasingAttributes()); - Construction::PlaceVolume(_steel, DetectorVolume, Construction::Transform(0.0, 0.0, half_detector_height - 0.5*steel_height)); - - // Construction::Export(DetectorVolume, folder, file, arg4 ); - - return Construction::PlaceVolume(DetectorVolume, world, - Construction::Transform(0.5L*x_edge_length + x_displacement, 0.5L*y_edge_length + y_displacement, -0.50*full_detector_height + floor_z_top + layer_spacing + 2.*layer_w_case + steel_height)); - -} - - -//---------------------------------------------------------------------------------------------- - -namespace CMS{ - - constexpr auto ____DEFINE_ME____ = 0.0*m; - - constexpr auto earth_total_depth = 4530.0L*cm + 1825.0L*cm + 3645.0L*cm; - - constexpr auto uxc55_cavern_length = 53.0*m; - constexpr auto uxc55_inner_radius = 13.250*m; - constexpr auto uxc55_outer_radius = 14.530*m; - constexpr auto IPDepth = earth_total_depth - uxc55_outer_radius; - constexpr auto _base_depth = earth_total_depth; - constexpr auto access_shaft_x = -14.00*m; - constexpr auto access_shaft_y = 0.00*m; - constexpr auto access_shaft_z = -1.0* uxc55_outer_radius; - - constexpr auto CMSSteelThickness = 1.48L*m; - constexpr auto CMSDetectorLength = 20.00L*m; - constexpr auto CMSDetectorRadius = 8.00L*m; - - constexpr auto AS_Depth = 20.50*m; - constexpr auto AS_Width = 20.50*m; - constexpr auto AS_Height = earth_total_depth - 2 * uxc55_outer_radius; - constexpr auto AS_Thickness = 1.5*m; - - - G4LogicalVolume* EarthVolume() { - using namespace Earth; - - auto earth_box = Construction::Box("", LayerWidthX(), LayerWidthY(), TotalDepth()); - - auto modified = new G4SubtractionSolid("", - earth_box, - Construction::Box("AirBox", x_edge_length + x_edge_increase, y_edge_length, air_gap), - Construction::Transform(0.5L*x_edge_length + x_displacement, - 0.5L*y_edge_length + y_displacement, - 0.5L*(air_gap-Earth::TotalDepth()) -9.50*m )); - - return Construction::Volume(modified); - } - - G4LogicalVolume* SandstoneVolume() { - using namespace Earth; - auto sandstone_box = Construction::Box("", LayerWidthX(), LayerWidthY(), SandstoneDepth()); - - return Construction::Volume(sandstone_box, Material::SiO2, Construction::BorderAttributes()); - } - - long double BaseDepth() { - return _base_depth - Earth::TotalShift(); - } - - long double TopDepth() { - return BaseDepth() - uxc55_outer_radius; - } - long double CenterDepth() { - return BaseDepth() - uxc55_outer_radius; - } - G4Translate3D IPTransform() { - return G4Translate3D(0.0, 0.0, IPDepth); - } - - G4Translate3D Access_Shaft_Transform(){ - return G4Translate3D(access_shaft_x, access_shaft_y, static_cast(access_shaft_z)); - } - - G4Translate3D Cavern_Transform(){ - return G4Translate3D(0, 0, -0.5 * Earth::TotalDepth() + IPDepth); - //* Construction::Rotate(0, 1, 0, 90*deg)) - } - - - G4LogicalVolume* CMSRingVolume() { - return Construction::Volume(Construction::Cylinder("DetectorRing", - CMSDetectorLength, CMSDetectorRadius - CMSSteelThickness, CMSDetectorRadius), - Construction::Material::Iron, - Construction::CasingAttributes()); - } - - G4LogicalVolume* CMSVolume(){ - using namespace Construction; - auto cavern = Cylinder("cavern", uxc55_cavern_length, 0.0*m, uxc55_outer_radius); - - auto access_shaft = Construction::Box("shaft", AS_Width, AS_Depth, AS_Height ); - - return Construction::Volume(new G4UnionSolid("fake_cms", - access_shaft, - cavern, - Construction::Rotate(1, 0, 0, 90*deg) - *G4Translate3D(0.0, -1.0*static_cast(uxc55_outer_radius + 0.5*AS_Height), -0.5*uxc55_cavern_length + 0.5*AS_Width) )); - - } - - //__Calculate Subtraction of Volumes____________________________________________________________ - G4LogicalVolume* _calculate_modification(const std::string& name, - G4LogicalVolume* earth_component, - const double base_depth, - const double top_depth) { - return Construction::Volume(new G4SubtractionSolid(name, - earth_component->GetSolid(), - CMSVolume()->GetSolid(), - Construction::Transform(-0.5*uxc55_cavern_length + 0.5*AS_Width, 0.0, 0.0) - *Construction::Rotate(0, 0, 1, 90*deg) - *Construction::Transform(0.0, 0.0, -0.5 * (base_depth - top_depth) + CenterDepth() - top_depth - uxc55_outer_radius - 0.5*AS_Height)), - earth_component->GetMaterial()); - } - -} // NAMESPACE CMS - - -//__Build Earth for Detector____________________________________________________________________ -G4VPhysicalVolume* Detector::ConstructEarth(G4LogicalVolume* world){ - - //EDIT EARTH TO ACCOMODATE VOLUME - using namespace Box::CavernConstruction; - using namespace Construction; - using namespace CMS; - - Earth::Material::Define(); - - auto earth = CMS::EarthVolume(); - - const auto mix_top = Earth::TotalDepth() - Earth::MixDepth(); - const auto marl_top = mix_top - Earth::MarlDepth(); - const auto sandstone_top = marl_top - Earth::SandstoneDepth(); - - Construction::PlaceVolume(CMS::_calculate_modification("modified_mix", Earth::MixVolume(), - mix_top + Earth::MixDepth(), mix_top), - earth, Earth::MixTransform()); - - Construction::PlaceVolume(CMS::_calculate_modification("modified_marl", Earth::MarlVolume(), - marl_top + Earth::MarlDepth(), marl_top), - earth, Earth::MarlTransform()); - - auto sandstone = CMS::_calculate_modification("modified_sandstone", CMS::SandstoneVolume(), - sandstone_top + Earth::SandstoneDepth(), sandstone_top); - - /////////////// UXC55 AND CMS DETECTOR CONSTRUCTION //////////////////////////////////////// - - auto UXC_55_cavern_solid = Cylinder("UXC55_outer", uxc55_cavern_length, 0.0*m, uxc55_outer_radius); - auto UXC55_outer_solid = Cylinder("UXC55_outer", uxc55_cavern_length, uxc55_inner_radius, uxc55_outer_radius); - auto UXC55_outer_logical = Volume(UXC55_outer_solid, Construction::Material::Concrete, Construction::CasingAttributes()); - auto CMS_Detector_logical = CMSRingVolume(); - auto UXC_55_air_v1 = new G4SubtractionSolid("UXC_55_air_v1", UXC_55_cavern_solid, UXC55_outer_solid); - auto UXC_55_air_v2 = new G4SubtractionSolid("UXC_55_air_v2", UXC_55_air_v1, CMS_Detector_logical->GetSolid()); - auto UXC55_air_logical = Volume("UXC55_air", UXC_55_air_v2, Construction::Material::Air, G4VisAttributes::GetInvisible()); - - Construction::PlaceVolume(UXC55_outer_logical, earth, Cavern_Transform()*Construction::Rotate(0, 1, 0, 90*deg) ); - Construction::PlaceVolume(CMS_Detector_logical, earth, Cavern_Transform()*Construction::Rotate(0, 1, 0, 90*deg) ); - Construction::PlaceVolume(UXC55_air_logical, earth, Cavern_Transform()*Construction::Rotate(0, 1, 0, 90*deg) ); - - - /////////////// ACCESS SHAFT CONSTRUCTION //////////////////////////////////////// - - //MAKE WHOLE SHAFT HERE - - auto Access_Shaft_outer_logical = OpenBoxVolume("Access_Shaft_outer", - AS_Width, - AS_Depth, - AS_Height, - AS_Thickness, - Construction::Material::Concrete, - Construction::CasingAttributes()); - - auto Access_Shaft_Air = BoxVolume("Acess_Shaft_Air", - AS_Width - 2* AS_Thickness, - AS_Depth - 2* AS_Thickness, - AS_Height - 2* AS_Thickness, - Construction::Material::Air, - G4VisAttributes::GetInvisible()); - - Construction::PlaceVolume(Access_Shaft_outer_logical, earth, Access_Shaft_Transform() ); - Construction::PlaceVolume(Access_Shaft_Air, earth, Access_Shaft_Transform()); - - auto modified = Construction::Volume(new G4SubtractionSolid("ModifiedSandstone", - sandstone->GetSolid(), - Construction::Box("AirBox", x_edge_length + x_edge_increase, y_edge_length, air_gap), - Construction::Transform(0.5L*x_edge_length + x_displacement, - 0.5L*y_edge_length + y_displacement, - 0.5L*(air_gap-Earth::SandstoneDepth()) - 9.50*m)), - Earth::Material::SiO2); - - Construction::PlaceVolume(modified, earth, Earth::SandstoneTransform()); - - //PLACE WHOLE THING IN WORLD - - //auto mod = CMS::_calculate_modification("modified_marl", Earth::MarlVolume(), - // marl_top + Earth::MarlDepth(), marl_top); - - - ////export geometry to gdml files - // Construction::Export(CMSVolume(), folder, file5, arg4 ); - // Construction::Export(earth, folder, file4, arg4 ); - - - //// Put Range Cuts on earth volume - // G4Region* cut_region = new G4Region("Earth_Cut_Region"); - // cut_region->AddRootLogicalVolume(earth); - // G4ProductionCuts* cuts = new G4ProductionCuts; - // cuts->SetProductionCut(1.0*m,G4ProductionCuts::GetIndex("gamma")); - // cuts->SetProductionCut(1.0*m,G4ProductionCuts::GetIndex("e-")); - // cuts->SetProductionCut(1.0*m,G4ProductionCuts::GetIndex("e+")); - // cuts->SetProductionCut(1.0*m,G4ProductionCuts::GetIndex("proton")); - // cut_region->SetProductionCuts(cuts); - - - - return Construction::PlaceVolume(earth, world, Earth::Transform()); - -} - - -//---------------------------------------------------------------------------------------------- - -} /* namespace Box */ ////////////////////////////////////////////////////////////////////////// - - - -} } /* namespace MATHUSLA::MU */ - From e8cf758a960b2f6db03fb3d97baea22bdd8b5ddb Mon Sep 17 00:00:00 2001 From: Abdulrahman Mohamed <88555369+Abdulrahma-M@users.noreply.github.com> Date: Mon, 8 May 2023 20:37:43 +0300 Subject: [PATCH 11/12] Delete Construction.hh --- Construction.hh | 294 ------------------------------------------------ 1 file changed, 294 deletions(-) delete mode 100644 Construction.hh diff --git a/Construction.hh b/Construction.hh deleted file mode 100644 index 58a9d3b..0000000 --- a/Construction.hh +++ /dev/null @@ -1,294 +0,0 @@ -#ifndef MU__GEOMETRY_CONSTRUCTION_HH -#define MU__GEOMETRY_CONSTRUCTION_HH -#pragma once - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "analysis.hh" -#include "ui.hh" - -namespace MATHUSLA { namespace MU { - -namespace Construction { /////////////////////////////////////////////////////////////////////// - -namespace Material { /////////////////////////////////////////////////////////////////////////// -extern G4Element* H; -extern G4Element* C; -extern G4Element* N; -extern G4Element* O; -extern G4Element* F; -extern G4Element* Al; -extern G4Element* Si; -extern G4Element* S; -extern G4Element* Ar; -extern G4Element* Ca; -extern G4Material* Air; -extern G4Material* Aluminum; -extern G4Material* Bakelite; -extern G4Material* Copper; -extern G4Material* Concrete; -extern G4Material* Iron; -extern G4Material* PolystyreneFoam; -extern G4Material* Polyvinyltoluene; -} /* namespace Material */ ///////////////////////////////////////////////////////////////////// - -//__Size of The World___________________________________________________________________________ -static constexpr const auto WorldLength = 1600*m; -//---------------------------------------------------------------------------------------------- - -//__Geometry Builder Class______________________________________________________________________ -class Builder : public G4VUserDetectorConstruction, public G4UImessenger { -public: - Builder(const std::string& detector, - const std::string& export_dir, - const bool save_option, - const bool cut_save_option); - G4VPhysicalVolume* Construct(); - void ConstructSDandField(); - - void SetNewValue(G4UIcommand* command, G4String value); - - static const std::string MessengerDirectory; - - static void SetDetector(const std::string& detector); - static void SetSaveOption(const bool save_option, const bool cut_save_option); - static void SaveInfo(std::string& prefix); - - static const std::string& GetDetectorName(); - static bool IsDetectorDataPerEvent(); - static const std::string& GetDetectorDataName(); - static const Analysis::ROOT::DataKeyList& GetDetectorDataKeys(); - static const Analysis::ROOT::DataKeyTypeList& GetDetectorDataKeyTypes(); - -private: - Command::NoArg* _list; - Command::NoArg* _current; - Command::StringArg* _select; -}; -//---------------------------------------------------------------------------------------------- - -//__Sensitive Material Attribute Definition_____________________________________________________ -const G4VisAttributes SensitiveAttributes(); -const G4VisAttributes SensitiveAttributes2(); -const G4VisAttributes HighlightRed(); - -//---------------------------------------------------------------------------------------------- - -//__Casing Material Attribute Definition________________________________________________________ -const G4VisAttributes CasingAttributes(); -//---------------------------------------------------------------------------------------------- - -//__Border Attribute Definition_________________________________________________________________ -const G4VisAttributes BorderAttributes(); -//---------------------------------------------------------------------------------------------- -//__Iron and Aluminum Material Attribute Defintions - -const G4VisAttributes AlAttributes(); -//---------------------------------------------------------------------------------------------- -const G4VisAttributes IronAttributes(); - -//---------------------------------------------------------------------------------------------- -//__Box Builder_________________________________________________________________________________ -G4Box* Box(const std::string& name, - const double width, - const double height, - const double depth); -//---------------------------------------------------------------------------------------------- - -//__Trapezoid Builder___________________________________________________________________________ -G4Trap* Trap(const std::string& name, - const double height, - const double minwidth, - const double maxwidth, - const double depth); -//---------------------------------------------------------------------------------------------- - -//__Cylinder Builder____________________________________________________________________________ -G4Tubs* Cylinder(const std::string& name, - const double height, - const double inner_radius, - const double outer_radius); -//---------------------------------------------------------------------------------------------- - -//__Volume Builder______________________________________________________________________________ -G4LogicalVolume* Volume(const std::string& name, - G4VSolid* solid, - G4Material* material=Material::Air, - const G4VisAttributes& attr=G4VisAttributes()); -//---------------------------------------------------------------------------------------------- - -//__Volume Builder______________________________________________________________________________ -G4LogicalVolume* Volume(G4VSolid* solid, - G4Material* material=Material::Air, - const G4VisAttributes& attr=G4VisAttributes()); -//---------------------------------------------------------------------------------------------- - -//__Volume Builder______________________________________________________________________________ -G4LogicalVolume* Volume(const std::string& name, - G4VSolid* solid, - const G4VisAttributes& attr); -//---------------------------------------------------------------------------------------------- - -//__Volume Builder______________________________________________________________________________ -G4LogicalVolume* Volume(G4VSolid* solid, - const G4VisAttributes& attr); -//---------------------------------------------------------------------------------------------- - -//__Box Volume Builder__________________________________________________________________________ -G4LogicalVolume* BoxVolume(const std::string& name, - const double width, - const double height, - const double depth, - G4Material* material=Material::Air, - const G4VisAttributes& attr=G4VisAttributes()); -//---------------------------------------------------------------------------------------------- - -//__Open Box Volume Builder_____________________________________________________________________ -G4LogicalVolume* OpenBoxVolume(const std::string& name, - const double width, - const double height, - const double depth, - const double thickness, - G4Material* material=Material::Air, - const G4VisAttributes& attr=G4VisAttributes()); -//---------------------------------------------------------------------------------------------- - -//__Physical Volume Placer______________________________________________________________________ -G4VPhysicalVolume* PlaceVolume(const std::string& name, - G4LogicalVolume* current, - G4LogicalVolume* parent, - const G4Transform3D& transform=G4Transform3D()); -//---------------------------------------------------------------------------------------------- - -//__Physical Volume Placer______________________________________________________________________ -G4VPhysicalVolume* PlaceVolume(G4LogicalVolume* current, - G4LogicalVolume* parent, - const G4Transform3D& transform=G4Transform3D()); -//---------------------------------------------------------------------------------------------- - -//__Physical Volume Placer______________________________________________________________________ -G4VPhysicalVolume* PlaceVolume(const std::string& name, - G4VSolid* solid, - G4Material* material, - G4LogicalVolume* parent, - const G4Transform3D& transform=G4Transform3D()); -//---------------------------------------------------------------------------------------------- - -//__Physical Volume Placer______________________________________________________________________ -G4VPhysicalVolume* PlaceVolume(G4VSolid* solid, - G4Material* material, - G4LogicalVolume* parent, - const G4Transform3D& transform=G4Transform3D()); -//---------------------------------------------------------------------------------------------- - -//__Physical Volume Placer______________________________________________________________________ -G4VPhysicalVolume* PlaceVolume(const std::string& name, - G4LogicalVolume* current, - const G4VisAttributes& attr, - G4LogicalVolume* parent, - const G4Transform3D& transform=G4Transform3D()); -//---------------------------------------------------------------------------------------------- - -//__Physical Volume Placer______________________________________________________________________ -G4VPhysicalVolume* PlaceVolume(G4LogicalVolume* current, - const G4VisAttributes& attr, - G4LogicalVolume* parent, - const G4Transform3D& transform=G4Transform3D()); -//---------------------------------------------------------------------------------------------- - -//__Physical Volume Placer______________________________________________________________________ -G4VPhysicalVolume* PlaceVolume(const std::string& name, - G4VSolid* solid, - G4Material* material, - const G4VisAttributes& attr, - G4LogicalVolume* parent, - const G4Transform3D& transform=G4Transform3D()); -//---------------------------------------------------------------------------------------------- - -//__Physical Volume Placer______________________________________________________________________ -G4VPhysicalVolume* PlaceVolume(G4VSolid* solid, - G4Material* material, - const G4VisAttributes& attr, - G4LogicalVolume* parent, - const G4Transform3D& transform=G4Transform3D()); -//---------------------------------------------------------------------------------------------- - -//__Translation Transformation Generator_______________________________________________________ -G4Transform3D Transform(const double x, - const double y, - const double z); -//---------------------------------------------------------------------------------------------- - -//__Rotation/Translation Transformation Generator_______________________________________________ -G4Transform3D Transform(const G4ThreeVector& translate, - const G4ThreeVector& axis, - const double angle); -//---------------------------------------------------------------------------------------------- - -//__Rotation/Translation Transformation Generator_______________________________________________ -G4Transform3D Transform(const double x, - const double y, - const double z, - const double axisx, - const double axisy, - const double axisz, - const double angle); -//---------------------------------------------------------------------------------------------- - -//__Rotation Transformation Generator___________________________________________________________ -G4Transform3D Rotate(const double axisx, - const double axisy, - const double axisz, - const double angle); -//---------------------------------------------------------------------------------------------- - -//__Matrix Transformation Generator_____________________________________________________________ -G4RotationMatrix Matrix(const double th1, - const double phi1, - const double th2, - const double phi2, - const double th3, - const double phi3); -//---------------------------------------------------------------------------------------------- - -//__Matrix Transformation Generator_____________________________________________________________ -G4RotationMatrix Matrix(const double mxx, - const double mxy, - const double mxz, - const double myx, - const double myy, - const double myz, - const double mzx, - const double mzy, - const double mzz); -//---------------------------------------------------------------------------------------------- - -//__GDML File Export____________________________________________________________________________ -void Export(const G4LogicalVolume* volume, - const std::string& dir, - const std::string& file, - const std::string& schema=""); -//---------------------------------------------------------------------------------------------- - -//__GDML File Export____________________________________________________________________________ -void Export(const G4VPhysicalVolume* volume, - const std::string& dir, - const std::string& file, - const std::string& schema=""); -//---------------------------------------------------------------------------------------------- - -} /* namespace Construction */ ///////////////////////////////////////////////////////////////// - -} } /* namespace MATHUSLA::MU */ - -#endif /* MU__GEOMETRY_CONSTRUCTION_HH */ From e2085c96b84127bbd5254f437df05bbf343f64b2 Mon Sep 17 00:00:00 2001 From: Abdulrahman Mohamed <88555369+Abdulrahma-M@users.noreply.github.com> Date: Mon, 8 May 2023 20:55:55 +0300 Subject: [PATCH 12/12] Update Box.cc --- src/geometry/box/Box.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/geometry/box/Box.cc b/src/geometry/box/Box.cc index df4fc0d..80818c0 100644 --- a/src/geometry/box/Box.cc +++ b/src/geometry/box/Box.cc @@ -259,7 +259,7 @@ void Detector::SaveInfo(std::string & prefix){ { //store array contents to text fill for (auto& layer_y : layer_z_world) - fw<< '[' << (-layer_y+Box_IP_Depth+scintillator_casing_thickness)/Units::Length << ", "<< (-layer_y+Box_IP_Depth+scintillator_casing_thickness+scintillator_height)/Units::Length<< "]\n"; + fw<< '[' << (-layer_y+Box_IP_Depth)/Units::Length << ", "<< (-layer_y+Box_IP_Depth+2*scintillator_casing_thickness+scintillator_height)/Units::Length<< "]\n"; fw.close(); } else std::cout << "Problem with opening file";