7 #include <phparameter/PHParameters.h>
23 #include <calobase/RawTowerDefs.h>
24 #include <calobase/RawTowerGeom.h>
25 #include <calobase/RawTowerGeomContainer.h>
26 #include <calobase/RawTowerGeomContainer_Cylinderv1.h>
27 #include <calobase/RawTowerGeomv1.h>
31 #include <Geant4/G4AssemblyVolume.hh>
32 #include <Geant4/G4Box.hh>
33 #include <Geant4/G4ExtrudedSolid.hh>
34 #include <Geant4/G4IntersectionSolid.hh>
35 #include <Geant4/G4LogicalVolume.hh>
36 #include <Geant4/G4Material.hh>
37 #include <Geant4/G4PVPlacement.hh>
38 #include <Geant4/G4RotationMatrix.hh>
39 #include <Geant4/G4String.hh>
40 #include <Geant4/G4SubtractionSolid.hh>
41 #include <Geant4/G4SystemOfUnits.hh>
42 #include <Geant4/G4ThreeVector.hh>
43 #include <Geant4/G4Transform3D.hh>
44 #include <Geant4/G4Tubs.hh>
45 #include <Geant4/G4TwoVector.hh>
46 #include <Geant4/G4UserLimits.hh>
47 #include <Geant4/G4VPhysicalVolume.hh>
48 #include <Geant4/G4VSolid.hh>
50 #pragma GCC diagnostic push
51 #pragma GCC diagnostic ignored "-Wshadow"
52 #pragma GCC diagnostic ignored "-Wpedantic"
53 #include <CGAL/Boolean_set_operations_2.h>
54 #include <CGAL/Circular_kernel_intersections.h>
55 #include <CGAL/Exact_circular_kernel_2.h>
56 #include <CGAL/Object.h>
57 #include <CGAL/point_generators_2.h>
58 #pragma GCC diagnostic pop
60 #include <boost/lexical_cast.hpp>
61 #include <boost/tokenizer.hpp>
72 using Circle_2 = CGAL::Circle_2<PHG4OuterHcalDetector::Circular_k>;
74 using Line_2 = CGAL::Line_2<PHG4OuterHcalDetector::Circular_k>;
75 using Segment_2 = CGAL::Segment_2<PHG4OuterHcalDetector::Circular_k>;
88 , m_InnerRadius(m_Params->get_double_param(
"inner_radius") *
cm)
89 , m_OuterRadius(m_Params->get_double_param(
"outer_radius") *
cm)
90 , m_SizeZ(m_Params->get_double_param(
"size_z") *
cm)
91 , m_ScintiTileZ(m_SizeZ)
92 , m_ScintiTileThickness(m_Params->get_double_param(
"scinti_tile_thickness") *
cm)
93 , m_ScintiGap(m_Params->get_double_param(
"scinti_gap") *
cm)
94 , m_ScintiInnerRadius(m_Params->get_double_param(
"scinti_inner_radius") *
cm)
95 , m_ScintiOuterRadius(m_Params->get_double_param(
"scinti_outer_radius") *
cm)
96 , m_TiltAngle(m_Params->get_double_param(
"tilt_angle") *
deg)
97 , m_EnvelopeInnerRadius(m_InnerRadius)
98 , m_EnvelopeOuterRadius(m_OuterRadius)
99 , m_EnvelopeZ(m_SizeZ)
100 , m_NumScintiPlates(m_Params->get_int_param(PHG4HcalDefs::
scipertwr) * m_Params->get_int_param(
"n_towers"))
101 , m_NumScintiTiles(m_Params->get_int_param(
"n_scinti_tiles"))
102 , m_ActiveFlag(m_Params->get_int_param(
"active"))
103 , m_AbsorberActiveFlag(m_Params->get_int_param(
"absorberactive"))
104 , m_ScintiLogicNamePrefix(
"HcalOuterScinti")
147 Line_2 s2(p_in_1, p_upperedge);
151 Circle_2 outer_circle(sc1, sc2, sc3);
152 std::vector<CGAL::Object> res;
153 CGAL::intersection(outer_circle, perp, std::back_inserter(res));
155 std::vector<CGAL::Object>::const_iterator iter;
156 for (iter = res.begin(); iter != res.end(); ++iter)
158 CGAL::Object obj = *iter;
159 if (
const std::pair<CGAL::Circular_arc_point_2<PHG4OuterHcalDetector::Circular_k>,
unsigned> *
point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4OuterHcalDetector::Circular_k>,
unsigned>>(&obj))
161 if (CGAL::to_double(
point->first.x()) > CGAL::to_double(p_upperedge.x()))
163 double deltax = CGAL::to_double(
point->first.x()) - CGAL::to_double(p_upperedge.x());
164 double deltay = CGAL::to_double(
point->first.y()) - CGAL::to_double(p_upperedge.y());
173 std::cout <<
"CGAL::Object type not pair..." << std::endl;
180 Line_2 s3(p_in_1, p_loweredge);
181 Line_2 l_lower = s3.perpendicular(p_loweredge);
183 Circle_2 inner_circle(ic1, ic2, ic3);
185 CGAL::intersection(inner_circle, l_lower, std::back_inserter(res));
190 for (iter = res.begin(); iter != res.end(); ++iter)
192 CGAL::Object obj = *iter;
193 if (
const std::pair<CGAL::Circular_arc_point_2<PHG4OuterHcalDetector::Circular_k>,
unsigned> *
point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4OuterHcalDetector::Circular_k>,
unsigned>>(&obj))
195 if (CGAL::to_double(
point->first.x()) > minx)
197 minx = CGAL::to_double(
point->first.x());
198 double deltax = CGAL::to_double(
point->first.x()) - CGAL::to_double(p_loweredge.x());
199 double deltay = CGAL::to_double(
point->first.y()) - CGAL::to_double(p_loweredge.y());
223 double xcoord =
m_ScintiGap / 2. * cos(angle_mid_scinti /
rad) + mid_radius;
224 double ycoord =
m_ScintiGap / 2. * sin(angle_mid_scinti /
rad) + 0;
226 Line_2 s2(p_in_1, p_loweredge);
229 Circle_2 inner_circle(sc1, sc2, sc3);
230 std::vector<CGAL::Object> res;
231 CGAL::intersection(inner_circle, perp, std::back_inserter(res));
233 std::vector<CGAL::Object>::const_iterator iter;
234 for (iter = res.begin(); iter != res.end(); ++iter)
236 CGAL::Object obj = *iter;
237 if (
const std::pair<CGAL::Circular_arc_point_2<PHG4OuterHcalDetector::Circular_k>,
unsigned> *
point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4OuterHcalDetector::Circular_k>,
unsigned>>(&obj))
239 if (CGAL::to_double(
point->first.x()) > 0)
247 std::cout <<
"CGAL::Object type not pair..." << std::endl;
251 Circle_2 outer_circle(so1, so2, so3);
253 CGAL::intersection(outer_circle, perp, std::back_inserter(res));
255 for (iter = res.begin(); iter != res.end(); ++iter)
257 CGAL::Object obj = *iter;
258 if (
const std::pair<CGAL::Circular_arc_point_2<PHG4OuterHcalDetector::Circular_k>,
unsigned> *
point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4OuterHcalDetector::Circular_k>,
unsigned>>(&obj))
260 if (CGAL::to_double(
point->first.x()) > CGAL::to_double(p_loweredge.x()))
268 std::cout <<
"CGAL::Object type not pair..." << std::endl;
275 double xmidpoint = cos(phi_midpoint) * mid_radius;
276 double ymidpoint = sin(phi_midpoint) * mid_radius;
278 angle_mid_scinti = (M_PI / 2. - phi_midpoint) - (M_PI / 2. +
m_TiltAngle /
rad);
279 double xcoordup = xmidpoint -
m_ScintiGap / 2. * sin(angle_mid_scinti /
rad);
280 double ycoordup = ymidpoint -
m_ScintiGap / 2. * cos(angle_mid_scinti /
rad);
285 Line_2 sup(mid_upperscint, p_upperedge);
286 Line_2 perpA = sup.perpendicular(p_upperedge);
288 Circle_2 inner_circleA(sc1A, sc2A, sc3A);
289 std::vector<CGAL::Object> resA;
290 CGAL::intersection(inner_circleA, perpA, std::back_inserter(resA));
291 std::vector<CGAL::Object>::const_iterator iterA;
293 for (iterA = resA.begin(); iterA != resA.end(); ++iterA)
295 CGAL::Object obj = *iterA;
296 if (
const std::pair<CGAL::Circular_arc_point_2<PHG4OuterHcalDetector::Circular_k>,
unsigned> *
point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4OuterHcalDetector::Circular_k>,
unsigned>>(&obj))
298 if (CGAL::to_double(
point->first.x()) > pxmax)
300 pxmax = CGAL::to_double(
point->first.x());
307 std::cout <<
"CGAL::Object type not pair..." << std::endl;
311 Circle_2 outer_circleA(so1A, so2A, so3A);
313 CGAL::intersection(outer_circleA, perpA, std::back_inserter(resA));
314 for (iterA = resA.begin(); iterA != resA.end(); ++iterA)
316 CGAL::Object obj = *iterA;
317 if (
const std::pair<CGAL::Circular_arc_point_2<PHG4OuterHcalDetector::Circular_k>,
unsigned> *
point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4OuterHcalDetector::Circular_k>,
unsigned>>(&obj))
319 if (CGAL::to_double(
point->first.x()) > CGAL::to_double(p_loweredge.x()))
327 std::cout <<
"CGAL::Object type not pair..." << std::endl;
333 G4TwoVector
v1(CGAL::to_double(upperleft.x()), CGAL::to_double(upperleft.y()));
334 G4TwoVector
v2(CGAL::to_double(upperright.x()), CGAL::to_double(upperright.y()));
335 G4TwoVector
v3(CGAL::to_double(lowerright.x()), CGAL::to_double(lowerright.y()));
336 G4TwoVector
v4(CGAL::to_double(lowerleft.x()), CGAL::to_double(lowerleft.y()));
337 std::vector<G4TwoVector> vertexes;
338 vertexes.push_back(v1);
339 vertexes.push_back(v2);
340 vertexes.push_back(v3);
341 vertexes.push_back(v4);
342 G4TwoVector
zero(0, 0);
343 G4VSolid *steel_plate_uncut =
new G4ExtrudedSolid(
"SteelPlateUnCut",
353 return steel_plate_uncut;
355 G4RotationMatrix *rotm =
new G4RotationMatrix();
356 rotm->rotateX(-90 *
deg);
357 G4VSolid *steel_firstcut_solid =
new G4SubtractionSolid(
"SteelPlateFirstCut", steel_plate_uncut,
m_SteelCutoutForMagnetG4Solid, rotm, G4ThreeVector(0, 0, 0));
363 rotm =
new G4RotationMatrix();
364 rotm->rotateX(90 *
deg);
365 G4VSolid *steel_cut_solid =
new G4SubtractionSolid(
"SteelPlateCut", steel_firstcut_solid,
m_SteelCutoutForMagnetG4Solid, rotm, G4ThreeVector(0, 0, 0));
368 return steel_cut_solid;
373 Line_2 secant(lowleft, upleft);
376 double xmid = (CGAL::to_double(lowleft.x()) + CGAL::to_double(upleft.x())) / 2.;
377 double ymid = (CGAL::to_double(lowleft.y()) + CGAL::to_double(upleft.y())) / 2.;
379 Line_2 sekperp = secant.perpendicular(midpoint);
381 Circle_2 inner_circle(sc1, sc2, sc3);
382 std::vector<CGAL::Object> res;
383 CGAL::intersection(inner_circle, sekperp, std::back_inserter(res));
384 std::vector<CGAL::Object>::const_iterator iter;
387 for (iter = res.begin(); iter != res.end(); ++iter)
389 CGAL::Object obj = *iter;
390 if (
const std::pair<CGAL::Circular_arc_point_2<PHG4OuterHcalDetector::Circular_k>,
unsigned> *
point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4OuterHcalDetector::Circular_k>,
unsigned>>(&obj))
392 if (CGAL::to_double(
point->first.x()) > pxmax)
394 pxmax = CGAL::to_double(
point->first.x());
401 std::cout <<
"CGAL::Object type not pair..." << std::endl;
404 Line_2 leftside = sekperp.perpendicular(tangtouch);
405 CGAL::Object result = CGAL::intersection(upedge, leftside);
410 result = CGAL::intersection(lowedge, leftside);
428 G4LogicalVolume *hcal_envelope_log =
new G4LogicalVolume(hcal_envelope_cylinder, Air, G4String(
"OuterHcal_envelope"),
nullptr,
nullptr,
nullptr);
429 G4RotationMatrix hcal_rotm;
457 boost::char_separator<char> sep(
"_");
458 boost::tokenizer<boost::char_separator<char>> tok((*it)->GetName(), sep);
459 boost::tokenizer<boost::char_separator<char>>::const_iterator tokeniter;
460 for (tokeniter = tok.begin(); tokeniter != tok.end(); ++tokeniter)
462 if (*tokeniter ==
"impr")
465 if (tokeniter != tok.end())
467 int layer_id = boost::lexical_cast<
int>(*tokeniter);
472 int tower_id = (*it)->GetCopyNo() - layer_id;
474 std::pair<int, int> layer_twr = std::make_pair(layer_id, tower_id);
478 std::cout <<
"invalid scintillator row " << layer_id
485 std::cout <<
PHWHERE <<
" Error parsing " << (*it)->GetName()
486 <<
" for mother volume number " << std::endl;
523 std::ostringstream
name;
533 sumshift = sumshift - 2 * shiftup;
545 double xp = cos(phi) * middlerad;
546 double yp = sin(phi) * middlerad;
553 phi = -atan(yo / xo);
559 phi = -atan(yo / xo);
564 G4RotationMatrix *Rot =
new G4RotationMatrix();
565 double ypos = sin(phi) * middlerad;
566 double xpos = cos(phi) * middlerad;
573 G4ThreeVector g4vec(xpos, ypos, 0);
581 Rot =
new G4RotationMatrix();
582 Rot->rotateZ(-phi *
rad);
584 name <<
"OuterHcalSteel_" <<
i;
603 std::ostringstream
name;
610 double inner_offset =
offset;
620 std::fill_n(zsteelcut, 4, NAN);
622 double steel_offset = 1 *
cm + steel_overhang;
625 double steel_inner_offset = steel_offset;
626 xsteelcut[0] = steel_x_inner + steel_magnet_cutout_x;
627 xsteelcut[1] = xsteelcut[0];
629 xsteelcut[3] = xsteelcut[2];
636 inner_offset = offset - magnet_cutout_x;
658 z[0] += scinti_gap_neighbor / 2.;
659 z[1] += scinti_gap_neighbor / 2.;
660 z[2] -= scinti_gap_neighbor / 2.;
661 z[3] -= scinti_gap_neighbor / 2.;
665 z[0] =
x_at_y(leftsidelow, leftsidehigh, x[0]);
667 z[1] =
x_at_y(leftsidelow, leftsidehigh, x[1]);
671 z[2] =
x_at_y(rightsidelow, rightsidehigh, x[2]);
673 z[3] =
x_at_y(rightsidelow, rightsidehigh, x[3]);
677 double x0 =
m_InnerRadius - (steel_inner_offset - steel_magnet_cutout_x);
678 double z0 =
x_at_y(leftsidelow, leftsidehigh, x0);
681 zsteelcut[3] =
x_at_y(leftsidelow, leftsidehigh, xpos);
684 double z2 =
x_at_y(rightsidelow, rightsidehigh, x2);
685 zsteelcut[1] = z2 + 1 *
cm;
686 zsteelcut[2] = z2 + 1 *
cm;
687 std::vector<G4TwoVector> vertexes;
688 for (
int j = 0;
j < 4;
j++)
690 G4TwoVector
v(x[
j], z[j]);
691 vertexes.push_back(v);
693 G4TwoVector
zero(0, 0);
695 G4VSolid *scinti =
new G4ExtrudedSolid(
"ScintillatorTile",
700 G4RotationMatrix *rotm =
new G4RotationMatrix();
701 rotm->rotateX(-90 *
deg);
703 name <<
"scintillator_" <<
i <<
"_left";
707 rotm =
new G4RotationMatrix();
708 rotm->rotateX(90 *
deg);
710 name <<
"scintillator_" << i <<
"_right";
725 std::vector<G4TwoVector> vertexes;
726 for (
int j = 0;
j < 4;
j++)
728 if (!isfinite(zsteelcut[
j]))
732 G4TwoVector
v(xsteelcut[j], zsteelcut[j]);
733 vertexes.push_back(v);
735 G4TwoVector
zero(0, 0);
749 x[0] = CGAL::to_double(p0.x());
750 x[1] = CGAL::to_double(p1.x());
752 double newx = fabs(x[0]) + fabs(x[1]);
756 CGAL::Object result = CGAL::intersection(l, seg);
759 xret = CGAL::to_double(ipoint->x());
763 std::cout <<
PHWHERE <<
" failed for y = " << yin << std::endl;
764 std::cout <<
"p0(x): " << CGAL::to_double(p0.x()) <<
", p0(y): " << CGAL::to_double(p0.y()) << std::endl;
765 std::cout <<
"p1(x): " << CGAL::to_double(p1.x()) <<
", p1(y): " << CGAL::to_double(p1.y()) << std::endl;
779 G4AssemblyVolume *assmeblyvol =
new G4AssemblyVolume();
780 std::ostringstream
name;
787 G4UserLimits *g4userlimits =
nullptr;
788 if (isfinite(steplimits))
790 g4userlimits =
new G4UserLimits(steplimits);
794 assmeblyvol->AddPlacedVolume(scinti_tile_logic, g4vec,
nullptr);
809 <<
" cm" << std::endl;
816 <<
" cm" << std::endl;
823 <<
" cm" << std::endl;
830 <<
" cm" << std::endl;
837 <<
" cm" << std::endl;
844 <<
" cm" << std::endl;
851 <<
" cm" << std::endl;
869 std::cout <<
"both number of crossings and tilt angle are set" << std::endl;
870 std::cout <<
"using number of crossings to determine tilt angle" << std::endl;
878 Circle_2 inner_circle(pin1, pin2, pin3);
880 Circle_2 mid_circle(pmid1, pmid2, pmid3);
882 Circle_2 outer_circle(pout1, pout2, pout3);
883 Line_2 l_up(pnull, phightmp);
884 std::vector<CGAL::Object> res;
885 CGAL::intersection(outer_circle, l_up, std::back_inserter(res));
887 std::vector<CGAL::Object>::const_iterator iter;
888 for (iter = res.begin(); iter != res.end(); ++iter)
890 CGAL::Object obj = *iter;
891 if (
const std::pair<CGAL::Circular_arc_point_2<PHG4OuterHcalDetector::Circular_k>,
unsigned> *
point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4OuterHcalDetector::Circular_k>,
unsigned>>(&obj))
893 if (CGAL::to_double(
point->first.x()) > 0)
901 std::cout <<
"CGAL::Object type not pair..." << std::endl;
905 Line_2 l_right(plow, upperright);
908 CGAL::intersection(mid_circle, l_right, std::back_inserter(res));
909 for (iter = res.begin(); iter != res.end(); ++iter)
911 CGAL::Object obj = *iter;
912 if (
const std::pair<CGAL::Circular_arc_point_2<PHG4OuterHcalDetector::Circular_k>,
unsigned> *
point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4OuterHcalDetector::Circular_k>,
unsigned>>(&obj))
914 if (CGAL::to_double(
point->first.x()) > 0)
922 std::cout <<
"CGAL::Object type not pair..." << std::endl;
927 double ll = sqrt((CGAL::to_double(midpoint.x()) -
m_InnerRadius) * (CGAL::to_double(midpoint.x()) -
m_InnerRadius) + CGAL::to_double(midpoint.y()) * CGAL::to_double(midpoint.y()));
928 double upside = sqrt(CGAL::to_double(midpoint.x()) * CGAL::to_double(midpoint.x()) + CGAL::to_double(midpoint.y()) * CGAL::to_double(midpoint.y()));
932 tiltangle = tiltangle *
rad;
955 Circle_2 inner_circle(sc1, sc2, sc3);
956 std::vector<CGAL::Object> res;
957 CGAL::intersection(inner_circle, s2, std::back_inserter(res));
961 <<
" too large, no intersection with inner radius" << std::endl;
969 std::cout <<
"Outer Hcal Detector:" << std::endl;
970 if (what ==
"ALL" || what ==
"VOLUME")
987 std::cout <<
"could not locate volume " << volume->GetName()
988 <<
" in Inner Hcal scintillator map" << std::endl;
1002 std::cout <<
PHWHERE <<
"Run Node missing, exiting." << std::endl;
1027 double geom_ref_radius = innerrad + thickness / 2.;
1029 if (!std::isfinite(phistart))
1031 std::cout <<
PHWHERE <<
" phistart is not finite: " << phistart
1032 <<
", exiting now (this will crash anyway)" << std::endl;
1038 std::pair<double, double> range = std::make_pair(phiend, phistart);
1046 double etahibound = etalowbound +
1048 std::pair<double, double> range = std::make_pair(etalowbound, etahibound);
1050 etalowbound = etahibound;
1067 std::cout <<
"IHCalDetector::InitRun - Tower geometry " << key <<
" already exists" << std::endl;
1070 if (fabs(tg->get_center_x() -
x) > 1
e-4)
1072 std::cout <<
"IHCalDetector::InitRun - Fatal Error - duplicated Tower geometry " << key <<
" with existing x = " << tg->
get_center_x() <<
" and expected x = " << x
1077 if (fabs(tg->get_center_y() -
y) > 1
e-4)
1079 std::cout <<
"IHCalDetector::InitRun - Fatal Error - duplicated Tower geometry " << key <<
" with existing y = " << tg->get_center_y() <<
" and expected y = " << y
1083 if (fabs(tg->get_center_z() -
z) > 1
e-4)
1085 std::cout <<
"IHCalDetector::InitRun - Fatal Error - duplicated Tower geometry " << key <<
" with existing z= " << tg->get_center_z() <<
" and expected z = " <<
z
1094 std::cout <<
"IHCalDetector::InitRun - building tower geometry " << key <<
"" << std::endl;
1099 tg->set_center_x(x);
1100 tg->set_center_y(y);
1101 tg->set_center_z(
z);