7 #include <phparameter/PHParameters.h>
19 #include <Geant4/G4AssemblyVolume.hh>
20 #include <Geant4/G4Box.hh>
21 #include <Geant4/G4ExtrudedSolid.hh>
22 #include <Geant4/G4IntersectionSolid.hh>
23 #include <Geant4/G4LogicalVolume.hh>
24 #include <Geant4/G4Material.hh>
25 #include <Geant4/G4PVPlacement.hh>
26 #include <Geant4/G4RotationMatrix.hh>
27 #include <Geant4/G4String.hh>
28 #include <Geant4/G4SystemOfUnits.hh>
29 #include <Geant4/G4ThreeVector.hh>
30 #include <Geant4/G4Transform3D.hh>
31 #include <Geant4/G4Tubs.hh>
32 #include <Geant4/G4TwoVector.hh>
33 #include <Geant4/G4Types.hh>
34 #include <Geant4/G4UserLimits.hh>
35 #include <Geant4/G4VSolid.hh>
37 #include <CGAL/Boolean_set_operations_2.h>
38 #include <CGAL/Circular_kernel_intersections.h>
39 #include <CGAL/Exact_circular_kernel_2.h>
40 #include <CGAL/Object.h>
41 #include <CGAL/point_generators_2.h>
43 #include <boost/lexical_cast.hpp>
44 #include <boost/tokenizer.hpp>
55 typedef CGAL::Circle_2<PHG4InnerHcalDetector::Circular_k>
Circle_2;
57 typedef CGAL::Line_2<PHG4InnerHcalDetector::Circular_k>
Line_2;
58 typedef CGAL::Segment_2<PHG4InnerHcalDetector::Circular_k>
Segment_2;
70 , m_Params(parameters)
71 , m_ScintiMotherAssembly(nullptr)
72 , m_InnerRadius(m_Params->get_double_param(
"inner_radius") *
cm)
73 , m_OuterRadius(m_Params->get_double_param(
"outer_radius") *
cm)
74 , m_SizeZ(m_Params->get_double_param(
"size_z") *
cm)
76 , m_ScintiTileXLower(NAN)
77 , m_ScintiTileXUpper(NAN)
78 , m_ScintiTileZ(m_SizeZ)
79 , m_ScintiTileThickness(m_Params->get_double_param(
"scinti_tile_thickness") *
cm)
80 , m_ScintiInnerGap(m_Params->get_double_param(
"scinti_inner_gap") *
cm)
81 , m_ScintiOuterGap(m_Params->get_double_param(
"scinti_outer_gap") *
cm)
82 , m_ScintiOuterRadius(m_Params->get_double_param(
"scinti_outer_radius") *
cm)
83 , m_TiltAngle(m_Params->get_double_param(
"tilt_angle") *
deg)
84 , m_EnvelopeInnerRadius(m_InnerRadius)
85 , m_EnvelopeOuterRadius(m_OuterRadius)
86 , m_EnvelopeZ(m_SizeZ)
87 , m_VolumeEnvelope(NAN)
89 , m_VolumeScintillator(NAN)
90 , m_NumScintiPlates(m_Params->get_int_param(PHG4HcalDefs::
scipertwr) * m_Params->get_int_param(
"n_towers"))
91 , m_NumScintiTilesPos(m_Params->get_int_param(
"n_scinti_tiles_pos"))
92 , m_NumScintiTilesNeg(m_Params->get_int_param(
"n_scinti_tiles_neg"))
93 , m_Active(m_Params->get_int_param(
"active"))
94 , m_AbsorberActive(m_Params->get_int_param(
"absorberactive"))
96 , m_ScintiLogicNamePrefix(
"HcalInnerScinti")
153 Point_2 p_upperedge(xcoord, ycoord);
154 Line_2 s2(p_in_1, p_upperedge);
158 Circle_2 outer_circle(sc1, sc2, sc3);
159 vector<CGAL::Object> res;
160 CGAL::intersection(outer_circle, perp, std::back_inserter(res));
162 vector<CGAL::Object>::const_iterator iter;
163 for (iter = res.begin(); iter != res.end(); ++iter)
165 CGAL::Object obj = *iter;
166 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))
168 if (CGAL::to_double(
point->first.x()) > CGAL::to_double(p_upperedge.x()))
170 double deltax = CGAL::to_double(
point->first.x()) - CGAL::to_double(p_upperedge.x());
171 double deltay = CGAL::to_double(
point->first.y()) - CGAL::to_double(p_upperedge.y());
174 Point_2 pntmp(CGAL::to_double(
point->first.x()), CGAL::to_double(
point->first.y()));
180 cout <<
"CGAL::Object type not pair..." << endl;
186 Point_2 p_loweredge(xcoord, ycoord);
187 Line_2 s3(p_in_1, p_loweredge);
188 Line_2 l_lower = s3.perpendicular(p_loweredge);
190 Circle_2 inner_circle(ic1, ic2, ic3);
192 CGAL::intersection(inner_circle, l_lower, std::back_inserter(res));
197 for (iter = res.begin(); iter != res.end(); ++iter)
199 CGAL::Object obj = *iter;
200 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))
202 if (CGAL::to_double(
point->first.x()) > minx)
204 minx = CGAL::to_double(
point->first.x());
205 double deltax = CGAL::to_double(
point->first.x()) - CGAL::to_double(p_loweredge.x());
206 double deltay = CGAL::to_double(
point->first.y()) - CGAL::to_double(p_loweredge.y());
208 Point_2 pntmp(CGAL::to_double(
point->first.x()), CGAL::to_double(
point->first.y()));
232 Point_2 p_loweredge(xcoord, ycoord);
233 Line_2 s2(p_in_1, p_loweredge);
236 Circle_2 inner_circle(sc1, sc2, sc3);
237 vector<CGAL::Object> res;
238 CGAL::intersection(inner_circle, perp, std::back_inserter(res));
240 vector<CGAL::Object>::const_iterator iter;
241 for (iter = res.begin(); iter != res.end(); ++iter)
243 CGAL::Object obj = *iter;
244 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))
246 if (CGAL::to_double(
point->first.x()) > 0)
248 Point_2 pntmp(CGAL::to_double(
point->first.x()), CGAL::to_double(
point->first.y()));
254 cout <<
"CGAL::Object type not pair..." << endl;
260 Point_2 p_loweredge2(xcoord2, ycoord2);
261 Line_2 s2_2(p_in_1, p_loweredge2);
262 Line_2 perp2 = s2_2.perpendicular(p_loweredge2);
265 Circle_2 outer_circle(so1, so2, so3);
267 CGAL::intersection(outer_circle, perp2, std::back_inserter(res));
269 for (iter = res.begin(); iter != res.end(); ++iter)
271 CGAL::Object obj = *iter;
272 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))
274 if (CGAL::to_double(
point->first.x()) > CGAL::to_double(p_loweredge.x()))
276 Point_2 pntmp(CGAL::to_double(
point->first.x()), CGAL::to_double(
point->first.y()));
282 cout <<
"CGAL::Object type not pair..." << endl;
289 double xmidpoint = cos(phi_midpoint) * mid_radius;
290 double ymidpoint = sin(phi_midpoint) * mid_radius;
297 Point_2 mid_upperscint(xmidpoint, ymidpoint);
298 Point_2 p_upperedge(xcoordup, ycoordup);
300 Line_2 sup(mid_upperscint, p_upperedge);
301 Line_2 perp = sup.perpendicular(p_upperedge);
303 Circle_2 inner_circle(sc1, sc2, sc3);
304 vector<CGAL::Object> res;
305 CGAL::intersection(inner_circle, perp, std::back_inserter(res));
306 vector<CGAL::Object>::const_iterator iter;
308 for (iter = res.begin(); iter != res.end(); ++iter)
310 CGAL::Object obj = *iter;
311 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))
313 if (CGAL::to_double(
point->first.x()) > pxmax)
315 pxmax = CGAL::to_double(
point->first.x());
316 Point_2 pntmp(CGAL::to_double(
point->first.x()), CGAL::to_double(
point->first.y()));
322 cout <<
"CGAL::Object type not pair..." << endl;
328 Point_2 p_upperedge2(xcoordup2, ycoordup2);
329 Line_2 sup2(mid_upperscint, p_upperedge2);
330 Line_2 perp2 = sup2.perpendicular(p_upperedge2);
333 Circle_2 outer_circle(so1, so2, so3);
335 CGAL::intersection(outer_circle, perp2, std::back_inserter(res));
336 for (iter = res.begin(); iter != res.end(); ++iter)
338 CGAL::Object obj = *iter;
339 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))
341 if (CGAL::to_double(
point->first.x()) > CGAL::to_double(p_loweredge.x()))
343 Point_2 pntmp(CGAL::to_double(
point->first.x()), CGAL::to_double(
point->first.y()));
349 cout <<
"CGAL::Object type not pair..." << endl;
356 G4TwoVector v1(CGAL::to_double(upperleft.x()), CGAL::to_double(upperleft.y()));
357 G4TwoVector v2(CGAL::to_double(upperright.x()), CGAL::to_double(upperright.y()));
358 G4TwoVector v3(CGAL::to_double(lowerright.x()), CGAL::to_double(lowerright.y()));
359 G4TwoVector v4(CGAL::to_double(lowerleft.x()), CGAL::to_double(lowerleft.y()));
360 std::vector<G4TwoVector> vertexes;
361 vertexes.push_back(v1);
362 vertexes.push_back(v2);
363 vertexes.push_back(v3);
364 vertexes.push_back(v4);
379 Line_2 secant(lowleft, upleft);
382 double xmid = (CGAL::to_double(lowleft.x()) + CGAL::to_double(upleft.x())) / 2.;
383 double ymid = (CGAL::to_double(lowleft.y()) + CGAL::to_double(upleft.y())) / 2.;
385 Line_2 sekperp = secant.perpendicular(midpoint);
387 Circle_2 inner_circle(sc1, sc2, sc3);
388 vector<CGAL::Object> res;
389 CGAL::intersection(inner_circle, sekperp, std::back_inserter(res));
390 vector<CGAL::Object>::const_iterator iter;
393 for (iter = res.begin(); iter != res.end(); ++iter)
395 CGAL::Object obj = *iter;
396 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))
398 if (CGAL::to_double(
point->first.x()) > pxmax)
400 pxmax = CGAL::to_double(
point->first.x());
401 Point_2 pntmp(CGAL::to_double(
point->first.x()), CGAL::to_double(
point->first.y()));
407 cout <<
"CGAL::Object type not pair..." << endl;
410 Line_2 leftside = sekperp.perpendicular(tangtouch);
411 CGAL::Object result = CGAL::intersection(upedge, leftside);
412 if (
const Point_2 *ipoint = CGAL::object_cast<Point_2>(&result))
416 result = CGAL::intersection(lowedge, leftside);
417 if (
const Point_2 *ipoint = CGAL::object_cast<Point_2>(&result))
462 boost::char_separator<char> sep(
"_");
463 boost::tokenizer<boost::char_separator<char>> tok((*it)->GetName(), sep);
464 boost::tokenizer<boost::char_separator<char>>::const_iterator tokeniter;
465 for (tokeniter = tok.begin(); tokeniter != tok.end(); ++tokeniter)
467 if (*tokeniter ==
"impr")
470 if (tokeniter != tok.end())
472 int layer_id = boost::lexical_cast<
int>(*tokeniter);
477 int tower_id = (*it)->GetCopyNo() - layer_id;
479 pair<int, int> layer_twr = make_pair(layer_id, tower_id);
483 cout <<
"invalid scintillator row " << layer_id
490 cout <<
PHWHERE <<
" Error parsing " << (*it)->GetName()
491 <<
" for mother volume number " << endl;
528 double xp = cos(phi) * middlerad;
529 double yp = sin(phi) * middlerad;
536 phi = -atan(yo / xo);
542 phi = -atan(yo / xo);
548 double ypos = sin(phi) * middlerad;
549 double xpos = cos(phi) * middlerad;
567 name <<
"InnerHcalSteel_" << i;
611 eta += delta_eta_pos;
618 z[0] += scinti_gap_neighbor / 2.;
619 z[1] += scinti_gap_neighbor / 2.;
620 z[2] -= scinti_gap_neighbor / 2.;
621 z[3] -= scinti_gap_neighbor / 2.;
622 Point_2 leftsidelow(z[0], x[0]);
623 Point_2 leftsidehigh(z[1], x[1]);
625 z[0] =
x_at_y(leftsidelow, leftsidehigh, x[0]);
627 z[1] =
x_at_y(leftsidelow, leftsidehigh, x[1]);
628 Point_2 rightsidelow(z[2], x[2]);
629 Point_2 rightsidehigh(z[3], x[3]);
631 z[2] =
x_at_y(rightsidelow, rightsidehigh, x[2]);
633 z[3] =
x_at_y(rightsidelow, rightsidehigh, x[3]);
635 vector<G4TwoVector> vertexes;
636 for (
int j = 0; j < 4; j++)
639 vertexes.push_back(v);
651 name <<
"scintillator_" << i <<
"_left";
666 eta += delta_eta_neg;
673 z[0] += scinti_gap_neighbor / 2.;
674 z[1] += scinti_gap_neighbor / 2.;
675 z[2] -= scinti_gap_neighbor / 2.;
676 z[3] -= scinti_gap_neighbor / 2.;
677 Point_2 leftsidelow(z[0], x[0]);
678 Point_2 leftsidehigh(z[1], x[1]);
680 z[0] =
x_at_y(leftsidelow, leftsidehigh, x[0]);
682 z[1] =
x_at_y(leftsidelow, leftsidehigh, x[1]);
683 Point_2 rightsidelow(z[2], x[2]);
684 Point_2 rightsidehigh(z[3], x[3]);
686 z[2] =
x_at_y(rightsidelow, rightsidehigh, x[2]);
688 z[3] =
x_at_y(rightsidelow, rightsidehigh, x[3]);
690 vector<G4TwoVector> vertexes;
691 for (
int j = 0; j < 4; j++)
694 vertexes.push_back(v);
707 name <<
"scintillator_" << i <<
"_right";
731 x[0] = CGAL::to_double(p0.x());
732 y[0] = CGAL::to_double(p0.y());
733 x[1] = CGAL::to_double(p1.x());
734 y[1] = CGAL::to_double(p1.y());
736 double newx = fabs(x[0]) + fabs(x[1]);
740 CGAL::Object result = CGAL::intersection(l, s);
741 if (
const Point_2 *ipoint = CGAL::object_cast<Point_2>(&result))
743 xret = CGAL::to_double(ipoint->x());
747 cout <<
PHWHERE <<
" failed for y = " << y << endl;
748 cout <<
"p0(x): " << CGAL::to_double(p0.x()) <<
", p0(y): " << CGAL::to_double(p0.y()) << endl;
749 cout <<
"p1(x): " << CGAL::to_double(p1.x()) <<
", p1(y): " << CGAL::to_double(p1.y()) << endl;
794 Point_2 pxnull(xcoord, ycoord);
797 Circle_2 inner_circle(sc1, sc2, sc3);
798 vector<CGAL::Object> res;
799 CGAL::intersection(inner_circle, s2, std::back_inserter(res));
803 <<
" too large, no intersection with inner radius" << endl;
850 cout <<
"both number of crossings and tilt angle are set" << endl;
851 cout <<
"using number of crossings to determine tilt angle" << endl;
857 Point_2 phightmp(1, tan(deltaphi));
859 Circle_2 inner_circle(pin1, pin2, pin3);
860 Point_2 pmid1(mid_radius, 0), pmid2(0, mid_radius), pmid3(-mid_radius, 0);
861 Circle_2 mid_circle(pmid1, pmid2, pmid3);
863 Circle_2 outer_circle(pout1, pout2, pout3);
864 Line_2 l_up(pnull, phightmp);
865 vector<CGAL::Object> res;
866 CGAL::intersection(outer_circle, l_up, std::back_inserter(res));
868 vector<CGAL::Object>::const_iterator iter;
869 for (iter = res.begin(); iter != res.end(); ++iter)
871 CGAL::Object obj = *iter;
872 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))
874 if (CGAL::to_double(
point->first.x()) > 0)
876 Point_2 pntmp(CGAL::to_double(
point->first.x()), CGAL::to_double(
point->first.y()));
882 cout <<
"CGAL::Object type not pair..." << endl;
886 Line_2 l_right(plow, upperright);
889 CGAL::intersection(mid_circle, l_right, std::back_inserter(res));
890 for (iter = res.begin(); iter != res.end(); ++iter)
892 CGAL::Object obj = *iter;
893 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))
895 if (CGAL::to_double(
point->first.x()) > 0)
897 Point_2 pntmp(CGAL::to_double(
point->first.x()), CGAL::to_double(
point->first.y()));
903 cout <<
"CGAL::Object type not pair..." << endl;
908 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()));
909 double upside = sqrt(CGAL::to_double(midpoint.x()) * CGAL::to_double(midpoint.x()) + CGAL::to_double(midpoint.y()) * CGAL::to_double(midpoint.y()));
913 tiltangle = tiltangle *
rad;
921 cout <<
"Inner Hcal Detector:" << endl;
922 if (what ==
"ALL" || what ==
"VOLUME")
939 cout <<
"could not locate volume " << volume->
GetName()
940 <<
" in Inner Hcal scintillator map" << endl;