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/G4SystemOfUnits.hh>
41 #include <Geant4/G4ThreeVector.hh>
42 #include <Geant4/G4Transform3D.hh>
43 #include <Geant4/G4Tubs.hh>
44 #include <Geant4/G4TwoVector.hh>
45 #include <Geant4/G4Types.hh>
46 #include <Geant4/G4UserLimits.hh>
47 #include <Geant4/G4VSolid.hh>
49 #pragma GCC diagnostic push
50 #pragma GCC diagnostic ignored "-Wshadow"
51 #pragma GCC diagnostic ignored "-Wpedantic"
52 #include <CGAL/Boolean_set_operations_2.h>
53 #include <CGAL/Circular_kernel_intersections.h>
54 #include <CGAL/Exact_circular_kernel_2.h>
55 #include <CGAL/Object.h>
56 #include <CGAL/point_generators_2.h>
57 #pragma GCC diagnostic pop
59 #include <boost/lexical_cast.hpp>
60 #include <boost/tokenizer.hpp>
71 using Circle_2 = CGAL::Circle_2<PHG4InnerHcalDetector::Circular_k>;
73 using Line_2 = CGAL::Line_2<PHG4InnerHcalDetector::Circular_k>;
74 using Segment_2 = CGAL::Segment_2<PHG4InnerHcalDetector::Circular_k>;
84 , m_Params(parameters)
85 , m_InnerRadius(m_Params->get_double_param(
"inner_radius") *
cm)
86 , m_OuterRadius(m_Params->get_double_param(
"outer_radius") *
cm)
87 , m_SizeZ(m_Params->get_double_param(
"size_z") *
cm)
88 , m_ScintiTileZ(m_SizeZ)
89 , m_ScintiTileThickness(m_Params->get_double_param(
"scinti_tile_thickness") *
cm)
90 , m_ScintiInnerGap(m_Params->get_double_param(
"scinti_inner_gap") *
cm)
91 , m_ScintiOuterGap(m_Params->get_double_param(
"scinti_outer_gap") *
cm)
92 , m_ScintiOuterRadius(m_Params->get_double_param(
"scinti_outer_radius") *
cm)
93 , m_TiltAngle(m_Params->get_double_param(
"tilt_angle") *
deg)
94 , m_EnvelopeInnerRadius(m_InnerRadius)
95 , m_EnvelopeOuterRadius(m_OuterRadius)
96 , m_EnvelopeZ(m_SizeZ)
97 , m_NumScintiPlates(m_Params->get_int_param(PHG4HcalDefs::
scipertwr) * m_Params->get_int_param(
"n_towers"))
98 , m_NumScintiTilesPos(m_Params->get_int_param(
"n_scinti_tiles_pos"))
99 , m_NumScintiTilesNeg(m_Params->get_int_param(
"n_scinti_tiles_neg"))
100 , m_Active(m_Params->get_int_param(
"active"))
101 , m_AbsorberActive(m_Params->get_int_param(
"absorberactive"))
102 , m_ScintiLogicNamePrefix(
"HcalInnerScinti")
159 Point_2 p_upperedge(xcoord, ycoord);
160 Line_2 s2(p_in_1, p_upperedge);
164 Circle_2 outer_circle(sc1, sc2, sc3);
165 std::vector<CGAL::Object> res;
166 CGAL::intersection(outer_circle, perp, std::back_inserter(res));
168 std::vector<CGAL::Object>::const_iterator iter;
169 for (iter = res.begin(); iter != res.end(); ++iter)
171 CGAL::Object obj = *iter;
172 if (
const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>,
unsigned> *
point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>,
unsigned>>(&obj))
174 if (CGAL::to_double(
point->first.x()) > CGAL::to_double(p_upperedge.x()))
176 double deltax = CGAL::to_double(
point->first.x()) - CGAL::to_double(p_upperedge.x());
177 double deltay = CGAL::to_double(
point->first.y()) - CGAL::to_double(p_upperedge.y());
180 Point_2 pntmp(CGAL::to_double(
point->first.x()), CGAL::to_double(
point->first.y()));
186 std::cout <<
"CGAL::Object type not pair..." << std::endl;
192 Point_2 p_loweredge(xcoord, ycoord);
193 Line_2 s3(p_in_1, p_loweredge);
194 Line_2 l_lower = s3.perpendicular(p_loweredge);
196 Circle_2 inner_circle(ic1, ic2, ic3);
198 CGAL::intersection(inner_circle, l_lower, std::back_inserter(res));
203 for (iter = res.begin(); iter != res.end(); ++iter)
205 CGAL::Object obj = *iter;
206 if (
const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>,
unsigned> *
point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>,
unsigned>>(&obj))
208 if (CGAL::to_double(
point->first.x()) > minx)
210 minx = CGAL::to_double(
point->first.x());
211 double deltax = CGAL::to_double(
point->first.x()) - CGAL::to_double(p_loweredge.x());
212 double deltay = CGAL::to_double(
point->first.y()) - CGAL::to_double(p_loweredge.y());
214 Point_2 pntmp(CGAL::to_double(
point->first.x()), CGAL::to_double(
point->first.y()));
238 Point_2 p_loweredge(xcoord, ycoord);
239 Line_2 s2(p_in_1, p_loweredge);
242 Circle_2 inner_circle(sc1, sc2, sc3);
243 std::vector<CGAL::Object> res;
244 CGAL::intersection(inner_circle, perp, std::back_inserter(res));
246 std::vector<CGAL::Object>::const_iterator iter;
247 for (iter = res.begin(); iter != res.end(); ++iter)
249 CGAL::Object obj = *iter;
250 if (
const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>,
unsigned> *
point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>,
unsigned>>(&obj))
252 if (CGAL::to_double(
point->first.x()) > 0)
254 Point_2 pntmp(CGAL::to_double(
point->first.x()), CGAL::to_double(
point->first.y()));
260 std::cout <<
"CGAL::Object type not pair..." << std::endl;
266 Point_2 p_loweredge2(xcoord2, ycoord2);
267 Line_2 s2_2(p_in_1, p_loweredge2);
268 Line_2 perp2 = s2_2.perpendicular(p_loweredge2);
271 Circle_2 outer_circle(so1, so2, so3);
273 CGAL::intersection(outer_circle, perp2, std::back_inserter(res));
275 for (iter = res.begin(); iter != res.end(); ++iter)
277 CGAL::Object obj = *iter;
278 if (
const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>,
unsigned> *
point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>,
unsigned>>(&obj))
280 if (CGAL::to_double(
point->first.x()) > CGAL::to_double(p_loweredge.x()))
282 Point_2 pntmp(CGAL::to_double(
point->first.x()), CGAL::to_double(
point->first.y()));
288 std::cout <<
"CGAL::Object type not pair..." << std::endl;
295 double xmidpoint = cos(phi_midpoint) * mid_radius;
296 double ymidpoint = sin(phi_midpoint) * mid_radius;
298 angle_mid_scinti = (M_PI / 2. - phi_midpoint) - (M_PI / 2. +
m_TiltAngle /
rad);
303 Point_2 mid_upperscint(xmidpoint, ymidpoint);
304 Point_2 p_upperedge(xcoordup, ycoordup);
305 Line_2 sup(mid_upperscint, p_upperedge);
306 Line_2 perpA = sup.perpendicular(p_upperedge);
308 Circle_2 inner_circleA(sc1A, sc2A, sc3A);
309 std::vector<CGAL::Object> resA;
310 CGAL::intersection(inner_circleA, perpA, std::back_inserter(resA));
311 std::vector<CGAL::Object>::const_iterator iterA;
313 for (iterA = resA.begin(); iterA != resA.end(); ++iterA)
315 CGAL::Object obj = *iterA;
316 if (
const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>,
unsigned> *
point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>,
unsigned>>(&obj))
318 if (CGAL::to_double(
point->first.x()) > pxmax)
320 pxmax = CGAL::to_double(
point->first.x());
321 Point_2 pntmp(CGAL::to_double(
point->first.x()), CGAL::to_double(
point->first.y()));
327 std::cout <<
"CGAL::Object type not pair..." << std::endl;
333 Point_2 p_upperedge2(xcoordup2, ycoordup2);
334 Line_2 sup2(mid_upperscint, p_upperedge2);
335 Line_2 perpA2 = sup2.perpendicular(p_upperedge2);
338 Circle_2 outer_circleA(so1A, so2A, so3A);
340 CGAL::intersection(outer_circleA, perpA2, std::back_inserter(resA));
341 for (iterA = resA.begin(); iterA != resA.end(); ++iterA)
343 CGAL::Object obj = *iterA;
344 if (
const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>,
unsigned> *
point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>,
unsigned>>(&obj))
346 if (CGAL::to_double(
point->first.x()) > CGAL::to_double(p_loweredge.x()))
348 Point_2 pntmp(CGAL::to_double(
point->first.x()), CGAL::to_double(
point->first.y()));
354 std::cout <<
"CGAL::Object type not pair..." << std::endl;
360 G4TwoVector
v1(CGAL::to_double(upperleft.x()), CGAL::to_double(upperleft.y()));
361 G4TwoVector
v2(CGAL::to_double(upperright.x()), CGAL::to_double(upperright.y()));
362 G4TwoVector
v3(CGAL::to_double(lowerright.x()), CGAL::to_double(lowerright.y()));
363 G4TwoVector
v4(CGAL::to_double(lowerleft.x()), CGAL::to_double(lowerleft.y()));
364 std::vector<G4TwoVector> vertexes;
365 vertexes.push_back(v1);
366 vertexes.push_back(v2);
367 vertexes.push_back(v3);
368 vertexes.push_back(v4);
369 G4TwoVector
zero(0, 0);
370 G4VSolid *steel_plate =
new G4ExtrudedSolid(
"SteelPlate",
383 Line_2 secant(lowleft, upleft);
386 double xmid = (CGAL::to_double(lowleft.x()) + CGAL::to_double(upleft.x())) / 2.;
387 double ymid = (CGAL::to_double(lowleft.y()) + CGAL::to_double(upleft.y())) / 2.;
389 Line_2 sekperp = secant.perpendicular(midpoint);
391 Circle_2 inner_circle(sc1, sc2, sc3);
392 std::vector<CGAL::Object> res;
393 CGAL::intersection(inner_circle, sekperp, std::back_inserter(res));
394 std::vector<CGAL::Object>::const_iterator iter;
397 for (iter = res.begin(); iter != res.end(); ++iter)
399 CGAL::Object obj = *iter;
400 if (
const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>,
unsigned> *
point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>,
unsigned>>(&obj))
402 if (CGAL::to_double(
point->first.x()) > pxmax)
404 pxmax = CGAL::to_double(
point->first.x());
405 Point_2 pntmp(CGAL::to_double(
point->first.x()), CGAL::to_double(
point->first.y()));
411 std::cout <<
"CGAL::Object type not pair..." << std::endl;
414 Line_2 leftside = sekperp.perpendicular(tangtouch);
415 CGAL::Object result = CGAL::intersection(upedge, leftside);
416 if (
const Point_2 *ipoint = CGAL::object_cast<Point_2>(&result))
420 result = CGAL::intersection(lowedge, leftside);
421 if (
const Point_2 *ipoint = CGAL::object_cast<Point_2>(&result))
436 G4LogicalVolume *hcal_envelope_log =
new G4LogicalVolume(hcal_envelope_cylinder, Air, G4String(
"Hcal_envelope"),
nullptr,
nullptr,
nullptr);
438 G4RotationMatrix hcal_rotm;
466 boost::char_separator<char> sep(
"_");
467 boost::tokenizer<boost::char_separator<char>> tok((*it)->GetName(), sep);
468 boost::tokenizer<boost::char_separator<char>>::const_iterator tokeniter;
469 for (tokeniter = tok.begin(); tokeniter != tok.end(); ++tokeniter)
471 if (*tokeniter ==
"impr")
474 if (tokeniter != tok.end())
476 int layer_id = boost::lexical_cast<
int>(*tokeniter);
481 int tower_id = (*it)->GetCopyNo() - layer_id;
483 std::pair<int, int> layer_twr = std::make_pair(layer_id, tower_id);
487 std::cout <<
"invalid scintillator row " << layer_id
494 std::cout <<
PHWHERE <<
" Error parsing " << (*it)->GetName()
495 <<
" for mother volume number " << std::endl;
518 std::ostringstream
name;
533 double xp = cos(phi) * middlerad;
534 double yp = sin(phi) * middlerad;
541 phi = -atan(yo / xo);
547 phi = -atan(yo / xo);
552 G4RotationMatrix *Rot =
new G4RotationMatrix();
553 double ypos = sin(phi) * middlerad;
554 double xpos = cos(phi) * middlerad;
561 G4ThreeVector g4vec(xpos, ypos, 0);
569 Rot =
new G4RotationMatrix();
570 Rot->rotateZ(-phi *
rad);
572 name <<
"InnerHcalSteel_" <<
i;
601 std::ostringstream
name;
617 eta += delta_eta_pos;
624 z[0] += scinti_gap_neighbor / 2.;
625 z[1] += scinti_gap_neighbor / 2.;
626 z[2] -= scinti_gap_neighbor / 2.;
627 z[3] -= scinti_gap_neighbor / 2.;
628 Point_2 leftsidelow(z[0], x[0]);
629 Point_2 leftsidehigh(z[1], x[1]);
631 z[0] =
x_at_y(leftsidelow, leftsidehigh, x[0]);
633 z[1] =
x_at_y(leftsidelow, leftsidehigh, x[1]);
634 Point_2 rightsidelow(z[2], x[2]);
635 Point_2 rightsidehigh(z[3], x[3]);
637 z[2] =
x_at_y(rightsidelow, rightsidehigh, x[2]);
639 z[3] =
x_at_y(rightsidelow, rightsidehigh, x[3]);
641 std::vector<G4TwoVector> vertexes;
642 for (
int j = 0;
j < 4;
j++)
644 G4TwoVector
v(x[
j], z[j]);
645 vertexes.push_back(v);
647 G4TwoVector
zero(0, 0);
649 G4VSolid *scinti =
new G4ExtrudedSolid(
"ScintillatorTile",
654 G4RotationMatrix *rotm =
new G4RotationMatrix();
655 rotm->rotateX(-90 *
deg);
657 name <<
"scintillator_" <<
i <<
"_left";
671 eta += delta_eta_neg;
678 z[0] += scinti_gap_neighbor / 2.;
679 z[1] += scinti_gap_neighbor / 2.;
680 z[2] -= scinti_gap_neighbor / 2.;
681 z[3] -= scinti_gap_neighbor / 2.;
682 Point_2 leftsidelow(z[0], x[0]);
683 Point_2 leftsidehigh(z[1], x[1]);
685 z[0] =
x_at_y(leftsidelow, leftsidehigh, x[0]);
687 z[1] =
x_at_y(leftsidelow, leftsidehigh, x[1]);
688 Point_2 rightsidelow(z[2], x[2]);
689 Point_2 rightsidehigh(z[3], x[3]);
691 z[2] =
x_at_y(rightsidelow, rightsidehigh, x[2]);
693 z[3] =
x_at_y(rightsidelow, rightsidehigh, x[3]);
695 std::vector<G4TwoVector> vertexes;
696 for (
int j = 0;
j < 4;
j++)
698 G4TwoVector
v(x[
j], z[j]);
699 vertexes.push_back(v);
701 G4TwoVector
zero(0, 0);
703 G4VSolid *scinti =
new G4ExtrudedSolid(
"ScintillatorTile",
709 G4RotationMatrix *rotm =
new G4RotationMatrix();
710 rotm->rotateX(90 *
deg);
712 name <<
"scintillator_" <<
i <<
"_right";
735 x[0] = CGAL::to_double(p0.x());
736 y[0] = CGAL::to_double(p0.y());
737 x[1] = CGAL::to_double(p1.x());
738 y[1] = CGAL::to_double(p1.y());
740 double newx = fabs(x[0]) + fabs(x[1]);
744 CGAL::Object result = CGAL::intersection(l, seg);
745 if (
const Point_2 *ipoint = CGAL::object_cast<Point_2>(&result))
747 xret = CGAL::to_double(ipoint->x());
751 std::cout <<
PHWHERE <<
" failed for y = " << y << std::endl;
752 std::cout <<
"p0(x): " << CGAL::to_double(p0.x()) <<
", p0(y): " << CGAL::to_double(p0.y()) << std::endl;
753 std::cout <<
"p1(x): " << CGAL::to_double(p1.x()) <<
", p1(y): " << CGAL::to_double(p1.y()) << std::endl;
763 G4AssemblyVolume *assmeblyvol =
new G4AssemblyVolume();
764 std::ostringstream
name;
772 G4UserLimits *g4userlimits =
nullptr;
773 if (isfinite(steplimits))
775 g4userlimits =
new G4UserLimits(steplimits);
779 assmeblyvol->AddPlacedVolume(scinti_tile_logic, g4vec,
nullptr);
798 Point_2 pxnull(xcoord, ycoord);
801 Circle_2 inner_circle(sc1, sc2, sc3);
802 std::vector<CGAL::Object> res;
803 CGAL::intersection(inner_circle, s2, std::back_inserter(res));
807 <<
" too large, no intersection with inner radius" << std::endl;
820 <<
" cm" << std::endl;
827 <<
" cm" << std::endl;
834 <<
" cm" << std::endl;
854 std::cout <<
"both number of crossings and tilt angle are set" << std::endl;
855 std::cout <<
"using number of crossings to determine tilt angle" << std::endl;
861 Point_2 phightmp(1, tan(deltaphi));
863 Circle_2 inner_circle(pin1, pin2, pin3);
864 Point_2 pmid1(mid_radius, 0), pmid2(0, mid_radius), pmid3(-mid_radius, 0);
865 Circle_2 mid_circle(pmid1, pmid2, pmid3);
867 Circle_2 outer_circle(pout1, pout2, pout3);
868 Line_2 l_up(pnull, phightmp);
869 std::vector<CGAL::Object> res;
870 CGAL::intersection(outer_circle, l_up, std::back_inserter(res));
872 std::vector<CGAL::Object>::const_iterator iter;
873 for (iter = res.begin(); iter != res.end(); ++iter)
875 CGAL::Object obj = *iter;
876 if (
const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>,
unsigned> *
point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>,
unsigned>>(&obj))
878 if (CGAL::to_double(
point->first.x()) > 0)
880 Point_2 pntmp(CGAL::to_double(
point->first.x()), CGAL::to_double(
point->first.y()));
886 std::cout <<
"CGAL::Object type not pair..." << std::endl;
890 Line_2 l_right(plow, upperright);
893 CGAL::intersection(mid_circle, l_right, std::back_inserter(res));
894 for (iter = res.begin(); iter != res.end(); ++iter)
896 CGAL::Object obj = *iter;
897 if (
const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>,
unsigned> *
point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>,
unsigned>>(&obj))
899 if (CGAL::to_double(
point->first.x()) > 0)
901 Point_2 pntmp(CGAL::to_double(
point->first.x()), CGAL::to_double(
point->first.y()));
907 std::cout <<
"CGAL::Object type not pair..." << std::endl;
912 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()));
913 double upside = sqrt(CGAL::to_double(midpoint.x()) * CGAL::to_double(midpoint.x()) + CGAL::to_double(midpoint.y()) * CGAL::to_double(midpoint.y()));
917 tiltangle = tiltangle *
rad;
925 std::cout <<
"Inner Hcal Detector:" << std::endl;
926 if (what ==
"ALL" || what ==
"VOLUME")
929 std::cout <<
"Volume Steel: " <<
m_VolumeSteel /
cm3 <<
" cm^3" << std::endl;
943 std::cout <<
"could not locate volume " << volume->GetName()
944 <<
" in Inner Hcal scintillator map" << std::endl;
958 std::cout <<
PHWHERE <<
"Run Node missing, exiting." << std::endl;
983 double geom_ref_radius = innerrad + thickness / 2.;
985 if (!std::isfinite(phistart))
987 std::cout <<
PHWHERE <<
" phistart is not finite: " << phistart
988 <<
", exiting now (this will crash anyway)" << std::endl;
994 std::pair<double, double> range = std::make_pair(phiend, phistart);
1002 double etahibound = etalowbound +
1004 std::pair<double, double> range = std::make_pair(etalowbound, etahibound);
1006 etalowbound = etahibound;
1023 std::cout <<
"IHCalDetector::InitRun - Tower geometry " << key <<
" already exists" << std::endl;
1026 if (fabs(tg->get_center_x() -
x) > 1
e-4)
1028 std::cout <<
"IHCalDetector::InitRun - Fatal Error - duplicated Tower geometry " << key <<
" with existing x = " << tg->
get_center_x() <<
" and expected x = " << x
1033 if (fabs(tg->get_center_y() -
y) > 1
e-4)
1035 std::cout <<
"IHCalDetector::InitRun - Fatal Error - duplicated Tower geometry " << key <<
" with existing y = " << tg->get_center_y() <<
" and expected y = " << y
1039 if (fabs(tg->get_center_z() -
z) > 1
e-4)
1041 std::cout <<
"IHCalDetector::InitRun - Fatal Error - duplicated Tower geometry " << key <<
" with existing z= " << tg->get_center_z() <<
" and expected z = " <<
z
1050 std::cout <<
"IHCalDetector::InitRun - building tower geometry " << key <<
"" << std::endl;
1055 tg->set_center_x(x);
1056 tg->set_center_y(y);
1057 tg->set_center_z(
z);