ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4InnerHcalDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4InnerHcalDetector.cc
2 
3 #include "PHG4HcalDefs.h"
6 
7 #include <phparameter/PHParameters.h>
8 
9 #include <g4main/PHG4Detector.h>
11 #include <g4main/PHG4Subsystem.h>
12 #include <g4main/PHG4Utils.h>
13 
14 #include <phool/phool.h>
15 #include <phool/recoConsts.h>
16 
17 #include <TSystem.h>
18 
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>
36 
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>
42 
43 #include <boost/lexical_cast.hpp>
44 #include <boost/tokenizer.hpp>
45 
46 #include <algorithm>
47 #include <cmath>
48 #include <cstdlib>
49 #include <iostream>
50 #include <iterator>
51 #include <sstream>
52 
54 
55 typedef CGAL::Circle_2<PHG4InnerHcalDetector::Circular_k> Circle_2;
56 typedef CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k> Circular_arc_point_2;
57 typedef CGAL::Line_2<PHG4InnerHcalDetector::Circular_k> Line_2;
58 typedef CGAL::Segment_2<PHG4InnerHcalDetector::Circular_k> Segment_2;
59 
60 using namespace std;
61 
62 // there is still a minute problem for very low tilt angles where the scintillator
63 // face touches the boundary instead of the corner, subtracting 1 permille from the total
64 // scintilator length takes care of this
65 static double subtract_from_scinti_x = 0.1 * mm;
66 
68  : PHG4Detector(subsys, Node, dnam)
69  , m_DisplayAction(dynamic_cast<PHG4InnerHcalDisplayAction *>(subsys->GetDisplayAction()))
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)
75  , m_ScintiTileX(NAN)
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)
88  , m_VolumeSteel(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"))
95  , m_Layer(0)
96  , m_ScintiLogicNamePrefix("HcalInnerScinti")
97 {
98 
99  // n_scinti_tiles takes precedence
100 
102  if(nTiles <= 0){
104  }
105  else{
106  m_NumScintiTilesPos = nTiles/2;
107  m_Params->set_int_param("n_scinti_tiles_pos", nTiles/2);
108  m_NumScintiTilesNeg = nTiles/2;
109  m_Params->set_int_param("n_scinti_tiles_neg", nTiles/2);
110  }
111 
112  // allocate memory for scintillator plates
113  m_ScintiTilesVec.assign(nTiles, static_cast<G4VSolid *>(nullptr));
114 
115 }
116 
118 {
119  delete m_ScintiMotherAssembly;
120 }
121 
122 //_______________________________________________________________
123 //_______________________________________________________________
125 {
126  if (m_AbsorberActive)
127  {
128  if (m_SteelAbsorberPhysVolSet.find(volume) != m_SteelAbsorberPhysVolSet.end())
129  {
130  return -1;
131  }
132  }
133  if (m_Active)
134  {
135  if (m_ScintiTilePhysVolMap.find(volume) != m_ScintiTilePhysVolMap.end())
136  {
137  return 1;
138  }
139  }
140  return 0;
141 }
142 
143 G4VSolid *
145 {
146  double mid_radius = m_InnerRadius + (m_OuterRadius - m_InnerRadius) / 2.;
147  Point_2 p_in_1(mid_radius, 0); // center of scintillator
148 
149  // length of upper edge (middle till outer circle intersect
150  // x/y coordinate of end of center vertical
151  double xcoord = m_ScintiTileThickness / 2. * sin(fabs(m_TiltAngle) / rad) + mid_radius;
152  double ycoord = m_ScintiTileThickness / 2. * cos(fabs(m_TiltAngle) / rad) + 0;
153  Point_2 p_upperedge(xcoord, ycoord);
154  Line_2 s2(p_in_1, p_upperedge); // center vertical
155 
156  Line_2 perp = s2.perpendicular(p_upperedge);
157  Point_2 sc1(m_OuterRadius, 0), sc2(0, m_OuterRadius), sc3(-m_OuterRadius, 0);
158  Circle_2 outer_circle(sc1, sc2, sc3);
159  vector<CGAL::Object> res;
160  CGAL::intersection(outer_circle, perp, std::back_inserter(res));
161  Point_2 upperright;
162  vector<CGAL::Object>::const_iterator iter;
163  for (iter = res.begin(); iter != res.end(); ++iter)
164  {
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))
167  {
168  if (CGAL::to_double(point->first.x()) > CGAL::to_double(p_upperedge.x()))
169  {
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());
172  // the scintillator is twice as long
173  m_ScintiTileXUpper = sqrt(deltax * deltax + deltay * deltay); //
174  Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
175  upperright = pntmp;
176  }
177  }
178  else
179  {
180  cout << "CGAL::Object type not pair..." << endl;
181  }
182  }
183  // length of lower edge (middle till inner circle intersect
184  xcoord = mid_radius - m_ScintiTileThickness / 2. * sin(fabs(m_TiltAngle) / rad);
185  ycoord = 0 - m_ScintiTileThickness / 2. * cos(fabs(m_TiltAngle) / rad);
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);
189  Point_2 ic1(m_InnerRadius, 0), ic2(0, m_InnerRadius), ic3(-m_InnerRadius, 0);
190  Circle_2 inner_circle(ic1, ic2, ic3);
191  res.clear();
192  CGAL::intersection(inner_circle, l_lower, std::back_inserter(res));
193  Point_2 lowerleft;
194  // we have 2 intersections - we want the one furthest to the right (largest x). The correct one is
195  // certainly > 0 but the second one depends on the tilt angle and might also be > 0
196  double minx = 0;
197  for (iter = res.begin(); iter != res.end(); ++iter)
198  {
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))
201  {
202  if (CGAL::to_double(point->first.x()) > minx)
203  {
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());
207  m_ScintiTileXLower = sqrt(deltax * deltax + deltay * deltay);
208  Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
209  lowerleft = pntmp;
210  }
211  }
212  }
215  G4VSolid *scintibox = new G4Box("ScintiTile", m_ScintiTileX / 2., m_ScintiTileThickness / 2., m_ScintiTileZ / 2.);
217  return scintibox;
218 }
219 
220 G4VSolid *
222 {
223  // calculate steel plate on top of the scinti box. Lower edge is the upper edge of
224  // the scintibox + 1/2 the airgap
225  double mid_radius = m_InnerRadius + (m_OuterRadius - m_InnerRadius) / 2.;
226  // first the lower edge, just like the scinti box, just add the air gap
227  // and calculate intersection of edge with inner and outer radius.
228  Point_2 p_in_1(mid_radius, 0); // center of lower scintillator
229  double angle_mid_scinti = M_PI / 2. + m_TiltAngle / rad;
230  double xcoord = m_ScintiInnerGap / 2. * cos(angle_mid_scinti / rad) + mid_radius;
231  double ycoord = m_ScintiInnerGap / 2. * sin(angle_mid_scinti / rad) + 0;
232  Point_2 p_loweredge(xcoord, ycoord);
233  Line_2 s2(p_in_1, p_loweredge); // center vertical
234  Line_2 perp = s2.perpendicular(p_loweredge); // that is the lower edge of the steel plate
235  Point_2 sc1(m_InnerRadius, 0), sc2(0, m_InnerRadius), sc3(-m_InnerRadius, 0);
236  Circle_2 inner_circle(sc1, sc2, sc3);
237  vector<CGAL::Object> res;
238  CGAL::intersection(inner_circle, perp, std::back_inserter(res));
239  Point_2 lowerleft;
240  vector<CGAL::Object>::const_iterator iter;
241  for (iter = res.begin(); iter != res.end(); ++iter)
242  {
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))
245  {
246  if (CGAL::to_double(point->first.x()) > 0)
247  {
248  Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
249  lowerleft = pntmp;
250  }
251  }
252  else
253  {
254  cout << "CGAL::Object type not pair..." << endl;
255  }
256  }
257 
258  double xcoord2 = m_ScintiOuterGap / 2. * cos(angle_mid_scinti / rad) + mid_radius;
259  double ycoord2 = m_ScintiOuterGap / 2. * sin(angle_mid_scinti / rad) + 0;
260  Point_2 p_loweredge2(xcoord2, ycoord2);
261  Line_2 s2_2(p_in_1, p_loweredge2); // center vertical
262  Line_2 perp2 = s2_2.perpendicular(p_loweredge2); // that is the lower edge of the steel plate
263 
264  Point_2 so1(m_OuterRadius, 0), so2(0, m_OuterRadius), so3(-m_OuterRadius, 0);
265  Circle_2 outer_circle(so1, so2, so3);
266  res.clear(); // just clear the content from the last intersection search
267  CGAL::intersection(outer_circle, perp2, std::back_inserter(res));
268  Point_2 lowerright;
269  for (iter = res.begin(); iter != res.end(); ++iter)
270  {
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))
273  {
274  if (CGAL::to_double(point->first.x()) > CGAL::to_double(p_loweredge.x()))
275  {
276  Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
277  lowerright = pntmp;
278  }
279  }
280  else
281  {
282  cout << "CGAL::Object type not pair..." << endl;
283  }
284  }
285  // now we have the lower left and rigth corner, now find the upper edge
286  // find the center of the upper scintilator
287 
288  double phi_midpoint = 2 * M_PI / m_NumScintiPlates;
289  double xmidpoint = cos(phi_midpoint) * mid_radius;
290  double ymidpoint = sin(phi_midpoint) * mid_radius;
291  // angle of perp line at center of scintillator
292  angle_mid_scinti = (M_PI / 2. - phi_midpoint) - (M_PI / 2. + m_TiltAngle / rad);
293  double xcoordup = xmidpoint - m_ScintiInnerGap / 2. * sin(angle_mid_scinti / rad);
294  double ycoordup = ymidpoint - m_ScintiInnerGap / 2. * cos(angle_mid_scinti / rad);
295  Point_2 upperleft;
296  Point_2 upperright;
297  Point_2 mid_upperscint(xmidpoint, ymidpoint);
298  Point_2 p_upperedge(xcoordup, ycoordup);
299  {
300  Line_2 sup(mid_upperscint, p_upperedge); // center vertical
301  Line_2 perp = sup.perpendicular(p_upperedge); // that is the upper edge of the steel plate
302  Point_2 sc1(m_InnerRadius, 0), sc2(0, m_InnerRadius), sc3(-m_InnerRadius, 0);
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;
307  double pxmax = 0.;
308  for (iter = res.begin(); iter != res.end(); ++iter)
309  {
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))
312  {
313  if (CGAL::to_double(point->first.x()) > pxmax)
314  {
315  pxmax = CGAL::to_double(point->first.x());
316  Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
317  upperleft = pntmp;
318  }
319  }
320  else
321  {
322  cout << "CGAL::Object type not pair..." << endl;
323  }
324  }
325 
326  double xcoordup2 = xmidpoint - m_ScintiOuterGap / 2. * sin(angle_mid_scinti / rad);
327  double ycoordup2 = ymidpoint - m_ScintiOuterGap / 2. * cos(angle_mid_scinti / rad);
328  Point_2 p_upperedge2(xcoordup2, ycoordup2);
329  Line_2 sup2(mid_upperscint, p_upperedge2); // center vertical
330  Line_2 perp2 = sup2.perpendicular(p_upperedge2); // that is the upper edge of the steel plate
331 
332  Point_2 so1(m_OuterRadius, 0), so2(0, m_OuterRadius), so3(-m_OuterRadius, 0);
333  Circle_2 outer_circle(so1, so2, so3);
334  res.clear(); // just clear the content from the last intersection search
335  CGAL::intersection(outer_circle, perp2, std::back_inserter(res));
336  for (iter = res.begin(); iter != res.end(); ++iter)
337  {
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))
340  {
341  if (CGAL::to_double(point->first.x()) > CGAL::to_double(p_loweredge.x()))
342  {
343  Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
344  upperright = pntmp;
345  }
346  }
347  else
348  {
349  cout << "CGAL::Object type not pair..." << endl;
350  }
351  }
352  }
353  // the left corners are on a secant with the inner boundary, they need to be shifted
354  // to be a tangent at the center
355  ShiftSecantToTangent(lowerleft, upperleft, upperright, lowerright);
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);
365  G4TwoVector zero(0, 0);
366  G4VSolid *steel_plate = new G4ExtrudedSolid("SteelPlate",
367  vertexes,
368  m_SizeZ / 2.0,
369  zero, 1.0,
370  zero, 1.0);
371 
372  // DisplayVolume(steel_plate, hcalenvelope);
374  return steel_plate;
375 }
376 
377 void PHG4InnerHcalDetector::ShiftSecantToTangent(Point_2 &lowleft, Point_2 &upleft, Point_2 &upright, Point_2 &lowright)
378 {
379  Line_2 secant(lowleft, upleft);
380  Segment_2 upedge(upleft, upright);
381  Segment_2 lowedge(lowleft, lowright);
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.;
384  Point_2 midpoint(xmid, ymid);
385  Line_2 sekperp = secant.perpendicular(midpoint);
386  Point_2 sc1(m_InnerRadius, 0), sc2(0, m_InnerRadius), sc3(-m_InnerRadius, 0);
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;
391  double pxmax = 0.;
392  Point_2 tangtouch;
393  for (iter = res.begin(); iter != res.end(); ++iter)
394  {
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))
397  {
398  if (CGAL::to_double(point->first.x()) > pxmax)
399  {
400  pxmax = CGAL::to_double(point->first.x());
401  Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
402  tangtouch = pntmp;
403  }
404  }
405  else
406  {
407  cout << "CGAL::Object type not pair..." << endl;
408  }
409  }
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))
413  {
414  upleft = *ipoint;
415  }
416  result = CGAL::intersection(lowedge, leftside);
417  if (const Point_2 *ipoint = CGAL::object_cast<Point_2>(&result))
418  {
419  lowleft = *ipoint;
420  }
421  return;
422 }
423 
424 // Construct the envelope and the call the
425 // actual inner hcal construction
427 {
429  G4Material *Air = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
430  G4VSolid *hcal_envelope_cylinder = new G4Tubs("InnerHcal_envelope_solid", m_EnvelopeInnerRadius, m_EnvelopeOuterRadius, m_EnvelopeZ / 2., 0, 2 * M_PI);
431  m_VolumeEnvelope = hcal_envelope_cylinder->GetCubicVolume();
432  G4LogicalVolume *hcal_envelope_log = new G4LogicalVolume(hcal_envelope_cylinder, Air, G4String("Hcal_envelope"), 0, 0, 0);
433 
434  G4RotationMatrix hcal_rotm;
435  hcal_rotm.rotateX(m_Params->get_double_param("rot_x") * deg);
436  hcal_rotm.rotateY(m_Params->get_double_param("rot_y") * deg);
437  hcal_rotm.rotateZ(m_Params->get_double_param("rot_z") * deg);
438  G4VPhysicalVolume *mothervol = new G4PVPlacement(G4Transform3D(hcal_rotm, G4ThreeVector(m_Params->get_double_param("place_x") * cm, m_Params->get_double_param("place_y") * cm, m_Params->get_double_param("place_z") * cm)), hcal_envelope_log, "InnerHcalEnvelope", logicWorld, 0, false, OverlapCheck());
439  m_DisplayAction->SetMyTopVolume(mothervol);
440  ConstructInnerHcal(hcal_envelope_log);
441  vector<G4VPhysicalVolume *>::iterator it = m_ScintiMotherAssembly->GetVolumesIterator();
442  for (unsigned int i = 0; i < m_ScintiMotherAssembly->TotalImprintedVolumes(); i++)
443  {
444  // G4AssemblyVolumes naming convention:
445  // av_WWW_impr_XXX_YYY_ZZZ
446  // where:
447 
448  // WWW - assembly volume instance number
449  // XXX - assembly volume imprint number
450  // YYY - the name of the placed logical volume
451  // ZZZ - the logical volume index inside the assembly volume
452  // e.g. av_1_impr_82_HcalInnerScinti_11_pv_11
453  // 82 the number of the scintillator mother volume
454  // HcalInnerScinti_11: name of scintillator slat
455  // 11: number of scintillator slat logical volume
456  // use boost tokenizer to separate the _, then take value
457  // after "impr" for mother volume and after "pv" for scintillator slat
458  // use boost lexical cast for string -> int conversion
459  // the CopyNo is the mother volume + scinti id
460  // so we can use the CopyNo rather than decoding the string further
461  // looking for "pv"
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)
466  {
467  if (*tokeniter == "impr")
468  {
469  ++tokeniter;
470  if (tokeniter != tok.end())
471  {
472  int layer_id = boost::lexical_cast<int>(*tokeniter);
473  // check detector description, for assemblyvolumes it is not possible
474  // to give the first volume id=0, so they go from id=1 to id=n.
475  // I am not going to start with fortran again - our indices start
476  // at zero, id=0 to id=n-1. So subtract one here
477  int tower_id = (*it)->GetCopyNo() - layer_id;
478  layer_id--;
479  pair<int, int> layer_twr = make_pair(layer_id, tower_id);
480  m_ScintiTilePhysVolMap.insert(pair<G4VPhysicalVolume *, pair<int, int>>(*it, layer_twr));
481  if (layer_id < 0 || layer_id >= m_NumScintiPlates)
482  {
483  cout << "invalid scintillator row " << layer_id
484  << ", valid range 0 < row < " << m_NumScintiPlates << endl;
485  gSystem->Exit(1);
486  }
487  }
488  else
489  {
490  cout << PHWHERE << " Error parsing " << (*it)->GetName()
491  << " for mother volume number " << endl;
492  gSystem->Exit(1);
493  }
494  break;
495  }
496  }
497  ++it;
498  }
499  return;
500 }
501 
503 {
505  SetTiltViaNcross(); // if number of crossings is set, use it to determine tilt
506  CheckTiltAngle(); // die if the tilt angle is out of range
507  G4VSolid *steel_plate = ConstructSteelPlate(hcalenvelope);
508  G4LogicalVolume *steel_logical = new G4LogicalVolume(steel_plate, GetDetectorMaterial(m_Params->get_string_param("material")), "HcalInnerSteelPlate", 0, 0, 0);
509  m_DisplayAction->AddSteelVolume(steel_logical);
511  double phi = 0;
512  double deltaphi = 2 * M_PI / m_NumScintiPlates;
513  ostringstream name;
514  double middlerad = m_OuterRadius - (m_OuterRadius - m_InnerRadius) / 2.;
515  // for shorter scintillators we have to shift by 1/2 the "missing length"
516  double scinti_tile_orig_length = m_ScintiTileXUpper + m_ScintiTileXLower - subtract_from_scinti_x;
517  double shiftslat = fabs(m_ScintiTileXLower - m_ScintiTileXUpper) / 2. + (scinti_tile_orig_length - m_ScintiTileX) / 2.;
518  // calculate phi offset (copied from code inside following loop):
519  // first get the center point (phi=0) so it's middlerad/0
520  // then shift the scintillator center as documented in loop
521  // then
522  // for positive tilt angles we need the lower left corner of the scintillator
523  // for negative tilt angles we need the upper right corner of the scintillator
524  // as it turns out the code uses the middle of the face of the scintillator
525  // as reference, if this is a problem the code needs to be modified to
526  // actually calculate the corner (but the math of the construction is that
527  // the middle of the scintillator sits at zero)
528  double xp = cos(phi) * middlerad;
529  double yp = sin(phi) * middlerad;
530  xp -= cos((-m_TiltAngle) / rad - phi) * shiftslat;
531  yp += sin((-m_TiltAngle) / rad - phi) * shiftslat;
532  if (m_TiltAngle > 0)
533  {
534  double xo = xp - (m_ScintiTileX / 2.) * cos(m_TiltAngle / rad);
535  double yo = yp - (m_ScintiTileX / 2.) * sin(m_TiltAngle / rad);
536  phi = -atan(yo / xo);
537  }
538  else if (m_TiltAngle < 0)
539  {
540  double xo = xp + (m_ScintiTileX / 2.) * cos(m_TiltAngle / rad);
541  double yo = yp + (m_ScintiTileX / 2.) * sin(m_TiltAngle / rad);
542  phi = -atan(yo / xo);
543  }
544  // else (for m_TiltAngle = 0) phi stays zero
545  for (int i = 0; i < m_NumScintiPlates; i++)
546  {
547  G4RotationMatrix *Rot = new G4RotationMatrix();
548  double ypos = sin(phi) * middlerad;
549  double xpos = cos(phi) * middlerad;
550  // the center of the scintillator is not the center of the inner hcal
551  // but depends on the tilt angle. Therefore we need to shift
552  // the center from the mid point
553  ypos += sin((-m_TiltAngle) / rad - phi) * shiftslat;
554  xpos -= cos((-m_TiltAngle) / rad - phi) * shiftslat;
555  Rot->rotateZ(phi * rad + m_TiltAngle);
556  G4ThreeVector g4vec(xpos, ypos, 0);
557  // great the MakeImprint always adds 1 to the copy number and 0 has a
558  // special meaning (which then also adds 1). Basically our volume names
559  // will start at 1 instead of 0 and there is nothing short of patching this
560  // method. I'll take care of this in the decoding of the volume name
561  // AAAAAAARHGS
562  m_ScintiMotherAssembly->MakeImprint(hcalenvelope, g4vec, Rot, i, OverlapCheck());
563  delete Rot;
564  Rot = new G4RotationMatrix();
565  Rot->rotateZ(-phi * rad);
566  name.str("");
567  name << "InnerHcalSteel_" << i;
568  m_SteelAbsorberPhysVolSet.insert(new G4PVPlacement(Rot, G4ThreeVector(0, 0, 0), steel_logical, name.str().c_str(), hcalenvelope, 0, i, OverlapCheck()));
569  phi += deltaphi;
570  }
571  return 0;
572 }
573 
574 // split the big scintillator into tiles covering eta ranges
575 // since they are tilted it is not a straightforward theta cut
576 // it starts at a given eta at the inner radius but the outer radius needs adjusting
578 {
579  G4VSolid *bigtile = ConstructScintillatorBox(hcalenvelope);
580 
581  // eta->theta
582  // scinti_eta_coverage takes precedence
583  G4double eta_cov = m_Params->get_double_param("scinti_eta_coverage");
584  if(eta_cov>0){
585  m_Params->set_double_param("scinti_eta_coverage_pos", eta_cov);
586  m_Params->set_double_param("scinti_eta_coverage_neg", eta_cov);
587  }
588  G4double delta_eta_pos = m_Params->get_double_param("scinti_eta_coverage_pos") / m_NumScintiTilesPos;
589  G4double delta_eta_neg = m_Params->get_double_param("scinti_eta_coverage_neg") / m_NumScintiTilesNeg;
590 
591  G4double eta = 0;
592  G4double theta;
593  G4double x[4];
594  G4double z[4];
595  ostringstream name;
596  double overhang = (m_ScintiTileX - (m_OuterRadius - m_InnerRadius)) / 2.;
597  double offset = 1 * cm + overhang; // add 1cm to make sure the G4ExtrudedSolid
598  // is larger than the tile so we do not have
599  // funny edge effects when overlapping vols
600  double scinti_gap_neighbor = m_Params->get_double_param("scinti_gap_neighbor") * cm;
601 
602  // Positive first, then negative
603 
604  for (int i = 0; i < m_NumScintiTilesPos; i++)
605  {
606  theta = M_PI / 2 - PHG4Utils::get_theta(eta); // theta = 90 for eta=0
607  x[0] = m_InnerRadius - overhang;
608  z[0] = tan(theta) * m_InnerRadius;
609  x[1] = m_OuterRadius + overhang; // since the tile is tilted, x is not at the outer radius but beyond
610  z[1] = tan(theta) * m_OuterRadius;
611  eta += delta_eta_pos;
612  theta = M_PI / 2 - PHG4Utils::get_theta(eta); // theta = 90 for eta=0
613  x[2] = m_InnerRadius - overhang;
614  z[2] = tan(theta) * m_InnerRadius;
615  x[3] = m_OuterRadius + overhang; // since the tile is tilted, x is not at the outer radius but beyond
616  z[3] = tan(theta) * m_OuterRadius;
617  // apply gap between scintillators
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]);
624  x[0] = m_InnerRadius - offset;
625  z[0] = x_at_y(leftsidelow, leftsidehigh, x[0]);
626  x[1] = m_OuterRadius + offset;
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]);
630  x[2] = m_OuterRadius + offset;
631  z[2] = x_at_y(rightsidelow, rightsidehigh, x[2]);
632  x[3] = m_InnerRadius - offset;
633  z[3] = x_at_y(rightsidelow, rightsidehigh, x[3]);
634 
635  vector<G4TwoVector> vertexes;
636  for (int j = 0; j < 4; j++)
637  {
638  G4TwoVector v(x[j], z[j]);
639  vertexes.push_back(v);
640  }
641  G4TwoVector zero(0, 0);
642 
643  G4VSolid *scinti = new G4ExtrudedSolid("ScintillatorTile",
644  vertexes,
645  m_ScintiTileThickness + 0.2 * mm,
646  zero, 1.0,
647  zero, 1.0);
648  G4RotationMatrix *rotm = new G4RotationMatrix();
649  rotm->rotateX(-90 * deg);
650  name.str("");
651  name << "scintillator_" << i << "_left";
652  G4VSolid *scinti_tile = new G4IntersectionSolid(name.str(), bigtile, scinti, rotm, G4ThreeVector(-(m_InnerRadius + m_OuterRadius) / 2., 0, -m_Params->get_double_param("place_z") * cm));
653  delete rotm;
654  m_ScintiTilesVec[i + m_NumScintiTilesNeg] = scinti_tile;
655 
656  }
657 
658  eta = 0.0; // reset
659  for (int i = 0; i < m_NumScintiTilesNeg; i++)
660  {
661  theta = M_PI / 2 - PHG4Utils::get_theta(eta); // theta = 90 for eta=0
662  x[0] = m_InnerRadius - overhang;
663  z[0] = tan(theta) * m_InnerRadius;
664  x[1] = m_OuterRadius + overhang; // since the tile is tilted, x is not at the outer radius but beyond
665  z[1] = tan(theta) * m_OuterRadius;
666  eta += delta_eta_neg;
667  theta = M_PI / 2 - PHG4Utils::get_theta(eta); // theta = 90 for eta=0
668  x[2] = m_InnerRadius - overhang;
669  z[2] = tan(theta) * m_InnerRadius;
670  x[3] = m_OuterRadius + overhang; // since the tile is tilted, x is not at the outer radius but beyond
671  z[3] = tan(theta) * m_OuterRadius;
672  // apply gap between scintillators
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]);
679  x[0] = m_InnerRadius - offset;
680  z[0] = x_at_y(leftsidelow, leftsidehigh, x[0]);
681  x[1] = m_OuterRadius + offset;
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]);
685  x[2] = m_OuterRadius + offset;
686  z[2] = x_at_y(rightsidelow, rightsidehigh, x[2]);
687  x[3] = m_InnerRadius - offset;
688  z[3] = x_at_y(rightsidelow, rightsidehigh, x[3]);
689 
690  vector<G4TwoVector> vertexes;
691  for (int j = 0; j < 4; j++)
692  {
693  G4TwoVector v(x[j], z[j]);
694  vertexes.push_back(v);
695  }
696  G4TwoVector zero(0, 0);
697 
698  G4VSolid *scinti = new G4ExtrudedSolid("ScintillatorTile",
699  vertexes,
700  m_ScintiTileThickness + 0.2 * mm,
701  zero, 1.0,
702  zero, 1.0);
703 
704  G4RotationMatrix *rotm = new G4RotationMatrix();
705  rotm->rotateX(90 * deg);
706  name.str("");
707  name << "scintillator_" << i << "_right";
708  G4VSolid *scinti_tile = new G4IntersectionSolid(name.str(), bigtile, scinti, rotm, G4ThreeVector(-(m_InnerRadius + m_OuterRadius) / 2., 0, -m_Params->get_double_param("place_z") * cm));
709  m_ScintiTilesVec[m_NumScintiTilesNeg - i - 1] = scinti_tile;
710  delete rotm;
711 
712  }
713 
714  // for (unsigned int i=0; i<m_ScintiTilesVec.size(); i++)
715  // {
716  // if (m_ScintiTilesVec[i])
717  // {
718  // DisplayVolume(m_ScintiTilesVec[i],hcalenvelope );
719  // }
720  // }
721 
722  return;
723 }
724 
725 double
727 {
728  double xret = NAN;
729  double x[2];
730  double y[2];
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());
735  Line_2 l(p0, p1);
736  double newx = fabs(x[0]) + fabs(x[1]);
737  Point_2 p0new(-newx, yin);
738  Point_2 p1new(newx, yin);
739  Segment_2 s(p0new, p1new);
740  CGAL::Object result = CGAL::intersection(l, s);
741  if (const Point_2 *ipoint = CGAL::object_cast<Point_2>(&result))
742  {
743  xret = CGAL::to_double(ipoint->x());
744  }
745  else
746  {
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;
750  exit(1);
751  }
752  return xret;
753 }
754 
757 {
758  ConstructHcalSingleScintillators(hcalenvelope);
759  G4AssemblyVolume *assmeblyvol = new G4AssemblyVolume();
760  ostringstream name;
761  G4ThreeVector g4vec;
762 
763  double steplimits = m_Params->get_double_param("steplimits") * cm;
764  for (unsigned int i = 0; i < m_ScintiTilesVec.size(); i++)
765  {
766  name.str("");
767  name << m_ScintiLogicNamePrefix << i;
768  G4UserLimits *g4userlimits = nullptr;
769  if (isfinite(steplimits))
770  {
771  g4userlimits = new G4UserLimits(steplimits);
772  }
773  G4LogicalVolume *scinti_tile_logic = new G4LogicalVolume(m_ScintiTilesVec[i], GetDetectorMaterial("G4_POLYSTYRENE"), name.str(), nullptr, nullptr, g4userlimits);
774  m_DisplayAction->AddScintiVolume(scinti_tile_logic);
775  assmeblyvol->AddPlacedVolume(scinti_tile_logic, g4vec, nullptr);
776  }
777  return assmeblyvol;
778 }
779 
780 // check if tilt angle is reasonable - too large, no intersections with inner radius
782 {
783  if (fabs(m_TiltAngle / rad) >= M_PI)
784  {
785  cout << PHWHERE << "invalid tilt angle, abs(tilt) >= 90 deg: " << (m_TiltAngle / deg)
786  << endl;
787  exit(1);
788  }
789 
790  double mid_radius = m_InnerRadius + (m_OuterRadius - m_InnerRadius) / 2.;
791  Point_2 pmid(mid_radius, 0); // center of scintillator
792  double xcoord = 0;
793  double ycoord = mid_radius * tan(m_TiltAngle / rad);
794  Point_2 pxnull(xcoord, ycoord);
795  Line_2 s2(pmid, pxnull);
796  Point_2 sc1(m_InnerRadius, 0), sc2(0, m_InnerRadius), sc3(-m_InnerRadius, 0);
797  Circle_2 inner_circle(sc1, sc2, sc3);
798  vector<CGAL::Object> res;
799  CGAL::intersection(inner_circle, s2, std::back_inserter(res));
800  if (res.size() == 0)
801  {
802  cout << PHWHERE << " Tilt angle " << (m_TiltAngle / deg)
803  << " too large, no intersection with inner radius" << endl;
804  exit(1);
805  }
806  return 0;
807 }
808 
810 {
811  // just make sure the parameters make a bit of sense
813  {
814  cout << PHWHERE << ": Inner Radius " << m_InnerRadius / cm
815  << " cm larger than Outer Radius " << m_OuterRadius / cm
816  << " cm" << endl;
817  gSystem->Exit(1);
818  }
820  {
821  cout << PHWHERE << "Scintillator thickness " << m_ScintiTileThickness / cm
822  << " cm larger than scintillator inner gap " << m_ScintiInnerGap / cm
823  << " cm" << endl;
824  gSystem->Exit(1);
825  }
827  {
828  cout << PHWHERE << "Scintillator outer radius " << m_ScintiOuterRadius / cm
829  << " cm smaller than inner radius " << m_InnerRadius / cm
830  << " cm" << endl;
831  gSystem->Exit(1);
832  }
833  return 0;
834 }
835 
837 {
838  // if tilt angle is set it takes precedence
839  // this happens during readback where the tilt angle was set
840  // in this method
841  int ncross = m_Params->get_int_param("ncross");
842  if (!ncross || isfinite(m_TiltAngle))
843  {
844  // flag ncrossing as not used
845  m_Params->set_int_param("ncross", 0);
846  return;
847  }
848  if ((isfinite(m_TiltAngle)) && (Verbosity() > 0))
849  {
850  cout << "both number of crossings and tilt angle are set" << endl;
851  cout << "using number of crossings to determine tilt angle" << endl;
852  }
853  double mid_radius = m_InnerRadius + (m_OuterRadius - m_InnerRadius) / 2.;
854  double deltaphi = (2 * M_PI / m_NumScintiPlates) * ncross;
855  Point_2 pnull(0, 0);
856  Point_2 plow(m_InnerRadius, 0);
857  Point_2 phightmp(1, tan(deltaphi));
858  Point_2 pin1(m_InnerRadius, 0), pin2(0, m_InnerRadius), pin3(-m_InnerRadius, 0);
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);
862  Point_2 pout1(m_OuterRadius, 0), pout2(0, m_OuterRadius), pout3(-m_OuterRadius, 0);
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));
867  Point_2 upperright;
868  vector<CGAL::Object>::const_iterator iter;
869  for (iter = res.begin(); iter != res.end(); ++iter)
870  {
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))
873  {
874  if (CGAL::to_double(point->first.x()) > 0)
875  {
876  Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
877  upperright = pntmp;
878  }
879  }
880  else
881  {
882  cout << "CGAL::Object type not pair..." << endl;
883  exit(1);
884  }
885  }
886  Line_2 l_right(plow, upperright);
887  res.clear();
888  Point_2 midpoint;
889  CGAL::intersection(mid_circle, l_right, std::back_inserter(res));
890  for (iter = res.begin(); iter != res.end(); ++iter)
891  {
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))
894  {
895  if (CGAL::to_double(point->first.x()) > 0)
896  {
897  Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
898  midpoint = pntmp;
899  }
900  }
901  else
902  {
903  cout << "CGAL::Object type not pair..." << endl;
904  exit(1);
905  }
906  }
907  // length left side
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()));
910  // c^2 = a^2+b^2 - 2ab*cos(gamma)
911  // gamma = acos((a^2+b^2=c^2)/2ab
912  double tiltangle = acos((ll * ll + upside * upside - m_InnerRadius * m_InnerRadius) / (2 * ll * upside));
913  tiltangle = tiltangle * rad;
914  m_TiltAngle = copysign(tiltangle, ncross);
915  m_Params->set_double_param("tilt_angle", m_TiltAngle / deg);
916  return;
917 }
918 
919 void PHG4InnerHcalDetector::Print(const string &what) const
920 {
921  cout << "Inner Hcal Detector:" << endl;
922  if (what == "ALL" || what == "VOLUME")
923  {
924  cout << "Volume Envelope: " << m_VolumeEnvelope / cm3 << " cm^3" << endl;
925  cout << "Volume Steel: " << m_VolumeSteel / cm3 << " cm^3" << endl;
926  cout << "Volume Scintillator: " << m_VolumeScintillator / cm3 << " cm^3" << endl;
927  cout << "Volume Air: " << (m_VolumeEnvelope - m_VolumeSteel - m_VolumeScintillator) / cm3 << " cm^3" << endl;
928  }
929  return;
930 }
931 
933 {
934  auto it = m_ScintiTilePhysVolMap.find(volume);
935  if (it != m_ScintiTilePhysVolMap.end())
936  {
937  return it->second;
938  }
939  cout << "could not locate volume " << volume->GetName()
940  << " in Inner Hcal scintillator map" << endl;
941  gSystem->Exit(1);
942  // that's dumb but code checkers do not know that gSystem->Exit()
943  // terminates, so using the standard exit() makes them happy
944  exit(1);
945 }