6 #include <phparameter/PHParameters.h>
10 #include <Geant4/G4Cons.hh>
11 #include <Geant4/G4GenericTrap.hh>
12 #include <Geant4/G4LogicalVolume.hh>
13 #include <Geant4/G4Material.hh>
14 #include <Geant4/G4PVPlacement.hh>
15 #include <Geant4/G4RotationMatrix.hh>
16 #include <Geant4/G4String.hh>
17 #include <Geant4/G4SubtractionSolid.hh>
18 #include <Geant4/G4SystemOfUnits.hh>
19 #include <Geant4/G4ThreeVector.hh>
20 #include <Geant4/G4Transform3D.hh>
21 #include <Geant4/G4Trd.hh>
22 #include <Geant4/G4TwoVector.hh>
23 #include <Geant4/G4Types.hh>
46 , _dx_front(41.44 *
mm)
47 , _dy_front(41.44 *
mm)
48 , _dx_back(48.97454545455 *
mm)
49 , _dy_back(48.97454545455 *
mm)
50 , _dz_crystal(90.000 *
mm)
51 , _crystallogicnameprefix(
"eEcalCrystal")
52 , m_IsActive(GetParams()->get_int_param(
"active"))
53 , m_AbsorberActive(GetParams()->get_int_param(
"absorberactive"))
85 cout <<
"PHG4ProjCrystalCalorimeterDetector: Begin Construction" << endl;
88 if (
GetParams()->get_string_param(
"mappingtower").empty())
90 cout <<
"ERROR in PHG4ProjCrystalCalorimeterDetector: No tower mapping file specified. Abort detector construction." << endl;
91 cout <<
"Please run set_string_param(\"mappingtower\", std::string filename ) first." << endl;
115 ecal_envelope_log,
"CrystalCalorimeter", logicWorld, 0,
false,
OverlapCheck());
125 adjust_width = 0.1258525627 *
mm;
126 adjust_length = 2.4824474402 *
mm;
132 CF_width = 0.18 *
mm;
156 G4double carbon_fiber_width, air_gap_carbon_fiber, air_gap_crystals;
160 G4double dx_front_small = (
_dx_front - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0;
161 G4double dy_front_small = (
_dy_front - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0;
162 G4double dx_back_small = (
_dx_back - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0;
163 G4double dy_back_small = (
_dy_back - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0;
167 std::vector<G4TwoVector> vertices;
169 vertices.push_back(
G4TwoVector(0, dy_front_small));
170 vertices.push_back(
G4TwoVector(dx_front_small, dy_front_small));
171 vertices.push_back(
G4TwoVector(dx_front_small, 0));
174 vertices.push_back(
G4TwoVector(dx_back_small, dy_back_small));
193 G4double TwoByTwo_dx1 = ((2.0 * dx_front_small) + (2.0 * air_gap_carbon_fiber) + air_gap_crystals) / 2.0;
194 G4double TwoByTwo_dx2 = ((2.0 * dx_back_small) + (2.0 * air_gap_carbon_fiber) + air_gap_crystals) / 2.0;
195 G4double TwoByTwo_dy1 = TwoByTwo_dx1;
196 G4double TwoByTwo_dy2 = TwoByTwo_dx2;
219 const int NumberOfIndices = 9;
222 int NumberOfLines = 0;
223 ifstream
in(
GetParams()->get_string_param(
"mapping4x4"));
227 <<
"*******************************************************************" << endl;
228 cout <<
"ERROR in 2 by 2 crystal mapping ";
230 cout <<
"*******************************************************************" << endl
235 while (std::getline(in, unused))
241 in.open(
GetParams()->get_string_param(
"mapping4x4"));
242 G4int j_cry = NumberOfLines;
243 G4int k_cry = NumberOfIndices;
245 double TwoByTwo[j_cry][k_cry];
254 in >> TwoByTwo[j][
k];
267 G4double x_cent, y_cent, z_cent, rot_x, rot_y, rot_z;
273 MappingIndex = TwoByTwo[j][8];
274 if (MappingIndex == 1)
276 j_idx = TwoByTwo[j][0];
277 k_idx = TwoByTwo[j][1];
278 x_cent = TwoByTwo[j][2];
279 y_cent = TwoByTwo[j][3];
280 z_cent = TwoByTwo[j][4];
283 rot_z = TwoByTwo[j][7];
293 int copyno = (j_idx << 16) + k_idx;
301 x_cent = y_cent = z_cent = rot_x = rot_y = rot_z = 0.0;
313 MappingIndex = TwoByTwo[j][8];
314 if (MappingIndex == 4)
316 j_idx = TwoByTwo[j][0];
317 k_idx = TwoByTwo[j][1];
318 x_cent = TwoByTwo[j][2];
319 y_cent = TwoByTwo[j][3];
320 z_cent = TwoByTwo[j][4];
321 rot_x = TwoByTwo[j][5];
322 rot_y = TwoByTwo[j][6];
323 rot_z = TwoByTwo[j][7];
334 int copyno = (j_idx << 16) + k_idx;
342 x_cent = y_cent = z_cent = rot_x = rot_y = 0.0;
354 G4double dx1, dy1, dx2, dy2, dz_whole, carbon_fiber_adjust_width, carbon_fiber_adjust_length;
358 dx1 =
_dx_front + carbon_fiber_adjust_width;
360 dx2 =
_dx_back - carbon_fiber_adjust_width;
362 dz_whole =
_dz_crystal - carbon_fiber_adjust_length;
378 while (j_cry > counter)
380 MappingIndex = TwoByTwo[j][8];
381 if (MappingIndex != 4) j++;
385 x_cent = TwoByTwo[j][2];
386 y_cent = TwoByTwo[j][3];
387 z_cent = TwoByTwo[j][4];
388 rot_x = TwoByTwo[j][5];
389 rot_y = TwoByTwo[j][6];
409 x_cent = TwoByTwo[j][2];
410 y_cent = TwoByTwo[j][3];
411 z_cent = TwoByTwo[j][4];
412 rot_x = TwoByTwo[j][5];
413 rot_y = TwoByTwo[j][6];
416 Crystal_Center =
G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
433 x_cent = TwoByTwo[j][2];
434 y_cent = TwoByTwo[j][3];
435 z_cent = TwoByTwo[j][4];
436 rot_x = TwoByTwo[j][5];
437 rot_y = TwoByTwo[j][6];
440 Crystal_Center =
G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
457 x_cent = TwoByTwo[j][2];
458 y_cent = TwoByTwo[j][3];
459 z_cent = TwoByTwo[j][4];
460 rot_x = TwoByTwo[j][5];
461 rot_y = TwoByTwo[j][6];
464 Crystal_Center =
G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
484 "Carbon_Fiber_logic",
502 "Carbon_Fiber_Shell",
533 G4double carbon_fiber_width, air_gap_carbon_fiber, air_gap_crystals;
537 G4double dx_front_small = (
_dx_front - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0;
538 G4double dy_front_small = (
_dy_front - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0;
539 G4double dx_back_small = (
_dx_back - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0;
540 G4double dy_back_small = (
_dy_back - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0;
544 std::vector<G4TwoVector> vertices;
546 vertices.push_back(
G4TwoVector(0, dy_front_small));
547 vertices.push_back(
G4TwoVector(dx_front_small, dy_front_small));
548 vertices.push_back(
G4TwoVector(dx_front_small, 0));
551 vertices.push_back(
G4TwoVector(dx_back_small, dy_back_small));
570 G4double TwoByTwo_dx1 = ((2.0 * dx_front_small) + (2.0 * air_gap_carbon_fiber) + air_gap_crystals) / 2.0;
571 G4double TwoByTwo_dx2 = ((2.0 * dx_back_small) + (2.0 * air_gap_carbon_fiber) + air_gap_crystals) / 2.0;
572 G4double TwoByTwo_dy1 = TwoByTwo_dx1;
573 G4double TwoByTwo_dy2 = TwoByTwo_dx2;
596 const int NumberOfIndices = 9;
599 int NumberOfLines = 0;
600 ifstream
in(
GetParams()->get_string_param(
"mapping4x4"));
604 <<
"*******************************************************************" << endl;
605 cout <<
"ERROR in 2 by 2 crystal mapping ";
607 cout <<
"*******************************************************************" << endl
612 while (std::getline(in, unused))
618 in.open(
GetParams()->get_string_param(
"mapping4x4"));
619 G4int j_cry = NumberOfLines;
620 G4int k_cry = NumberOfIndices;
622 double TwoByTwo[j_cry][k_cry];
631 in >> TwoByTwo[j][
k];
643 G4double x_cent, y_cent, z_cent, rot_x, rot_y, rot_z;
649 MappingIndex = TwoByTwo[j][8];
650 if (MappingIndex == 1)
652 j_idx = TwoByTwo[j][0];
653 k_idx = TwoByTwo[j][1];
654 x_cent = TwoByTwo[j][2];
655 y_cent = TwoByTwo[j][3];
656 z_cent = TwoByTwo[j][4];
657 rot_x = TwoByTwo[j][5];
658 rot_y = TwoByTwo[j][6];
659 rot_z = TwoByTwo[j][7];
669 int copyno = (j_idx << 16) + k_idx;
679 x_cent = y_cent = z_cent = rot_z = 0.0;
691 j_idx = TwoByTwo[j][0];
692 k_idx = TwoByTwo[j][1];
693 x_cent = TwoByTwo[j][2];
694 y_cent = TwoByTwo[j][3];
695 z_cent = TwoByTwo[j][4];
696 rot_x = TwoByTwo[j][5];
697 rot_y = TwoByTwo[j][6];
698 rot_z = TwoByTwo[j][7];
699 MappingIndex = TwoByTwo[j][8];
710 string Two_by_Two_name =
"TwoByTwo_j_0_k_0";
718 else if (ident == 22)
720 if (MappingIndex == 22)
730 int copyno = (j_idx << 16) + k_idx;
738 else if (ident == 32)
740 if (MappingIndex == 32)
750 int copyno = (j_idx << 16) + k_idx;
761 <<
"Invalid Mapping Type: " << ident << endl;
766 x_cent = y_cent = z_cent = rot_x = rot_y = rot_z = 0.0;
780 (dx1 / 2.00 + 0.0209292929),
781 (dx2 / 2.00 - 0.0209292929),
782 (dy1 / 2.00 + 0.0209292929),
783 (dy2 / 2.0 - 0.0209292929),
808 "Carbon_Fiber_logic",
826 "Carbon_Fiber_Shell",
831 else if (ident == 22)
834 G4double carbon_fiber_adjust_length;
839 (dx1 / 2.0 + (x_adjust)),
840 (dx2 / 2.0 - (x_adjust)),
841 (dy1 + carbon_fiber_adjust_width),
842 (dy2 - carbon_fiber_adjust_width),
843 (dz - carbon_fiber_adjust_length));
847 while (j_cry > counter)
849 MappingIndex = TwoByTwo[j][8];
850 if (MappingIndex != 22) j++;
856 x_cent = TwoByTwo[j][2];
857 y_cent = TwoByTwo[j][3];
858 z_cent = TwoByTwo[j][4];
859 rot_x = TwoByTwo[j][5];
860 rot_y = TwoByTwo[j][6];
861 rot_z = TwoByTwo[j][7];
880 x_cent = TwoByTwo[j][2];
881 y_cent = TwoByTwo[j][3];
882 z_cent = TwoByTwo[j][4];
883 rot_x = TwoByTwo[j][5];
887 Crystal_Center =
G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
902 "Carbon_Fiber_logic",
920 "Carbon_Fiber_Shell",
925 else if (ident == 32)
927 x_cent = 22.6036363636 + 0.18000 * 2.0;
948 G4double carbon_fiber_width, air_gap_carbon_fiber, air_gap_crystals;
952 G4double dx_front_small = (
_dx_front - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0;
954 G4double dx_back_small = (
_dx_back - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0;
957 G4double TwoByTwo_dx1 = ((2.0 * dx_front_small) + (2.0 * air_gap_carbon_fiber) + air_gap_crystals) / 2.0;
958 G4double TwoByTwo_dx2 = ((2.0 * dx_back_small) + (2.0 * air_gap_carbon_fiber) + air_gap_crystals) / 2.0;
959 G4double TwoByTwo_dy1 = TwoByTwo_dx1;
960 G4double TwoByTwo_dy2 = TwoByTwo_dx2;
968 TwoByTwo_dz + 2 * mm);
971 FourByFour_hunk_solid,
981 while (j_cry > counter)
983 MappingIndex = TwoByTwo[j][8];
984 if (MappingIndex != 32) j++;
988 x_cent = TwoByTwo[j][2];
989 y_cent = TwoByTwo[j][3];
990 z_cent = TwoByTwo[j][4];
991 rot_x = TwoByTwo[j][5];
992 rot_y = TwoByTwo[j][6];
993 rot_z = TwoByTwo[j][7];
995 Crystal_Center =
G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
1012 x_cent = TwoByTwo[j][2];
1013 y_cent = TwoByTwo[j][3];
1014 z_cent = TwoByTwo[j][4];
1015 rot_x = TwoByTwo[j][5];
1016 rot_y = TwoByTwo[j][6];
1019 Crystal_Center =
G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
1036 x_cent = TwoByTwo[j][2];
1037 y_cent = TwoByTwo[j][3];
1038 z_cent = TwoByTwo[j][4];
1039 rot_x = TwoByTwo[j][5];
1040 rot_y = TwoByTwo[j][6];
1043 Crystal_Center =
G4ThreeVector(x_cent * mm, y_cent * mm, z_cent * mm);
1058 "Carbon_Fiber_logic",
1076 "Carbon_Fiber_Shell",
1084 <<
"This is an error message which should never be shown. You may have disabled an important portion of this code upsream. Sorry :)" << endl;
1094 G4int NumberOfLines;
1095 const G4int NumberOfIndices = 9;
1100 G4double x_cent, y_cent, z_cent, r_theta, r_phi;
1124 "eEcal_crystal_unit",
1157 x_cent = 22.6036363636 + 0.18000 * 2.0;
1184 G4double dx_front_small = (
_dx_front - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0;
1186 G4double dx_back_small = (
_dx_back - (2.0 * carbon_fiber_width) - (2.0 * air_gap_carbon_fiber) - air_gap_crystals) / 2.0;
1189 G4double TwoByTwo_dx1 = ((2.0 * dx_front_small) + (2.0 * air_gap_carbon_fiber) + air_gap_crystals) / 2.0;
1190 G4double TwoByTwo_dx2 = ((2.0 * dx_back_small) + (2.0 * air_gap_carbon_fiber) + air_gap_crystals) / 2.0;
1191 G4double TwoByTwo_dy1 = TwoByTwo_dx1;
1192 G4double TwoByTwo_dy2 = TwoByTwo_dx2;
1200 TwoByTwo_dz + 5 * mm);
1203 FourByFour_hunk_solid,
1227 if (!datafile.is_open())
1229 datafile.open(
GetParams()->get_string_param(
"mappingtower"));
1233 <<
"*******************************************************************" << endl;
1235 cout <<
"*******************************************************************" << endl
1243 ifstream
in(
GetParams()->get_string_param(
"mappingtower"));
1245 while (std::getline(in, unused))
1248 j_cry = NumberOfLines;
1249 k_cry = NumberOfIndices;
1251 double Crystals[j_cry][k_cry];
1261 datafile >> Crystals[j][
k];
1274 MappingIndex = Crystals[j][8];
1275 if (Crystals[j][1] > k_max) k_max = Crystals[j][1];
1283 MappingIndex = Crystals[j][8];
1285 j_idx = Crystals[j][0];
1286 k_idx = Crystals[j][1];
1290 r_theta = Crystals[j][5];
1291 r_phi = Crystals[j][6];
1292 rot_z = Crystals[j][7];
1294 if (MappingIndex == 16)
1303 string FourByFourName =
"FourByFour_j_" +
to_string(j_idx) +
"_k_" +
to_string(k_idx);
1304 int copyno = (j_idx << 16) + k_idx;
1311 else if (MappingIndex == 32)
1315 if (j_idx == (k_max - 1) / 2 && k_idx == (k_max + 1) / 2) rot_z = rot_z / 2.0;
1322 string FourByFourName =
"FourByFour_j_" +
to_string(j_idx) +
"_k_" +
to_string(k_idx);
1323 int copyno = (j_idx << 16) + k_idx;
1331 else if (MappingIndex == 22)
1340 string FourByFourName =
"FourByFour_j_" +
to_string(j_idx) +
"_k_" +
to_string(k_idx);
1341 int copyno = (j_idx << 16) + k_idx;
1348 else if (MappingIndex == 12)
1357 string FourByFourName =
"FourByFour_j_" +
to_string(j_idx) +
"_k_" +
to_string(k_idx);
1358 int copyno = (j_idx << 16) + k_idx;
1374 MappingIndex = Crystals[j][8];
1376 j_idx = k_max - Crystals[j][0];
1377 k_idx = Crystals[j][1];
1381 r_theta = -1.0 * Crystals[j][5];
1382 r_phi = Crystals[j][6];
1383 rot_z = Crystals[j][7] / (-4.0);
1385 if (MappingIndex == 16)
1394 string FourByFourName =
"FourByFour_j_" +
to_string(j_idx) +
"_k_" +
to_string(k_idx);
1395 int copyno = (j_idx << 16) + k_idx;
1403 else if (MappingIndex == 32)
1407 if (j_idx == (k_max + 1) / 2 && k_idx == (k_max + 1) / 2) rot_z = rot_z * 3.0;
1414 string FourByFourName =
"FourByFour_j_" +
to_string(j_idx) +
"_k_" +
to_string(k_idx);
1415 int copyno = (j_idx << 16) + k_idx;
1422 else if (MappingIndex == 22)
1426 rot_z = Crystals[j][7] * (-1.0);
1433 string FourByFourName =
"FourByFour_j_" +
to_string(j_idx) +
"_k_" +
to_string(k_idx);
1434 int copyno = (j_idx << 16) + k_idx;
1442 else if (MappingIndex == 12)
1446 rot_z = rot_z * -4.0000;
1453 string FourByFourName =
"FourByFour_j_" +
to_string(j_idx) +
"_k_" +
to_string(k_idx);
1454 int copyno = (j_idx << 16) + k_idx;
1470 MappingIndex = Crystals[j][8];
1471 j_idx = k_max - Crystals[j][0];
1472 k_idx = k_max - Crystals[j][1];
1476 r_theta = -1.0 * Crystals[j][5];
1477 r_phi = -1.0 * Crystals[j][6];
1478 rot_z = Crystals[j][7] / 2.0;
1480 if (MappingIndex == 16)
1489 string FourByFourName =
"FourByFour_j_" +
to_string(j_idx) +
"_k_" +
to_string(k_idx);
1490 int copyno = (j_idx << 16) + k_idx;
1498 else if (MappingIndex == 32)
1502 if (j_idx == (k_max + 1) / 2 && k_idx == (k_max - 1) / 2) rot_z = rot_z * -2.0;
1509 string FourByFourName =
"FourByFour_j_" +
to_string(j_idx) +
"_k_" +
to_string(k_idx);
1510 int copyno = (j_idx << 16) + k_idx;
1518 else if (MappingIndex == 22)
1522 rot_z = Crystals[j][7] * (-1.0);
1529 string FourByFourName =
"FourByFour_j_" +
to_string(j_idx) +
"_k_" +
to_string(k_idx);
1530 int copyno = (j_idx << 16) + k_idx;
1538 else if (MappingIndex == 12)
1547 string FourByFourName =
"FourByFour_j_" +
to_string(j_idx) +
"_k_" +
to_string(k_idx);
1548 int copyno = (j_idx << 16) + k_idx;
1563 MappingIndex = Crystals[j][8];
1565 j_idx = Crystals[j][0];
1566 k_idx = k_max - Crystals[j][1];
1570 r_theta = Crystals[j][5];
1571 r_phi = -1.0 * Crystals[j][6];
1572 rot_z = Crystals[j][7] / 4.0;
1574 if (MappingIndex == 16)
1583 string FourByFourName =
"FourByFour_j_" +
to_string(j_idx) +
"_k_" +
to_string(k_idx);
1584 int copyno = (j_idx << 16) + k_idx;
1592 else if (MappingIndex == 32)
1596 if (j_idx == (k_max - 1) / 2 && k_idx == (k_max - 1) / 2) rot_z = rot_z * 3.0;
1603 string FourByFourName =
"FourByFour_j_" +
to_string(j_idx) +
"_k_" +
to_string(k_idx);
1604 int copyno = (j_idx << 16) + k_idx;
1612 else if (MappingIndex == 22)
1616 rot_z = Crystals[j][7] * (1.0);
1623 string FourByFourName =
"FourByFour_j_" +
to_string(j_idx) +
"_k_" +
to_string(k_idx);
1624 int copyno = (j_idx << 16) + k_idx;
1631 else if (MappingIndex == 12)
1635 rot_z = rot_z * -4.0000;
1642 string FourByFourName =
"FourByFour_j_" +
to_string(j_idx) +
"_k_" +
to_string(k_idx);
1643 int copyno = (j_idx << 16) + k_idx;