Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4OuterHcalDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4OuterHcalDetector.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 #include <phool/PHCompositeNode.h>
17 #include <phool/PHIODataNode.h>
18 #include <phool/PHNode.h> // for PHNode
19 #include <phool/PHNodeIterator.h>
20 #include <phool/PHObject.h> // for PHObject
21 #include <phool/getClass.h>
22 
23 #include <calobase/RawTowerDefs.h> // for convert_name_...
24 #include <calobase/RawTowerGeom.h> // for RawTowerGeom
25 #include <calobase/RawTowerGeomContainer.h> // for RawTowerGeomC...
26 #include <calobase/RawTowerGeomContainer_Cylinderv1.h>
27 #include <calobase/RawTowerGeomv1.h>
28 
29 #include <TSystem.h>
30 
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>
49 
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
59 
60 #include <boost/lexical_cast.hpp>
61 #include <boost/tokenizer.hpp>
62 
63 #include <algorithm>
64 #include <cmath>
65 #include <cstdlib>
66 #include <iostream>
67 #include <iterator>
68 #include <sstream>
69 
70 class PHCompositeNode;
71 
72 using Circle_2 = CGAL::Circle_2<PHG4OuterHcalDetector::Circular_k>;
73 using Circular_arc_point_2 = CGAL::Circular_arc_point_2<PHG4OuterHcalDetector::Circular_k>;
74 using Line_2 = CGAL::Line_2<PHG4OuterHcalDetector::Circular_k>;
75 using Segment_2 = CGAL::Segment_2<PHG4OuterHcalDetector::Circular_k>;
76 
77 // just for debugging if you want a single layer of scintillators at the center of the world
78 //#define SCINTITEST
79 
80 // face touches the boundary instead of the corner, subtracting 1 permille from the total
81 // scintilator length takes care of this
82 static double subtract_from_scinti_x = 0.1 * mm;
83 
85  : PHG4Detector(subsys, Node, dnam)
86  , m_DisplayAction(dynamic_cast<PHG4OuterHcalDisplayAction *>(subsys->GetDisplayAction()))
87  , m_Params(parames)
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")
105 {
106  m_ScintiTilesVec.assign(2 * m_NumScintiTiles, static_cast<G4VSolid *>(nullptr));
107 }
108 
110 {
111  delete m_ScintiMotherAssembly;
112  delete m_FieldSetup;
113 }
114 
115 //_______________________________________________________________
116 //_______________________________________________________________
118 {
120  {
121  if (m_SteelAbsorberVec.find(volume) != m_SteelAbsorberVec.end())
122  {
123  return -1;
124  }
125  }
126  if (m_ActiveFlag)
127  {
128  if (m_ScintiTilePhysVolMap.find(volume) != m_ScintiTilePhysVolMap.end())
129  {
130  return 1;
131  }
132  }
133  return 0;
134 }
135 
136 G4VSolid *
137 PHG4OuterHcalDetector::ConstructScintillatorBox(G4LogicalVolume * /*hcalenvelope*/)
138 {
139  double mid_radius = m_InnerRadius + (m_OuterRadius - m_InnerRadius) / 2.;
140  PHG4OuterHcalDetector::Point_2 p_in_1(mid_radius, 0); // center of scintillator
141 
142  // length of upper edge (middle till outer circle intersect
143  // x/y coordinate of end of center vertical
144  double xcoord = m_ScintiTileThickness / 2. * sin(fabs(m_TiltAngle) / rad) + mid_radius;
145  double ycoord = m_ScintiTileThickness / 2. * cos(fabs(m_TiltAngle) / rad) + 0;
146  PHG4OuterHcalDetector::Point_2 p_upperedge(xcoord, ycoord);
147  Line_2 s2(p_in_1, p_upperedge); // center vertical
148 
149  Line_2 perp = s2.perpendicular(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)
157  {
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))
160  {
161  if (CGAL::to_double(point->first.x()) > CGAL::to_double(p_upperedge.x()))
162  {
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());
165  // the scintillator is twice as long
166  m_ScintiTileXUpper = sqrt(deltax * deltax + deltay * deltay); //
167  PHG4OuterHcalDetector::Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
168  upperright = pntmp;
169  }
170  }
171  else
172  {
173  std::cout << "CGAL::Object type not pair..." << std::endl;
174  }
175  }
176  // length of lower edge (middle till inner circle intersect
177  xcoord = mid_radius - m_ScintiTileThickness / 2. * sin(fabs(m_TiltAngle) / rad);
178  ycoord = 0 - m_ScintiTileThickness / 2. * cos(fabs(m_TiltAngle) / rad);
179  PHG4OuterHcalDetector::Point_2 p_loweredge(xcoord, ycoord);
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);
184  res.clear();
185  CGAL::intersection(inner_circle, l_lower, std::back_inserter(res));
187  // we have 2 intersections - we want the one furthest to the right (largest x). The correct one is
188  // certainly > 0 but the second one depends on the tilt angle and might also be > 0
189  double minx = 0;
190  for (iter = res.begin(); iter != res.end(); ++iter)
191  {
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))
194  {
195  if (CGAL::to_double(point->first.x()) > minx)
196  {
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());
200  m_ScintiTileXLower = sqrt(deltax * deltax + deltay * deltay);
201  PHG4OuterHcalDetector::Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
202  lowerleft = pntmp;
203  }
204  }
205  }
208  G4VSolid *scintibox = new G4Box("ScintiTile", m_ScintiTileX / 2., m_ScintiTileThickness / 2., m_ScintiTileZ / 2.);
209  m_VolumeScintillator = scintibox->GetCubicVolume() * m_NumScintiPlates;
210  return scintibox;
211 }
212 
213 G4VSolid *
214 PHG4OuterHcalDetector::ConstructSteelPlate(G4LogicalVolume * /*hcalenvelope*/)
215 {
216  // calculate steel plate on top of the scinti box. Lower edge is the upper edge of
217  // the scintibox + 1/2 the airgap
218  double mid_radius = m_InnerRadius + (m_OuterRadius - m_InnerRadius) / 2.;
219  // first the lower edge, just like the scinti box, just add the air gap
220  // and calculate intersection of edge with inner and outer radius.
221  PHG4OuterHcalDetector::Point_2 p_in_1(mid_radius, 0); // center of lower scintillator
222  double angle_mid_scinti = M_PI / 2. + m_TiltAngle / rad;
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;
225  PHG4OuterHcalDetector::Point_2 p_loweredge(xcoord, ycoord);
226  Line_2 s2(p_in_1, p_loweredge); // center vertical
227  Line_2 perp = s2.perpendicular(p_loweredge); // that is the lower edge of the steel plate
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)
235  {
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))
238  {
239  if (CGAL::to_double(point->first.x()) > 0)
240  {
241  PHG4OuterHcalDetector::Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
242  lowerleft = pntmp;
243  }
244  }
245  else
246  {
247  std::cout << "CGAL::Object type not pair..." << std::endl;
248  }
249  }
251  Circle_2 outer_circle(so1, so2, so3);
252  res.clear(); // just clear the content from the last intersection search
253  CGAL::intersection(outer_circle, perp, std::back_inserter(res));
255  for (iter = res.begin(); iter != res.end(); ++iter)
256  {
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))
259  {
260  if (CGAL::to_double(point->first.x()) > CGAL::to_double(p_loweredge.x()))
261  {
262  PHG4OuterHcalDetector::Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
263  lowerright = pntmp;
264  }
265  }
266  else
267  {
268  std::cout << "CGAL::Object type not pair..." << std::endl;
269  }
270  }
271  // now we have the lower left and rigth corner, now find the upper edge
272  // find the center of the upper scintilator
273 
274  double phi_midpoint = 2 * M_PI / m_NumScintiPlates;
275  double xmidpoint = cos(phi_midpoint) * mid_radius;
276  double ymidpoint = sin(phi_midpoint) * mid_radius;
277  // angle of perp line at center of scintillator
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);
283  PHG4OuterHcalDetector::Point_2 mid_upperscint(xmidpoint, ymidpoint);
284  PHG4OuterHcalDetector::Point_2 p_upperedge(xcoordup, ycoordup);
285  Line_2 sup(mid_upperscint, p_upperedge); // center vertical
286  Line_2 perpA = sup.perpendicular(p_upperedge); // that is the upper edge of the steel plate
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;
292  double pxmax = 0.;
293  for (iterA = resA.begin(); iterA != resA.end(); ++iterA)
294  {
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))
297  {
298  if (CGAL::to_double(point->first.x()) > pxmax)
299  {
300  pxmax = CGAL::to_double(point->first.x());
301  PHG4OuterHcalDetector::Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
302  upperleft = pntmp;
303  }
304  }
305  else
306  {
307  std::cout << "CGAL::Object type not pair..." << std::endl;
308  }
309  }
311  Circle_2 outer_circleA(so1A, so2A, so3A);
312  resA.clear(); // just clear the content from the last intersection search
313  CGAL::intersection(outer_circleA, perpA, std::back_inserter(resA));
314  for (iterA = resA.begin(); iterA != resA.end(); ++iterA)
315  {
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))
318  {
319  if (CGAL::to_double(point->first.x()) > CGAL::to_double(p_loweredge.x()))
320  {
321  PHG4OuterHcalDetector::Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
322  upperright = pntmp;
323  }
324  }
325  else
326  {
327  std::cout << "CGAL::Object type not pair..." << std::endl;
328  }
329  }
330  // the left corners are on a secant with the inner boundary, they need to be shifted
331  // to be a tangent at the center
332  ShiftSecantToTangent(lowerleft, upperleft, upperright, lowerright);
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",
344  vertexes,
345  m_SizeZ / 2.0,
346  zero, 1.0,
347  zero, 1.0);
348 
349  m_VolumeSteel = steel_plate_uncut->GetCubicVolume() * m_NumScintiPlates;
350  // now cut out space for magnet at the ends
352  {
353  return steel_plate_uncut;
354  }
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));
358  delete rotm;
359  // DisplayVolume(steel_plate_uncut, hcalenvelope);
360  // DisplayVolume(m_SteelCutoutForMagnetG4Solid, hcalenvelope);
361  // DisplayVolume(m_SteelCutoutForMagnetG4Solid, hcalenvelope,rotm);
362  // DisplayVolume(steel_firstcut_solid, hcalenvelope);
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));
366  // DisplayVolume(steel_cut_solid, hcalenvelope);
367  delete rotm;
368  return steel_cut_solid;
369 }
370 
372 {
373  Line_2 secant(lowleft, upleft);
374  Segment_2 upedge(upleft, upright);
375  Segment_2 lowedge(lowleft, lowright);
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.;
378  PHG4OuterHcalDetector::Point_2 midpoint(xmid, ymid);
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;
385  double pxmax = 0.;
387  for (iter = res.begin(); iter != res.end(); ++iter)
388  {
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))
391  {
392  if (CGAL::to_double(point->first.x()) > pxmax)
393  {
394  pxmax = CGAL::to_double(point->first.x());
395  PHG4OuterHcalDetector::Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
396  tangtouch = pntmp;
397  }
398  }
399  else
400  {
401  std::cout << "CGAL::Object type not pair..." << std::endl;
402  }
403  }
404  Line_2 leftside = sekperp.perpendicular(tangtouch);
405  CGAL::Object result = CGAL::intersection(upedge, leftside);
406  if (const PHG4OuterHcalDetector::Point_2 *ipoint = CGAL::object_cast<PHG4OuterHcalDetector::Point_2>(&result))
407  {
408  upleft = *ipoint;
409  }
410  result = CGAL::intersection(lowedge, leftside);
411  if (const PHG4OuterHcalDetector::Point_2 *ipoint = CGAL::object_cast<PHG4OuterHcalDetector::Point_2>(&result))
412  {
413  lowleft = *ipoint;
414  }
415  return;
416 }
417 
418 void PHG4OuterHcalDetector::ConstructMe(G4LogicalVolume *logicWorld)
419 {
420 #ifdef SCINTITEST
421  ConstructOuterHcal(logicWorld);
422  return;
423 #endif
425  G4Material *Air = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
426  G4VSolid *hcal_envelope_cylinder = new G4Tubs("OuterHcal_envelope_solid", m_EnvelopeInnerRadius, m_EnvelopeOuterRadius, m_EnvelopeZ / 2., 0, 2 * M_PI);
427  m_VolumeEnvelope = hcal_envelope_cylinder->GetCubicVolume();
428  G4LogicalVolume *hcal_envelope_log = new G4LogicalVolume(hcal_envelope_cylinder, Air, G4String("OuterHcal_envelope"), nullptr, nullptr, nullptr);
429  G4RotationMatrix hcal_rotm;
430  hcal_rotm.rotateX(m_Params->get_double_param("rot_x") * deg);
431  hcal_rotm.rotateY(m_Params->get_double_param("rot_y") * deg);
432  hcal_rotm.rotateZ(m_Params->get_double_param("rot_z") * deg);
433  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, "OuterHcal", logicWorld, false, false, OverlapCheck());
434  m_DisplayAction->SetMyTopVolume(mothervol);
435  ConstructOuterHcal(hcal_envelope_log);
436  std::vector<G4VPhysicalVolume *>::iterator it = m_ScintiMotherAssembly->GetVolumesIterator();
437  for (unsigned int i = 0; i < m_ScintiMotherAssembly->TotalImprintedVolumes(); i++)
438  {
439  // G4AssemblyVolumes naming convention:
440  // av_WWW_impr_XXX_YYY_ZZZ
441  // where:
442 
443  // WWW - assembly volume instance number
444  // XXX - assembly volume imprint number
445  // YYY - the name of the placed logical volume
446  // ZZZ - the logical volume index inside the assembly volume
447  // e.g. av_1_impr_82_HcalInnerScinti_11_pv_11
448  // 82 the number of the scintillator mother volume
449  // HcalInnerScinti_11: name of scintillator slat
450  // 11: number of scintillator slat logical volume
451  // use boost tokenizer to separate the _, then take value
452  // after "impr" for mother volume and after "pv" for scintillator slat
453  // use boost lexical cast for string -> int conversion
454  // the CopyNo is the mother volume + scinti id
455  // so we can use the CopyNo rather than decoding the string further
456  // looking for "pv"
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)
461  {
462  if (*tokeniter == "impr")
463  {
464  ++tokeniter;
465  if (tokeniter != tok.end())
466  {
467  int layer_id = boost::lexical_cast<int>(*tokeniter);
468  // check detector description, for assemblyvolumes it is not possible
469  // to give the first volume id=0, so they go from id=1 to id=n.
470  // I am not going to start with fortran again - our indices start
471  // at zero, id=0 to id=n-1. So subtract one here
472  int tower_id = (*it)->GetCopyNo() - layer_id;
473  layer_id--;
474  std::pair<int, int> layer_twr = std::make_pair(layer_id, tower_id);
475  m_ScintiTilePhysVolMap.insert(std::pair<G4VPhysicalVolume *, std::pair<int, int>>(*it, layer_twr));
476  if (layer_id < 0 || layer_id >= m_NumScintiPlates)
477  {
478  std::cout << "invalid scintillator row " << layer_id
479  << ", valid range 0 < row < " << m_NumScintiPlates << std::endl;
480  gSystem->Exit(1);
481  }
482  }
483  else
484  {
485  std::cout << PHWHERE << " Error parsing " << (*it)->GetName()
486  << " for mother volume number " << std::endl;
487  gSystem->Exit(1);
488  }
489  break;
490  }
491  }
492  ++it;
493  }
494  if(!m_Params->get_int_param("saveg4hit")) AddGeometryNode();
495  return;
496 }
497 
498 int PHG4OuterHcalDetector::ConstructOuterHcal(G4LogicalVolume *hcalenvelope)
499 {
501  SetTiltViaNcross(); // if number of crossings is set, use it to determine tilt
502  CheckTiltAngle(); // die if the tilt angle is out of range
503  // the needed steel cutout volume for the magnet is constructed with
504  // the scintillators since we have the theta angle at that point
505 
506  // call field setup here where we have the calculated tilt angle if number
507  // of crossings is given
509  m_NumScintiPlates, /*G4int steelPlates*/
510  m_ScintiGap, /*G4double scintiGap*/
511  m_TiltAngle); /*G4double tiltAngle*/
512 
514 #ifdef SCINTITEST
515  return 0;
516 #endif
517  G4VSolid *steel_plate = ConstructSteelPlate(hcalenvelope);
518  // DisplayVolume(steel_plate_4 ,hcalenvelope);
519  G4LogicalVolume *steel_logical = new G4LogicalVolume(steel_plate, GetDetectorMaterial(m_Params->get_string_param("material")), "HcalOuterSteelPlate", nullptr, nullptr, nullptr);
520  m_DisplayAction->AddSteelVolume(steel_logical);
521  double phi = 0;
522  double deltaphi = 2 * M_PI / m_NumScintiPlates;
523  std::ostringstream name;
524  double middlerad = m_OuterRadius - (m_OuterRadius - m_InnerRadius) / 2.;
525  // okay this is crude. Since the inner and outer radius of the scintillator is different from the inner/outer
526  // radius of the steel so the scintillator needs some shifting to get the gaps right
527  // basically the first shift (sumshift) puts it at the inner radius of the steel
528  // then it needs to be shifted up by twice the difference between the inner steel and inner scintillator radius
529  // since sumshift is div by 2 for G4 but we need to shift the full delta(inner radius)
530  double scinti_tile_orig_length = m_ScintiTileXUpper + m_ScintiTileXLower - subtract_from_scinti_x;
531  double shiftup = (m_ScintiInnerRadius - m_InnerRadius) / cos(m_TiltAngle / rad);
532  double sumshift = scinti_tile_orig_length - m_ScintiTileX;
533  sumshift = sumshift - 2 * shiftup;
534  double shiftslat = fabs(m_ScintiTileXLower - m_ScintiTileXUpper) / 2. + sumshift / 2.;
535  // calculate phi offset (copied from code inside following loop):
536  // first get the center point (phi=0) so it's middlerad/0
537  // then shift the scintillator center as documented in loop
538  // then
539  // for positive tilt angles we need the lower left corner of the scintillator
540  // for negative tilt angles we nee the upper right corner of the scintillator
541  // as it turns out the code uses the middle of the face of the scintillator
542  // as reference, if this is a problem the code needs to be modified to
543  // actually calculate the corner (but the math of the construction is that
544  // the middle of the scintillator sits at zero)
545  double xp = cos(phi) * middlerad;
546  double yp = sin(phi) * middlerad;
547  xp -= cos((-m_TiltAngle) / rad - phi) * shiftslat;
548  yp += sin((-m_TiltAngle) / rad - phi) * shiftslat;
549  if (m_TiltAngle > 0)
550  {
551  double xo = xp - (m_ScintiTileX / 2.) * cos(m_TiltAngle / rad);
552  double yo = yp - (m_ScintiTileX / 2.) * sin(m_TiltAngle / rad);
553  phi = -atan(yo / xo);
554  }
555  else if (m_TiltAngle < 0)
556  {
557  double xo = xp + (m_ScintiTileX / 2.) * cos(m_TiltAngle / rad);
558  double yo = yp + (m_ScintiTileX / 2.) * sin(m_TiltAngle / rad);
559  phi = -atan(yo / xo);
560  }
561  // else (for m_TiltAngle = 0) phi stays zero
562  for (int i = 0; i < m_NumScintiPlates; i++)
563  {
564  G4RotationMatrix *Rot = new G4RotationMatrix();
565  double ypos = sin(phi) * middlerad;
566  double xpos = cos(phi) * middlerad;
567  // the center of the scintillator is not the center of the inner hcal
568  // but depends on the tilt angle. Therefore we need to shift
569  // the center from the mid point
570  ypos += sin((-m_TiltAngle) / rad - phi) * shiftslat;
571  xpos -= cos((-m_TiltAngle) / rad - phi) * shiftslat;
572  Rot->rotateZ(phi * rad + m_TiltAngle);
573  G4ThreeVector g4vec(xpos, ypos, 0);
574  // great the MakeImprint always adds 1 to the copy number and 0 has a
575  // special meaning (which then also adds 1). Basically our volume names
576  // will start at 1 instead of 0 and there is nothing short of patching this
577  // method. I'll take care of this in the decoding of the volume name
578  // AAAAAAARHGS
579  m_ScintiMotherAssembly->MakeImprint(hcalenvelope, g4vec, Rot, i, OverlapCheck());
580  delete Rot;
581  Rot = new G4RotationMatrix();
582  Rot->rotateZ(-phi * rad);
583  name.str("");
584  name << "OuterHcalSteel_" << i;
585  m_SteelAbsorberVec.insert(new G4PVPlacement(Rot, G4ThreeVector(0, 0, 0), steel_logical, name.str(), hcalenvelope, false, i, OverlapCheck()));
586  phi += deltaphi;
587  }
588  hcalenvelope->SetFieldManager(m_FieldSetup->get_Field_Manager_Gap(), false);
589 
590  steel_logical->SetFieldManager(m_FieldSetup->get_Field_Manager_Iron(), true);
591  return 0;
592 }
593 
595 {
596  G4VSolid *bigtile = ConstructScintillatorBox(hcalenvelope);
597  // eta->theta
598  G4double delta_eta = m_Params->get_double_param("scinti_eta_coverage") / m_NumScintiTiles;
599  G4double eta = 0;
600  G4double theta;
601  G4double x[4];
602  G4double z[4];
603  std::ostringstream name;
604  double overhang = (m_ScintiTileX - (m_ScintiOuterRadius - m_ScintiInnerRadius)) / 2.;
605  double offset = 1 * cm + overhang; // add 1cm to make sure the G4ExtrudedSolid
606  // is larger than the tile so we do not have
607  // funny edge effects when overlapping vols
608  double magnet_cutout_x = (m_Params->get_double_param("magnet_cutout_scinti_radius") * cm - m_ScintiInnerRadius) / cos(m_TiltAngle / rad);
609  double x_inner = m_ScintiInnerRadius - overhang;
610  double inner_offset = offset;
611  // coordinates like the steel plates:
612  // 0/0 upper left
613  // 1/1 upper right
614  // 2/2 lower right
615  // 3/3 lower left
616  // sorry they are different than the coordinates used for the scintilators
617  // here, this is why the indices are seemingly mixed up
618  double xsteelcut[4];
619  double zsteelcut[4];
620  std::fill_n(zsteelcut, 4, NAN);
622  double steel_offset = 1 * cm + steel_overhang; // add 1cm to make sure the G4ExtrudedSolid
623  double steel_x_inner = m_InnerRadius - steel_overhang;
624  double steel_magnet_cutout_x = (m_Params->get_double_param("magnet_cutout_radius") * cm - m_InnerRadius) / cos(m_TiltAngle / rad);
625  double steel_inner_offset = steel_offset;
626  xsteelcut[0] = steel_x_inner + steel_magnet_cutout_x;
627  xsteelcut[1] = xsteelcut[0];
628  xsteelcut[2] = m_InnerRadius - steel_offset;
629  xsteelcut[3] = xsteelcut[2];
630  double scinti_gap_neighbor = m_Params->get_double_param("scinti_gap_neighbor") * cm;
631  for (int i = 0; i < m_NumScintiTiles; i++)
632  {
633  if (i >= m_Params->get_int_param("magnet_cutout_first_scinti"))
634  {
635  x_inner = m_ScintiInnerRadius - overhang + magnet_cutout_x;
636  inner_offset = offset - magnet_cutout_x;
637  }
638  theta = M_PI / 2 - PHG4Utils::get_theta(eta); // theta = 90 for eta=0
639  x[0] = x_inner;
640  z[0] = tan(theta) * m_ScintiInnerRadius;
641  x[1] = m_ScintiOuterRadius + overhang; // since the tile is tilted, x is not at the outer radius but beyond
642  z[1] = tan(theta) * m_ScintiOuterRadius;
643  if (i >= m_Params->get_int_param("magnet_cutout_first_scinti"))
644  {
645  z[0] = tan(theta) * (m_ScintiInnerRadius + (m_Params->get_double_param("magnet_cutout_scinti_radius") * cm - m_ScintiInnerRadius));
646  }
647  eta += delta_eta;
648  theta = M_PI / 2 - PHG4Utils::get_theta(eta); // theta = 90 for eta=0
649  x[2] = x_inner;
650  z[2] = tan(theta) * m_ScintiInnerRadius;
651  if (i >= m_Params->get_int_param("magnet_cutout_first_scinti"))
652  {
653  z[2] = tan(theta) * (m_ScintiInnerRadius + (m_Params->get_double_param("magnet_cutout_scinti_radius") * cm - m_ScintiInnerRadius));
654  }
655  x[3] = m_ScintiOuterRadius + overhang; // since the tile is tilted, x is not at the outer radius but beyond
656  z[3] = tan(theta) * m_ScintiOuterRadius;
657  // apply gap between scintillators
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.;
662  PHG4OuterHcalDetector::Point_2 leftsidelow(z[0], x[0]);
663  PHG4OuterHcalDetector::Point_2 leftsidehigh(z[1], x[1]);
664  x[0] = m_ScintiInnerRadius - inner_offset;
665  z[0] = x_at_y(leftsidelow, leftsidehigh, x[0]);
666  x[1] = m_ScintiOuterRadius + offset;
667  z[1] = x_at_y(leftsidelow, leftsidehigh, x[1]);
668  PHG4OuterHcalDetector::Point_2 rightsidelow(z[2], x[2]);
669  PHG4OuterHcalDetector::Point_2 rightsidehigh(z[3], x[3]);
670  x[2] = m_ScintiOuterRadius + offset;
671  z[2] = x_at_y(rightsidelow, rightsidehigh, x[2]);
672  x[3] = m_ScintiInnerRadius - inner_offset;
673  z[3] = x_at_y(rightsidelow, rightsidehigh, x[3]);
674  // store corner points of extruded solid we need to subtract from steel
675  if (i == m_Params->get_int_param("magnet_cutout_first_scinti"))
676  {
677  double x0 = m_InnerRadius - (steel_inner_offset - steel_magnet_cutout_x);
678  double z0 = x_at_y(leftsidelow, leftsidehigh, x0);
679  double xpos = m_InnerRadius - steel_offset;
680  zsteelcut[0] = z0;
681  zsteelcut[3] = x_at_y(leftsidelow, leftsidehigh, xpos);
682  }
683  double x2 = m_OuterRadius + steel_offset;
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++)
689  {
690  G4TwoVector v(x[j], z[j]);
691  vertexes.push_back(v);
692  }
693  G4TwoVector zero(0, 0);
694 
695  G4VSolid *scinti = new G4ExtrudedSolid("ScintillatorTile",
696  vertexes,
697  m_ScintiTileThickness + 0.2 * mm,
698  zero, 1.0,
699  zero, 1.0);
700  G4RotationMatrix *rotm = new G4RotationMatrix();
701  rotm->rotateX(-90 * deg);
702  name.str("");
703  name << "scintillator_" << i << "_left";
704  G4VSolid *scinti_tile = new G4IntersectionSolid(name.str(), bigtile, scinti, rotm, G4ThreeVector(-(m_ScintiInnerRadius + m_ScintiOuterRadius) / 2., 0, 0));
705  delete rotm;
706  m_ScintiTilesVec[i + m_NumScintiTiles] = scinti_tile;
707  rotm = new G4RotationMatrix();
708  rotm->rotateX(90 * deg);
709  name.str("");
710  name << "scintillator_" << i << "_right";
711  scinti_tile = new G4IntersectionSolid(name.str(), bigtile, scinti, rotm, G4ThreeVector(-(m_ScintiInnerRadius + m_ScintiOuterRadius) / 2., 0, 0));
712  delete rotm;
713  m_ScintiTilesVec[m_NumScintiTiles - i - 1] = scinti_tile;
714  }
715 #ifdef SCINTITEST
716  for (unsigned int i = 0; i < m_ScintiTilesVec.size(); i++)
717  {
718  if (m_ScintiTilesVec[i])
719  {
720  DisplayVolume(m_ScintiTilesVec[i], hcalenvelope);
721  }
722  }
723 #endif
724 
725  std::vector<G4TwoVector> vertexes;
726  for (int j = 0; j < 4; j++)
727  {
728  if (!isfinite(zsteelcut[j]))
729  {
730  return;
731  }
732  G4TwoVector v(xsteelcut[j], zsteelcut[j]);
733  vertexes.push_back(v);
734  }
735  G4TwoVector zero(0, 0);
736  m_SteelCutoutForMagnetG4Solid = new G4ExtrudedSolid("ScintillatorTile",
737  vertexes,
738  m_ScintiTileThickness + 20 * cm,
739  zero, 1.0,
740  zero, 1.0);
741  return;
742 }
743 
744 G4double
746 {
747  double xret = NAN;
748  double x[2];
749  x[0] = CGAL::to_double(p0.x());
750  x[1] = CGAL::to_double(p1.x());
751  Line_2 l(p0, p1);
752  double newx = fabs(x[0]) + fabs(x[1]);
753  PHG4OuterHcalDetector::Point_2 p0new(-newx, yin);
754  PHG4OuterHcalDetector::Point_2 p1new(newx, yin);
755  Segment_2 seg(p0new, p1new);
756  CGAL::Object result = CGAL::intersection(l, seg);
757  if (const PHG4OuterHcalDetector::Point_2 *ipoint = CGAL::object_cast<PHG4OuterHcalDetector::Point_2>(&result))
758  {
759  xret = CGAL::to_double(ipoint->x());
760  }
761  else
762  {
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;
766  exit(1);
767  }
768  return xret;
769 }
770 
771 G4AssemblyVolume *
773 {
774 #ifdef SCINTITEST
775  ConstructHcalSingleScintillators(hcalenvelope);
776  return nullptr;
777 #endif
778  ConstructHcalSingleScintillators(hcalenvelope);
779  G4AssemblyVolume *assmeblyvol = new G4AssemblyVolume();
780  std::ostringstream name;
781  G4ThreeVector g4vec;
782  double steplimits = m_Params->get_double_param("steplimits") * cm;
783  for (unsigned int i = 0; i < m_ScintiTilesVec.size(); i++)
784  {
785  name.str("");
786  name << m_ScintiLogicNamePrefix << i;
787  G4UserLimits *g4userlimits = nullptr;
788  if (isfinite(steplimits))
789  {
790  g4userlimits = new G4UserLimits(steplimits);
791  }
792  G4LogicalVolume *scinti_tile_logic = new G4LogicalVolume(m_ScintiTilesVec[i], GetDetectorMaterial("G4_POLYSTYRENE"), name.str(), nullptr, nullptr, g4userlimits);
793  m_DisplayAction->AddScintiVolume(scinti_tile_logic);
794  assmeblyvol->AddPlacedVolume(scinti_tile_logic, g4vec, nullptr);
795 
796  //field after burner
797  scinti_tile_logic->SetFieldManager(m_FieldSetup->get_Field_Manager_Gap(), true);
798  }
799  return assmeblyvol;
800 }
801 
803 {
804  // just make sure the parameters make a bit of sense
806  {
807  std::cout << PHWHERE << ": Inner Radius " << m_InnerRadius / cm
808  << " cm larger than Outer Radius " << m_OuterRadius / cm
809  << " cm" << std::endl;
810  gSystem->Exit(1);
811  }
813  {
814  std::cout << PHWHERE << "Scintillator thickness " << m_ScintiTileThickness / cm
815  << " cm larger than scintillator gap " << m_ScintiGap / cm
816  << " cm" << std::endl;
817  gSystem->Exit(1);
818  }
820  {
821  std::cout << PHWHERE << "Scintillator outer radius " << m_ScintiOuterRadius / cm
822  << " cm smaller than scintillator inner radius " << m_ScintiInnerRadius / cm
823  << " cm" << std::endl;
824  gSystem->Exit(1);
825  }
827  {
828  std::cout << PHWHERE << "Scintillator outer radius " << m_ScintiOuterRadius / cm
829  << " cm smaller than inner radius " << m_InnerRadius / cm
830  << " cm" << std::endl;
831  gSystem->Exit(1);
832  }
834  {
835  std::cout << PHWHERE << "Scintillator inner radius " << m_ScintiInnerRadius / cm
836  << " cm larger than inner radius " << m_InnerRadius / cm
837  << " cm" << std::endl;
838  gSystem->Exit(1);
839  }
840  if (m_Params->get_double_param("magnet_cutout_scinti_radius") * cm < m_ScintiInnerRadius)
841  {
842  std::cout << PHWHERE << "Magnet scintillator cutout radius " << m_Params->get_double_param("magnet_cutout_scinti_radius")
843  << " cm smaller than inner scintillator radius " << m_ScintiInnerRadius / cm
844  << " cm" << std::endl;
845  gSystem->Exit(1);
846  }
847  if (m_Params->get_double_param("magnet_cutout_radius") * cm < m_InnerRadius)
848  {
849  std::cout << PHWHERE << "Magnet steel cutout radius " << m_Params->get_double_param("magnet_cutout_radius")
850  << " cm smaller than inner radius " << m_InnerRadius / cm
851  << " cm" << std::endl;
852  gSystem->Exit(1);
853  }
854 
855  return 0;
856 }
857 
859 {
860  int ncross = m_Params->get_int_param("ncross");
861  if (!ncross || isfinite(m_TiltAngle))
862  {
863  // mark ncross parameter as not used
864  m_Params->set_int_param("ncross", 0);
865  return;
866  }
867  if ((isfinite(m_TiltAngle)) && (Verbosity() > 0))
868  {
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;
871  }
872  double mid_radius = m_InnerRadius + (m_OuterRadius - m_InnerRadius) / 2.;
873  double deltaphi = (2 * M_PI / m_NumScintiPlates) * ncross;
874  PHG4OuterHcalDetector::Point_2 pnull(0, 0);
876  PHG4OuterHcalDetector::Point_2 phightmp(1, tan(deltaphi));
878  Circle_2 inner_circle(pin1, pin2, pin3);
879  PHG4OuterHcalDetector::Point_2 pmid1(mid_radius, 0), pmid2(0, mid_radius), pmid3(-mid_radius, 0);
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)
889  {
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))
892  {
893  if (CGAL::to_double(point->first.x()) > 0)
894  {
895  PHG4OuterHcalDetector::Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
896  upperright = pntmp;
897  }
898  }
899  else
900  {
901  std::cout << "CGAL::Object type not pair..." << std::endl;
902  exit(1);
903  }
904  }
905  Line_2 l_right(plow, upperright);
906  res.clear();
908  CGAL::intersection(mid_circle, l_right, std::back_inserter(res));
909  for (iter = res.begin(); iter != res.end(); ++iter)
910  {
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))
913  {
914  if (CGAL::to_double(point->first.x()) > 0)
915  {
916  PHG4OuterHcalDetector::Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
917  midpoint = pntmp;
918  }
919  }
920  else
921  {
922  std::cout << "CGAL::Object type not pair..." << std::endl;
923  exit(1);
924  }
925  }
926  // length left side
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()));
929  // c^2 = a^2+b^2 - 2ab*cos(gamma)
930  // gamma = acos((a^2+b^2=c^2)/2ab
931  double tiltangle = acos((ll * ll + upside * upside - m_InnerRadius * m_InnerRadius) / (2 * ll * upside));
932  tiltangle = tiltangle * rad;
933  m_TiltAngle = copysign(tiltangle, ncross);
934  m_Params->set_double_param("tilt_angle", m_TiltAngle / deg);
935  return;
936 }
937 
938 // check if tilt angle is reasonable - too large, no intersections with inner radius
940 {
941  if (fabs(m_TiltAngle) >= M_PI)
942  {
943  std::cout << PHWHERE << "invalid tilt angle, abs(tilt) >= 90 deg: " << (m_TiltAngle / deg)
944  << std::endl;
945  exit(1);
946  }
947 
948  double mid_radius = m_InnerRadius + (m_OuterRadius - m_InnerRadius) / 2.;
949  PHG4OuterHcalDetector::Point_2 pmid(mid_radius, 0); // center of scintillator
950  double xcoord = 0;
951  double ycoord = mid_radius * tan(m_TiltAngle / rad);
952  PHG4OuterHcalDetector::Point_2 pxnull(xcoord, ycoord);
953  Line_2 s2(pmid, pxnull);
955  Circle_2 inner_circle(sc1, sc2, sc3);
956  std::vector<CGAL::Object> res;
957  CGAL::intersection(inner_circle, s2, std::back_inserter(res));
958  if (res.size() == 0)
959  {
960  std::cout << PHWHERE << " Tilt angle " << (m_TiltAngle / deg)
961  << " too large, no intersection with inner radius" << std::endl;
962  exit(1);
963  }
964  return 0;
965 }
966 
968 {
969  std::cout << "Outer Hcal Detector:" << std::endl;
970  if (what == "ALL" || what == "VOLUME")
971  {
972  std::cout << "Volume Envelope: " << m_VolumeEnvelope / cm / cm / cm << " cm^3" << std::endl;
973  std::cout << "Volume Steel: " << m_VolumeSteel / cm / cm / cm << " cm^3" << std::endl;
974  std::cout << "Volume Scintillator: " << m_VolumeScintillator / cm / cm / cm << " cm^3" << std::endl;
975  std::cout << "Volume Air: " << (m_VolumeEnvelope - m_VolumeSteel - m_VolumeScintillator) / cm / cm / cm << " cm^3" << std::endl;
976  }
977  return;
978 }
979 
981 {
982  auto it = m_ScintiTilePhysVolMap.find(volume);
983  if (it != m_ScintiTilePhysVolMap.end())
984  {
985  return it->second;
986  }
987  std::cout << "could not locate volume " << volume->GetName()
988  << " in Inner Hcal scintillator map" << std::endl;
989  gSystem->Exit(1);
990  // that's dumb but code checkers do not know that gSystem->Exit()
991  // terminates, so using the standard exit() makes them happy
992  exit(1);
993 }
994 
995 // This is dulplicated code, we can get rid of it when we have the code to make towergeom for real data reco.
997 {
998  PHNodeIterator iter(topNode());
999  PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
1000  if (!runNode)
1001  {
1002  std::cout << PHWHERE << "Run Node missing, exiting." << std::endl;
1003  gSystem->Exit(1);
1004  exit(1);
1005  }
1006  PHNodeIterator runIter(runNode);
1007  PHCompositeNode *RunDetNode = dynamic_cast<PHCompositeNode *>(runIter.findFirst("PHCompositeNode", m_SuperDetector));
1008  if (!RunDetNode)
1009  {
1010  RunDetNode = new PHCompositeNode(m_SuperDetector);
1011  runNode->addNode(RunDetNode);
1012  }
1013  m_TowerGeomNodeName = "TOWERGEOM_" + m_SuperDetector;
1014  m_RawTowerGeom = findNode::getClass<RawTowerGeomContainer>(topNode(), m_TowerGeomNodeName);
1015  if (!m_RawTowerGeom)
1016  {
1019  RunDetNode->addNode(newNode);
1020  }
1023  m_RawTowerGeom->set_radius(innerrad);
1024  m_RawTowerGeom->set_thickness(thickness);
1027  double geom_ref_radius = innerrad + thickness / 2.;
1028  double phistart = m_Params->get_double_param("phistart");
1029  if (!std::isfinite(phistart))
1030  {
1031  std::cout << PHWHERE << " phistart is not finite: " << phistart
1032  << ", exiting now (this will crash anyway)" << std::endl;
1033  gSystem->Exit(1);
1034  }
1035  for (int i = 0; i < m_Params->get_int_param(PHG4HcalDefs::n_towers); i++)
1036  {
1037  double phiend = phistart + 2. * M_PI / m_Params->get_int_param(PHG4HcalDefs::n_towers);
1038  std::pair<double, double> range = std::make_pair(phiend, phistart);
1039  phistart = phiend;
1040  m_RawTowerGeom->set_phibounds(i, range);
1041  }
1042  double etalowbound = -m_Params->get_double_param("scinti_eta_coverage_neg");
1043  for (int i = 0; i < m_Params->get_int_param("etabins"); i++)
1044  {
1045  // double etahibound = etalowbound + 2.2 / get_int_param("etabins");
1046  double etahibound = etalowbound +
1047  (m_Params->get_double_param("scinti_eta_coverage_neg") + m_Params->get_double_param("scinti_eta_coverage_pos")) / m_Params->get_int_param("etabins");
1048  std::pair<double, double> range = std::make_pair(etalowbound, etahibound);
1049  m_RawTowerGeom->set_etabounds(i, range);
1050  etalowbound = etahibound;
1051  }
1052  for (int iphi = 0; iphi < m_RawTowerGeom->get_phibins(); iphi++)
1053  {
1054  for (int ieta = 0; ieta < m_RawTowerGeom->get_etabins(); ieta++)
1055  {
1057 
1058  const double x(geom_ref_radius * cos(m_RawTowerGeom->get_phicenter(iphi)));
1059  const double y(geom_ref_radius * sin(m_RawTowerGeom->get_phicenter(iphi)));
1060  const double z(geom_ref_radius / tan(PHG4Utils::get_theta(m_RawTowerGeom->get_etacenter(ieta))));
1061 
1063  if (tg)
1064  {
1065  if (Verbosity() > 0)
1066  {
1067  std::cout << "IHCalDetector::InitRun - Tower geometry " << key << " already exists" << std::endl;
1068  }
1069 
1070  if (fabs(tg->get_center_x() - x) > 1e-4)
1071  {
1072  std::cout << "IHCalDetector::InitRun - Fatal Error - duplicated Tower geometry " << key << " with existing x = " << tg->get_center_x() << " and expected x = " << x
1073  << std::endl;
1074 
1075  return;
1076  }
1077  if (fabs(tg->get_center_y() - y) > 1e-4)
1078  {
1079  std::cout << "IHCalDetector::InitRun - Fatal Error - duplicated Tower geometry " << key << " with existing y = " << tg->get_center_y() << " and expected y = " << y
1080  << std::endl;
1081  return;
1082  }
1083  if (fabs(tg->get_center_z() - z) > 1e-4)
1084  {
1085  std::cout << "IHCalDetector::InitRun - Fatal Error - duplicated Tower geometry " << key << " with existing z= " << tg->get_center_z() << " and expected z = " << z
1086  << std::endl;
1087  return;
1088  }
1089  }
1090  else
1091  {
1092  if (Verbosity() > 0)
1093  {
1094  std::cout << "IHCalDetector::InitRun - building tower geometry " << key << "" << std::endl;
1095  }
1096 
1097  tg = new RawTowerGeomv1(key);
1098 
1099  tg->set_center_x(x);
1100  tg->set_center_y(y);
1101  tg->set_center_z(z);
1103  }
1104  }
1105  }
1106  if (Verbosity() > 0)
1107  {
1109  }
1110 }