Analysis Software
Documentation for sPHENIX simulation 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/PHCompositeNode.h>
15 #include <phool/PHIODataNode.h>
16 #include <phool/PHNode.h> // for PHNode
17 #include <phool/PHNodeIterator.h>
18 #include <phool/PHObject.h> // for PHObject
19 #include <phool/getClass.h>
20 #include <phool/phool.h>
21 #include <phool/recoConsts.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/G4SystemOfUnits.hh>
41 #include <Geant4/G4ThreeVector.hh>
42 #include <Geant4/G4Transform3D.hh>
43 #include <Geant4/G4Tubs.hh>
44 #include <Geant4/G4TwoVector.hh>
45 #include <Geant4/G4Types.hh>
46 #include <Geant4/G4UserLimits.hh>
47 #include <Geant4/G4VSolid.hh>
48 
49 #pragma GCC diagnostic push
50 #pragma GCC diagnostic ignored "-Wshadow"
51 #pragma GCC diagnostic ignored "-Wpedantic"
52 #include <CGAL/Boolean_set_operations_2.h>
53 #include <CGAL/Circular_kernel_intersections.h>
54 #include <CGAL/Exact_circular_kernel_2.h>
55 #include <CGAL/Object.h>
56 #include <CGAL/point_generators_2.h>
57 #pragma GCC diagnostic pop
58 
59 #include <boost/lexical_cast.hpp>
60 #include <boost/tokenizer.hpp>
61 
62 #include <algorithm>
63 #include <cmath>
64 #include <cstdlib>
65 #include <iostream>
66 #include <iterator>
67 #include <sstream>
68 
69 class PHCompositeNode;
70 
71 using Circle_2 = CGAL::Circle_2<PHG4InnerHcalDetector::Circular_k>;
72 using Circular_arc_point_2 = CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>;
73 using Line_2 = CGAL::Line_2<PHG4InnerHcalDetector::Circular_k>;
74 using Segment_2 = CGAL::Segment_2<PHG4InnerHcalDetector::Circular_k>;
75 
76 // there is still a minute problem for very low tilt angles where the scintillator
77 // face touches the boundary instead of the corner, subtracting 1 permille from the total
78 // scintilator length takes care of this
79 static double subtract_from_scinti_x = 0.1 * mm;
80 
82  : PHG4Detector(subsys, Node, dnam)
83  , m_DisplayAction(dynamic_cast<PHG4InnerHcalDisplayAction *>(subsys->GetDisplayAction()))
84  , m_Params(parameters)
85  , m_InnerRadius(m_Params->get_double_param("inner_radius") * cm)
86  , m_OuterRadius(m_Params->get_double_param("outer_radius") * cm)
87  , m_SizeZ(m_Params->get_double_param("size_z") * cm)
88  , m_ScintiTileZ(m_SizeZ)
89  , m_ScintiTileThickness(m_Params->get_double_param("scinti_tile_thickness") * cm)
90  , m_ScintiInnerGap(m_Params->get_double_param("scinti_inner_gap") * cm)
91  , m_ScintiOuterGap(m_Params->get_double_param("scinti_outer_gap") * cm)
92  , m_ScintiOuterRadius(m_Params->get_double_param("scinti_outer_radius") * cm)
93  , m_TiltAngle(m_Params->get_double_param("tilt_angle") * deg)
94  , m_EnvelopeInnerRadius(m_InnerRadius)
95  , m_EnvelopeOuterRadius(m_OuterRadius)
96  , m_EnvelopeZ(m_SizeZ)
97  , m_NumScintiPlates(m_Params->get_int_param(PHG4HcalDefs::scipertwr) * m_Params->get_int_param("n_towers"))
98  , m_NumScintiTilesPos(m_Params->get_int_param("n_scinti_tiles_pos"))
99  , m_NumScintiTilesNeg(m_Params->get_int_param("n_scinti_tiles_neg"))
100  , m_Active(m_Params->get_int_param("active"))
101  , m_AbsorberActive(m_Params->get_int_param("absorberactive"))
102  , m_ScintiLogicNamePrefix("HcalInnerScinti")
103 {
104  // n_scinti_tiles takes precedence
105 
107  if (nTiles <= 0)
108  {
110  }
111  else
112  {
113  m_NumScintiTilesPos = nTiles / 2;
114  m_Params->set_int_param("n_scinti_tiles_pos", nTiles / 2);
115  m_NumScintiTilesNeg = nTiles / 2;
116  m_Params->set_int_param("n_scinti_tiles_neg", nTiles / 2);
117  }
118 
119  // allocate memory for scintillator plates
120  m_ScintiTilesVec.assign(nTiles, static_cast<G4VSolid *>(nullptr));
121 }
122 
124 {
125  delete m_ScintiMotherAssembly;
126 }
127 
128 //_______________________________________________________________
129 //_______________________________________________________________
131 {
132  if (m_AbsorberActive)
133  {
134  if (m_SteelAbsorberPhysVolSet.find(volume) != m_SteelAbsorberPhysVolSet.end())
135  {
136  return -1;
137  }
138  }
139  if (m_Active)
140  {
141  if (m_ScintiTilePhysVolMap.find(volume) != m_ScintiTilePhysVolMap.end())
142  {
143  return 1;
144  }
145  }
146  return 0;
147 }
148 
149 G4VSolid *
150 PHG4InnerHcalDetector::ConstructScintillatorBox(G4LogicalVolume * /*hcalenvelope*/)
151 {
152  double mid_radius = m_InnerRadius + (m_OuterRadius - m_InnerRadius) / 2.;
153  Point_2 p_in_1(mid_radius, 0); // center of scintillator
154 
155  // length of upper edge (middle till outer circle intersect
156  // x/y coordinate of end of center vertical
157  double xcoord = m_ScintiTileThickness / 2. * sin(fabs(m_TiltAngle) / rad) + mid_radius;
158  double ycoord = m_ScintiTileThickness / 2. * cos(fabs(m_TiltAngle) / rad) + 0;
159  Point_2 p_upperedge(xcoord, ycoord);
160  Line_2 s2(p_in_1, p_upperedge); // center vertical
161 
162  Line_2 perp = s2.perpendicular(p_upperedge);
163  Point_2 sc1(m_OuterRadius, 0), sc2(0, m_OuterRadius), sc3(-m_OuterRadius, 0);
164  Circle_2 outer_circle(sc1, sc2, sc3);
165  std::vector<CGAL::Object> res;
166  CGAL::intersection(outer_circle, perp, std::back_inserter(res));
167  Point_2 upperright;
168  std::vector<CGAL::Object>::const_iterator iter;
169  for (iter = res.begin(); iter != res.end(); ++iter)
170  {
171  CGAL::Object obj = *iter;
172  if (const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned> *point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned>>(&obj))
173  {
174  if (CGAL::to_double(point->first.x()) > CGAL::to_double(p_upperedge.x()))
175  {
176  double deltax = CGAL::to_double(point->first.x()) - CGAL::to_double(p_upperedge.x());
177  double deltay = CGAL::to_double(point->first.y()) - CGAL::to_double(p_upperedge.y());
178  // the scintillator is twice as long
179  m_ScintiTileXUpper = sqrt(deltax * deltax + deltay * deltay); //
180  Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
181  upperright = pntmp;
182  }
183  }
184  else
185  {
186  std::cout << "CGAL::Object type not pair..." << std::endl;
187  }
188  }
189  // length of lower edge (middle till inner circle intersect
190  xcoord = mid_radius - m_ScintiTileThickness / 2. * sin(fabs(m_TiltAngle) / rad);
191  ycoord = 0 - m_ScintiTileThickness / 2. * cos(fabs(m_TiltAngle) / rad);
192  Point_2 p_loweredge(xcoord, ycoord);
193  Line_2 s3(p_in_1, p_loweredge);
194  Line_2 l_lower = s3.perpendicular(p_loweredge);
195  Point_2 ic1(m_InnerRadius, 0), ic2(0, m_InnerRadius), ic3(-m_InnerRadius, 0);
196  Circle_2 inner_circle(ic1, ic2, ic3);
197  res.clear();
198  CGAL::intersection(inner_circle, l_lower, std::back_inserter(res));
199  Point_2 lowerleft;
200  // we have 2 intersections - we want the one furthest to the right (largest x). The correct one is
201  // certainly > 0 but the second one depends on the tilt angle and might also be > 0
202  double minx = 0;
203  for (iter = res.begin(); iter != res.end(); ++iter)
204  {
205  CGAL::Object obj = *iter;
206  if (const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned> *point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned>>(&obj))
207  {
208  if (CGAL::to_double(point->first.x()) > minx)
209  {
210  minx = CGAL::to_double(point->first.x());
211  double deltax = CGAL::to_double(point->first.x()) - CGAL::to_double(p_loweredge.x());
212  double deltay = CGAL::to_double(point->first.y()) - CGAL::to_double(p_loweredge.y());
213  m_ScintiTileXLower = sqrt(deltax * deltax + deltay * deltay);
214  Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
215  lowerleft = pntmp;
216  }
217  }
218  }
221  G4VSolid *scintibox = new G4Box("ScintiTile", m_ScintiTileX / 2., m_ScintiTileThickness / 2., m_ScintiTileZ / 2.);
222  m_VolumeScintillator = scintibox->GetCubicVolume() * m_NumScintiPlates;
223  return scintibox;
224 }
225 
226 G4VSolid *
227 PHG4InnerHcalDetector::ConstructSteelPlate(G4LogicalVolume * /*hcalenvelope*/)
228 {
229  // calculate steel plate on top of the scinti box. Lower edge is the upper edge of
230  // the scintibox + 1/2 the airgap
231  double mid_radius = m_InnerRadius + (m_OuterRadius - m_InnerRadius) / 2.;
232  // first the lower edge, just like the scinti box, just add the air gap
233  // and calculate intersection of edge with inner and outer radius.
234  Point_2 p_in_1(mid_radius, 0); // center of lower scintillator
235  double angle_mid_scinti = M_PI / 2. + m_TiltAngle / rad;
236  double xcoord = m_ScintiInnerGap / 2. * cos(angle_mid_scinti / rad) + mid_radius;
237  double ycoord = m_ScintiInnerGap / 2. * sin(angle_mid_scinti / rad) + 0;
238  Point_2 p_loweredge(xcoord, ycoord);
239  Line_2 s2(p_in_1, p_loweredge); // center vertical
240  Line_2 perp = s2.perpendicular(p_loweredge); // that is the lower edge of the steel plate
241  Point_2 sc1(m_InnerRadius, 0), sc2(0, m_InnerRadius), sc3(-m_InnerRadius, 0);
242  Circle_2 inner_circle(sc1, sc2, sc3);
243  std::vector<CGAL::Object> res;
244  CGAL::intersection(inner_circle, perp, std::back_inserter(res));
245  Point_2 lowerleft;
246  std::vector<CGAL::Object>::const_iterator iter;
247  for (iter = res.begin(); iter != res.end(); ++iter)
248  {
249  CGAL::Object obj = *iter;
250  if (const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned> *point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned>>(&obj))
251  {
252  if (CGAL::to_double(point->first.x()) > 0)
253  {
254  Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
255  lowerleft = pntmp;
256  }
257  }
258  else
259  {
260  std::cout << "CGAL::Object type not pair..." << std::endl;
261  }
262  }
263 
264  double xcoord2 = m_ScintiOuterGap / 2. * cos(angle_mid_scinti / rad) + mid_radius;
265  double ycoord2 = m_ScintiOuterGap / 2. * sin(angle_mid_scinti / rad) + 0;
266  Point_2 p_loweredge2(xcoord2, ycoord2);
267  Line_2 s2_2(p_in_1, p_loweredge2); // center vertical
268  Line_2 perp2 = s2_2.perpendicular(p_loweredge2); // that is the lower edge of the steel plate
269 
270  Point_2 so1(m_OuterRadius, 0), so2(0, m_OuterRadius), so3(-m_OuterRadius, 0);
271  Circle_2 outer_circle(so1, so2, so3);
272  res.clear(); // just clear the content from the last intersection search
273  CGAL::intersection(outer_circle, perp2, std::back_inserter(res));
274  Point_2 lowerright;
275  for (iter = res.begin(); iter != res.end(); ++iter)
276  {
277  CGAL::Object obj = *iter;
278  if (const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned> *point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned>>(&obj))
279  {
280  if (CGAL::to_double(point->first.x()) > CGAL::to_double(p_loweredge.x()))
281  {
282  Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
283  lowerright = pntmp;
284  }
285  }
286  else
287  {
288  std::cout << "CGAL::Object type not pair..." << std::endl;
289  }
290  }
291  // now we have the lower left and rigth corner, now find the upper edge
292  // find the center of the upper scintilator
293 
294  double phi_midpoint = 2 * M_PI / m_NumScintiPlates;
295  double xmidpoint = cos(phi_midpoint) * mid_radius;
296  double ymidpoint = sin(phi_midpoint) * mid_radius;
297  // angle of perp line at center of scintillator
298  angle_mid_scinti = (M_PI / 2. - phi_midpoint) - (M_PI / 2. + m_TiltAngle / rad);
299  double xcoordup = xmidpoint - m_ScintiInnerGap / 2. * sin(angle_mid_scinti / rad);
300  double ycoordup = ymidpoint - m_ScintiInnerGap / 2. * cos(angle_mid_scinti / rad);
301  Point_2 upperleft;
302  Point_2 upperright;
303  Point_2 mid_upperscint(xmidpoint, ymidpoint);
304  Point_2 p_upperedge(xcoordup, ycoordup);
305  Line_2 sup(mid_upperscint, p_upperedge); // center vertical
306  Line_2 perpA = sup.perpendicular(p_upperedge); // that is the upper edge of the steel plate
307  Point_2 sc1A(m_InnerRadius, 0), sc2A(0, m_InnerRadius), sc3A(-m_InnerRadius, 0);
308  Circle_2 inner_circleA(sc1A, sc2A, sc3A);
309  std::vector<CGAL::Object> resA;
310  CGAL::intersection(inner_circleA, perpA, std::back_inserter(resA));
311  std::vector<CGAL::Object>::const_iterator iterA;
312  double pxmax = 0.;
313  for (iterA = resA.begin(); iterA != resA.end(); ++iterA)
314  {
315  CGAL::Object obj = *iterA;
316  if (const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned> *point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned>>(&obj))
317  {
318  if (CGAL::to_double(point->first.x()) > pxmax)
319  {
320  pxmax = CGAL::to_double(point->first.x());
321  Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
322  upperleft = pntmp;
323  }
324  }
325  else
326  {
327  std::cout << "CGAL::Object type not pair..." << std::endl;
328  }
329  }
330 
331  double xcoordup2 = xmidpoint - m_ScintiOuterGap / 2. * sin(angle_mid_scinti / rad);
332  double ycoordup2 = ymidpoint - m_ScintiOuterGap / 2. * cos(angle_mid_scinti / rad);
333  Point_2 p_upperedge2(xcoordup2, ycoordup2);
334  Line_2 sup2(mid_upperscint, p_upperedge2); // center vertical
335  Line_2 perpA2 = sup2.perpendicular(p_upperedge2); // that is the upper edge of the steel plate
336 
337  Point_2 so1A(m_OuterRadius, 0), so2A(0, m_OuterRadius), so3A(-m_OuterRadius, 0);
338  Circle_2 outer_circleA(so1A, so2A, so3A);
339  resA.clear(); // just clear the content from the last intersection search
340  CGAL::intersection(outer_circleA, perpA2, std::back_inserter(resA));
341  for (iterA = resA.begin(); iterA != resA.end(); ++iterA)
342  {
343  CGAL::Object obj = *iterA;
344  if (const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned> *point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned>>(&obj))
345  {
346  if (CGAL::to_double(point->first.x()) > CGAL::to_double(p_loweredge.x()))
347  {
348  Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
349  upperright = pntmp;
350  }
351  }
352  else
353  {
354  std::cout << "CGAL::Object type not pair..." << std::endl;
355  }
356  }
357  // the left corners are on a secant with the inner boundary, they need to be shifted
358  // to be a tangent at the center
359  ShiftSecantToTangent(lowerleft, upperleft, upperright, lowerright);
360  G4TwoVector v1(CGAL::to_double(upperleft.x()), CGAL::to_double(upperleft.y()));
361  G4TwoVector v2(CGAL::to_double(upperright.x()), CGAL::to_double(upperright.y()));
362  G4TwoVector v3(CGAL::to_double(lowerright.x()), CGAL::to_double(lowerright.y()));
363  G4TwoVector v4(CGAL::to_double(lowerleft.x()), CGAL::to_double(lowerleft.y()));
364  std::vector<G4TwoVector> vertexes;
365  vertexes.push_back(v1);
366  vertexes.push_back(v2);
367  vertexes.push_back(v3);
368  vertexes.push_back(v4);
369  G4TwoVector zero(0, 0);
370  G4VSolid *steel_plate = new G4ExtrudedSolid("SteelPlate",
371  vertexes,
372  m_SizeZ / 2.0,
373  zero, 1.0,
374  zero, 1.0);
375 
376  // DisplayVolume(steel_plate, hcalenvelope);
377  m_VolumeSteel = steel_plate->GetCubicVolume() * m_NumScintiPlates;
378  return steel_plate;
379 }
380 
381 void PHG4InnerHcalDetector::ShiftSecantToTangent(Point_2 &lowleft, Point_2 &upleft, Point_2 &upright, Point_2 &lowright)
382 {
383  Line_2 secant(lowleft, upleft);
384  Segment_2 upedge(upleft, upright);
385  Segment_2 lowedge(lowleft, lowright);
386  double xmid = (CGAL::to_double(lowleft.x()) + CGAL::to_double(upleft.x())) / 2.;
387  double ymid = (CGAL::to_double(lowleft.y()) + CGAL::to_double(upleft.y())) / 2.;
388  Point_2 midpoint(xmid, ymid);
389  Line_2 sekperp = secant.perpendicular(midpoint);
390  Point_2 sc1(m_InnerRadius, 0), sc2(0, m_InnerRadius), sc3(-m_InnerRadius, 0);
391  Circle_2 inner_circle(sc1, sc2, sc3);
392  std::vector<CGAL::Object> res;
393  CGAL::intersection(inner_circle, sekperp, std::back_inserter(res));
394  std::vector<CGAL::Object>::const_iterator iter;
395  double pxmax = 0.;
396  Point_2 tangtouch;
397  for (iter = res.begin(); iter != res.end(); ++iter)
398  {
399  CGAL::Object obj = *iter;
400  if (const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned> *point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned>>(&obj))
401  {
402  if (CGAL::to_double(point->first.x()) > pxmax)
403  {
404  pxmax = CGAL::to_double(point->first.x());
405  Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
406  tangtouch = pntmp;
407  }
408  }
409  else
410  {
411  std::cout << "CGAL::Object type not pair..." << std::endl;
412  }
413  }
414  Line_2 leftside = sekperp.perpendicular(tangtouch);
415  CGAL::Object result = CGAL::intersection(upedge, leftside);
416  if (const Point_2 *ipoint = CGAL::object_cast<Point_2>(&result))
417  {
418  upleft = *ipoint;
419  }
420  result = CGAL::intersection(lowedge, leftside);
421  if (const Point_2 *ipoint = CGAL::object_cast<Point_2>(&result))
422  {
423  lowleft = *ipoint;
424  }
425  return;
426 }
427 
428 // Construct the envelope and the call the
429 // actual inner hcal construction
430 void PHG4InnerHcalDetector::ConstructMe(G4LogicalVolume *logicWorld)
431 {
433  G4Material *Air = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
434  G4VSolid *hcal_envelope_cylinder = new G4Tubs("InnerHcal_envelope_solid", m_EnvelopeInnerRadius, m_EnvelopeOuterRadius, m_EnvelopeZ / 2., 0, 2 * M_PI);
435  m_VolumeEnvelope = hcal_envelope_cylinder->GetCubicVolume();
436  G4LogicalVolume *hcal_envelope_log = new G4LogicalVolume(hcal_envelope_cylinder, Air, G4String("Hcal_envelope"), nullptr, nullptr, nullptr);
437 
438  G4RotationMatrix hcal_rotm;
439  hcal_rotm.rotateX(m_Params->get_double_param("rot_x") * deg);
440  hcal_rotm.rotateY(m_Params->get_double_param("rot_y") * deg);
441  hcal_rotm.rotateZ(m_Params->get_double_param("rot_z") * deg);
442  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, false, false, OverlapCheck());
443  m_DisplayAction->SetMyTopVolume(mothervol);
444  ConstructInnerHcal(hcal_envelope_log);
445  std::vector<G4VPhysicalVolume *>::iterator it = m_ScintiMotherAssembly->GetVolumesIterator();
446  for (unsigned int i = 0; i < m_ScintiMotherAssembly->TotalImprintedVolumes(); i++)
447  {
448  // G4AssemblyVolumes naming convention:
449  // av_WWW_impr_XXX_YYY_ZZZ
450  // where:
451 
452  // WWW - assembly volume instance number
453  // XXX - assembly volume imprint number
454  // YYY - the name of the placed logical volume
455  // ZZZ - the logical volume index inside the assembly volume
456  // e.g. av_1_impr_82_HcalInnerScinti_11_pv_11
457  // 82 the number of the scintillator mother volume
458  // HcalInnerScinti_11: name of scintillator slat
459  // 11: number of scintillator slat logical volume
460  // use boost tokenizer to separate the _, then take value
461  // after "impr" for mother volume and after "pv" for scintillator slat
462  // use boost lexical cast for string -> int conversion
463  // the CopyNo is the mother volume + scinti id
464  // so we can use the CopyNo rather than decoding the string further
465  // looking for "pv"
466  boost::char_separator<char> sep("_");
467  boost::tokenizer<boost::char_separator<char>> tok((*it)->GetName(), sep);
468  boost::tokenizer<boost::char_separator<char>>::const_iterator tokeniter;
469  for (tokeniter = tok.begin(); tokeniter != tok.end(); ++tokeniter)
470  {
471  if (*tokeniter == "impr")
472  {
473  ++tokeniter;
474  if (tokeniter != tok.end())
475  {
476  int layer_id = boost::lexical_cast<int>(*tokeniter);
477  // check detector description, for assemblyvolumes it is not possible
478  // to give the first volume id=0, so they go from id=1 to id=n.
479  // I am not going to start with fortran again - our indices start
480  // at zero, id=0 to id=n-1. So subtract one here
481  int tower_id = (*it)->GetCopyNo() - layer_id;
482  layer_id--;
483  std::pair<int, int> layer_twr = std::make_pair(layer_id, tower_id);
484  m_ScintiTilePhysVolMap.insert(std::pair<G4VPhysicalVolume *, std::pair<int, int>>(*it, layer_twr));
485  if (layer_id < 0 || layer_id >= m_NumScintiPlates)
486  {
487  std::cout << "invalid scintillator row " << layer_id
488  << ", valid range 0 < row < " << m_NumScintiPlates << std::endl;
489  gSystem->Exit(1);
490  }
491  }
492  else
493  {
494  std::cout << PHWHERE << " Error parsing " << (*it)->GetName()
495  << " for mother volume number " << std::endl;
496  gSystem->Exit(1);
497  }
498  break;
499  }
500  }
501  ++it;
502  }
503  if (!m_Params->get_int_param("saveg4hit")) AddGeometryNode();
504  return;
505 }
506 
507 int PHG4InnerHcalDetector::ConstructInnerHcal(G4LogicalVolume *hcalenvelope)
508 {
510  SetTiltViaNcross(); // if number of crossings is set, use it to determine tilt
511  CheckTiltAngle(); // die if the tilt angle is out of range
512  G4VSolid *steel_plate = ConstructSteelPlate(hcalenvelope);
513  G4LogicalVolume *steel_logical = new G4LogicalVolume(steel_plate, GetDetectorMaterial(m_Params->get_string_param("material")), "HcalInnerSteelPlate", nullptr, nullptr, nullptr);
514  m_DisplayAction->AddSteelVolume(steel_logical);
516  double phi = 0;
517  double deltaphi = 2 * M_PI / m_NumScintiPlates;
518  std::ostringstream name;
519  double middlerad = m_OuterRadius - (m_OuterRadius - m_InnerRadius) / 2.;
520  // for shorter scintillators we have to shift by 1/2 the "missing length"
521  double scinti_tile_orig_length = m_ScintiTileXUpper + m_ScintiTileXLower - subtract_from_scinti_x;
522  double shiftslat = fabs(m_ScintiTileXLower - m_ScintiTileXUpper) / 2. + (scinti_tile_orig_length - m_ScintiTileX) / 2.;
523  // calculate phi offset (copied from code inside following loop):
524  // first get the center point (phi=0) so it's middlerad/0
525  // then shift the scintillator center as documented in loop
526  // then
527  // for positive tilt angles we need the lower left corner of the scintillator
528  // for negative tilt angles we need the upper right corner of the scintillator
529  // as it turns out the code uses the middle of the face of the scintillator
530  // as reference, if this is a problem the code needs to be modified to
531  // actually calculate the corner (but the math of the construction is that
532  // the middle of the scintillator sits at zero)
533  double xp = cos(phi) * middlerad;
534  double yp = sin(phi) * middlerad;
535  xp -= cos((-m_TiltAngle) / rad - phi) * shiftslat;
536  yp += sin((-m_TiltAngle) / rad - phi) * shiftslat;
537  if (m_TiltAngle > 0)
538  {
539  double xo = xp - (m_ScintiTileX / 2.) * cos(m_TiltAngle / rad);
540  double yo = yp - (m_ScintiTileX / 2.) * sin(m_TiltAngle / rad);
541  phi = -atan(yo / xo);
542  }
543  else if (m_TiltAngle < 0)
544  {
545  double xo = xp + (m_ScintiTileX / 2.) * cos(m_TiltAngle / rad);
546  double yo = yp + (m_ScintiTileX / 2.) * sin(m_TiltAngle / rad);
547  phi = -atan(yo / xo);
548  }
549  // else (for m_TiltAngle = 0) phi stays zero
550  for (int i = 0; i < m_NumScintiPlates; i++)
551  {
552  G4RotationMatrix *Rot = new G4RotationMatrix();
553  double ypos = sin(phi) * middlerad;
554  double xpos = cos(phi) * middlerad;
555  // the center of the scintillator is not the center of the inner hcal
556  // but depends on the tilt angle. Therefore we need to shift
557  // the center from the mid point
558  ypos += sin((-m_TiltAngle) / rad - phi) * shiftslat;
559  xpos -= cos((-m_TiltAngle) / rad - phi) * shiftslat;
560  Rot->rotateZ(phi * rad + m_TiltAngle);
561  G4ThreeVector g4vec(xpos, ypos, 0);
562  // great the MakeImprint always adds 1 to the copy number and 0 has a
563  // special meaning (which then also adds 1). Basically our volume names
564  // will start at 1 instead of 0 and there is nothing short of patching this
565  // method. I'll take care of this in the decoding of the volume name
566  // AAAAAAARHGS
567  m_ScintiMotherAssembly->MakeImprint(hcalenvelope, g4vec, Rot, i, OverlapCheck());
568  delete Rot;
569  Rot = new G4RotationMatrix();
570  Rot->rotateZ(-phi * rad);
571  name.str("");
572  name << "InnerHcalSteel_" << i;
573  m_SteelAbsorberPhysVolSet.insert(new G4PVPlacement(Rot, G4ThreeVector(0, 0, 0), steel_logical, name.str(), hcalenvelope, false, i, OverlapCheck()));
574  phi += deltaphi;
575  }
576  return 0;
577 }
578 
579 // split the big scintillator into tiles covering eta ranges
580 // since they are tilted it is not a straightforward theta cut
581 // it starts at a given eta at the inner radius but the outer radius needs adjusting
583 {
584  G4VSolid *bigtile = ConstructScintillatorBox(hcalenvelope);
585 
586  // eta->theta
587  // scinti_eta_coverage takes precedence
588  G4double eta_cov = m_Params->get_double_param("scinti_eta_coverage");
589  if (eta_cov > 0)
590  {
591  m_Params->set_double_param("scinti_eta_coverage_pos", eta_cov);
592  m_Params->set_double_param("scinti_eta_coverage_neg", eta_cov);
593  }
594  G4double delta_eta_pos = m_Params->get_double_param("scinti_eta_coverage_pos") / m_NumScintiTilesPos;
595  G4double delta_eta_neg = m_Params->get_double_param("scinti_eta_coverage_neg") / m_NumScintiTilesNeg;
596 
597  G4double eta = 0;
598  G4double theta;
599  G4double x[4];
600  G4double z[4];
601  std::ostringstream name;
602  double overhang = (m_ScintiTileX - (m_OuterRadius - m_InnerRadius)) / 2.;
603  double offset = 1 * cm + overhang; // add 1cm to make sure the G4ExtrudedSolid
604  // is larger than the tile so we do not have
605  // funny edge effects when overlapping vols
606  double scinti_gap_neighbor = m_Params->get_double_param("scinti_gap_neighbor") * cm;
607 
608  // Positive first, then negative
609 
610  for (int i = 0; i < m_NumScintiTilesPos; i++)
611  {
612  theta = M_PI / 2 - PHG4Utils::get_theta(eta); // theta = 90 for eta=0
613  x[0] = m_InnerRadius - overhang;
614  z[0] = tan(theta) * m_InnerRadius;
615  x[1] = m_OuterRadius + overhang; // since the tile is tilted, x is not at the outer radius but beyond
616  z[1] = tan(theta) * m_OuterRadius;
617  eta += delta_eta_pos;
618  theta = M_PI / 2 - PHG4Utils::get_theta(eta); // theta = 90 for eta=0
619  x[2] = m_InnerRadius - overhang;
620  z[2] = tan(theta) * m_InnerRadius;
621  x[3] = m_OuterRadius + overhang; // since the tile is tilted, x is not at the outer radius but beyond
622  z[3] = tan(theta) * m_OuterRadius;
623  // apply gap between scintillators
624  z[0] += scinti_gap_neighbor / 2.;
625  z[1] += scinti_gap_neighbor / 2.;
626  z[2] -= scinti_gap_neighbor / 2.;
627  z[3] -= scinti_gap_neighbor / 2.;
628  Point_2 leftsidelow(z[0], x[0]);
629  Point_2 leftsidehigh(z[1], x[1]);
630  x[0] = m_InnerRadius - offset;
631  z[0] = x_at_y(leftsidelow, leftsidehigh, x[0]);
632  x[1] = m_OuterRadius + offset;
633  z[1] = x_at_y(leftsidelow, leftsidehigh, x[1]);
634  Point_2 rightsidelow(z[2], x[2]);
635  Point_2 rightsidehigh(z[3], x[3]);
636  x[2] = m_OuterRadius + offset;
637  z[2] = x_at_y(rightsidelow, rightsidehigh, x[2]);
638  x[3] = m_InnerRadius - offset;
639  z[3] = x_at_y(rightsidelow, rightsidehigh, x[3]);
640 
641  std::vector<G4TwoVector> vertexes;
642  for (int j = 0; j < 4; j++)
643  {
644  G4TwoVector v(x[j], z[j]);
645  vertexes.push_back(v);
646  }
647  G4TwoVector zero(0, 0);
648 
649  G4VSolid *scinti = new G4ExtrudedSolid("ScintillatorTile",
650  vertexes,
651  m_ScintiTileThickness + 0.2 * mm,
652  zero, 1.0,
653  zero, 1.0);
654  G4RotationMatrix *rotm = new G4RotationMatrix();
655  rotm->rotateX(-90 * deg);
656  name.str("");
657  name << "scintillator_" << i << "_left";
658  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));
659  delete rotm;
660  m_ScintiTilesVec[i + m_NumScintiTilesNeg] = scinti_tile;
661  }
662 
663  eta = 0.0; // reset
664  for (int i = 0; i < m_NumScintiTilesNeg; i++)
665  {
666  theta = M_PI / 2 - PHG4Utils::get_theta(eta); // theta = 90 for eta=0
667  x[0] = m_InnerRadius - overhang;
668  z[0] = tan(theta) * m_InnerRadius;
669  x[1] = m_OuterRadius + overhang; // since the tile is tilted, x is not at the outer radius but beyond
670  z[1] = tan(theta) * m_OuterRadius;
671  eta += delta_eta_neg;
672  theta = M_PI / 2 - PHG4Utils::get_theta(eta); // theta = 90 for eta=0
673  x[2] = m_InnerRadius - overhang;
674  z[2] = tan(theta) * m_InnerRadius;
675  x[3] = m_OuterRadius + overhang; // since the tile is tilted, x is not at the outer radius but beyond
676  z[3] = tan(theta) * m_OuterRadius;
677  // apply gap between scintillators
678  z[0] += scinti_gap_neighbor / 2.;
679  z[1] += scinti_gap_neighbor / 2.;
680  z[2] -= scinti_gap_neighbor / 2.;
681  z[3] -= scinti_gap_neighbor / 2.;
682  Point_2 leftsidelow(z[0], x[0]);
683  Point_2 leftsidehigh(z[1], x[1]);
684  x[0] = m_InnerRadius - offset;
685  z[0] = x_at_y(leftsidelow, leftsidehigh, x[0]);
686  x[1] = m_OuterRadius + offset;
687  z[1] = x_at_y(leftsidelow, leftsidehigh, x[1]);
688  Point_2 rightsidelow(z[2], x[2]);
689  Point_2 rightsidehigh(z[3], x[3]);
690  x[2] = m_OuterRadius + offset;
691  z[2] = x_at_y(rightsidelow, rightsidehigh, x[2]);
692  x[3] = m_InnerRadius - offset;
693  z[3] = x_at_y(rightsidelow, rightsidehigh, x[3]);
694 
695  std::vector<G4TwoVector> vertexes;
696  for (int j = 0; j < 4; j++)
697  {
698  G4TwoVector v(x[j], z[j]);
699  vertexes.push_back(v);
700  }
701  G4TwoVector zero(0, 0);
702 
703  G4VSolid *scinti = new G4ExtrudedSolid("ScintillatorTile",
704  vertexes,
705  m_ScintiTileThickness + 0.2 * mm,
706  zero, 1.0,
707  zero, 1.0);
708 
709  G4RotationMatrix *rotm = new G4RotationMatrix();
710  rotm->rotateX(90 * deg);
711  name.str("");
712  name << "scintillator_" << i << "_right";
713  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));
714  m_ScintiTilesVec[m_NumScintiTilesNeg - i - 1] = scinti_tile;
715  delete rotm;
716  }
717 
718  // for (unsigned int i=0; i<m_ScintiTilesVec.size(); i++)
719  // {
720  // if (m_ScintiTilesVec[i])
721  // {
722  // DisplayVolume(m_ScintiTilesVec[i],hcalenvelope );
723  // }
724  // }
725 
726  return;
727 }
728 
729 double
731 {
732  double xret = NAN;
733  double x[2];
734  double y[2];
735  x[0] = CGAL::to_double(p0.x());
736  y[0] = CGAL::to_double(p0.y());
737  x[1] = CGAL::to_double(p1.x());
738  y[1] = CGAL::to_double(p1.y());
739  Line_2 l(p0, p1);
740  double newx = fabs(x[0]) + fabs(x[1]);
741  Point_2 p0new(-newx, yin);
742  Point_2 p1new(newx, yin);
743  Segment_2 seg(p0new, p1new);
744  CGAL::Object result = CGAL::intersection(l, seg);
745  if (const Point_2 *ipoint = CGAL::object_cast<Point_2>(&result))
746  {
747  xret = CGAL::to_double(ipoint->x());
748  }
749  else
750  {
751  std::cout << PHWHERE << " failed for y = " << y << std::endl;
752  std::cout << "p0(x): " << CGAL::to_double(p0.x()) << ", p0(y): " << CGAL::to_double(p0.y()) << std::endl;
753  std::cout << "p1(x): " << CGAL::to_double(p1.x()) << ", p1(y): " << CGAL::to_double(p1.y()) << std::endl;
754  exit(1);
755  }
756  return xret;
757 }
758 
759 G4AssemblyVolume *
761 {
762  ConstructHcalSingleScintillators(hcalenvelope);
763  G4AssemblyVolume *assmeblyvol = new G4AssemblyVolume();
764  std::ostringstream name;
765  G4ThreeVector g4vec;
766 
767  double steplimits = m_Params->get_double_param("steplimits") * cm;
768  for (unsigned int i = 0; i < m_ScintiTilesVec.size(); i++)
769  {
770  name.str("");
771  name << m_ScintiLogicNamePrefix << i;
772  G4UserLimits *g4userlimits = nullptr;
773  if (isfinite(steplimits))
774  {
775  g4userlimits = new G4UserLimits(steplimits);
776  }
777  G4LogicalVolume *scinti_tile_logic = new G4LogicalVolume(m_ScintiTilesVec[i], GetDetectorMaterial("G4_POLYSTYRENE"), name.str(), nullptr, nullptr, g4userlimits);
778  m_DisplayAction->AddScintiVolume(scinti_tile_logic);
779  assmeblyvol->AddPlacedVolume(scinti_tile_logic, g4vec, nullptr);
780  }
781  return assmeblyvol;
782 }
783 
784 // check if tilt angle is reasonable - too large, no intersections with inner radius
786 {
787  if (fabs(m_TiltAngle / rad) >= M_PI)
788  {
789  std::cout << PHWHERE << "invalid tilt angle, abs(tilt) >= 90 deg: " << (m_TiltAngle / deg)
790  << std::endl;
791  exit(1);
792  }
793 
794  double mid_radius = m_InnerRadius + (m_OuterRadius - m_InnerRadius) / 2.;
795  Point_2 pmid(mid_radius, 0); // center of scintillator
796  double xcoord = 0;
797  double ycoord = mid_radius * tan(m_TiltAngle / rad);
798  Point_2 pxnull(xcoord, ycoord);
799  Line_2 s2(pmid, pxnull);
800  Point_2 sc1(m_InnerRadius, 0), sc2(0, m_InnerRadius), sc3(-m_InnerRadius, 0);
801  Circle_2 inner_circle(sc1, sc2, sc3);
802  std::vector<CGAL::Object> res;
803  CGAL::intersection(inner_circle, s2, std::back_inserter(res));
804  if (res.size() == 0)
805  {
806  std::cout << PHWHERE << " Tilt angle " << (m_TiltAngle / deg)
807  << " too large, no intersection with inner radius" << std::endl;
808  exit(1);
809  }
810  return 0;
811 }
812 
814 {
815  // just make sure the parameters make a bit of sense
817  {
818  std::cout << PHWHERE << ": Inner Radius " << m_InnerRadius / cm
819  << " cm larger than Outer Radius " << m_OuterRadius / cm
820  << " cm" << std::endl;
821  gSystem->Exit(1);
822  }
824  {
825  std::cout << PHWHERE << "Scintillator thickness " << m_ScintiTileThickness / cm
826  << " cm larger than scintillator inner gap " << m_ScintiInnerGap / cm
827  << " cm" << std::endl;
828  gSystem->Exit(1);
829  }
831  {
832  std::cout << PHWHERE << "Scintillator outer radius " << m_ScintiOuterRadius / cm
833  << " cm smaller than inner radius " << m_InnerRadius / cm
834  << " cm" << std::endl;
835  gSystem->Exit(1);
836  }
837  return 0;
838 }
839 
841 {
842  // if tilt angle is set it takes precedence
843  // this happens during readback where the tilt angle was set
844  // in this method
845  int ncross = m_Params->get_int_param("ncross");
846  if (!ncross || isfinite(m_TiltAngle))
847  {
848  // flag ncrossing as not used
849  m_Params->set_int_param("ncross", 0);
850  return;
851  }
852  if ((isfinite(m_TiltAngle)) && (Verbosity() > 0))
853  {
854  std::cout << "both number of crossings and tilt angle are set" << std::endl;
855  std::cout << "using number of crossings to determine tilt angle" << std::endl;
856  }
857  double mid_radius = m_InnerRadius + (m_OuterRadius - m_InnerRadius) / 2.;
858  double deltaphi = (2 * M_PI / m_NumScintiPlates) * ncross;
859  Point_2 pnull(0, 0);
860  Point_2 plow(m_InnerRadius, 0);
861  Point_2 phightmp(1, tan(deltaphi));
862  Point_2 pin1(m_InnerRadius, 0), pin2(0, m_InnerRadius), pin3(-m_InnerRadius, 0);
863  Circle_2 inner_circle(pin1, pin2, pin3);
864  Point_2 pmid1(mid_radius, 0), pmid2(0, mid_radius), pmid3(-mid_radius, 0);
865  Circle_2 mid_circle(pmid1, pmid2, pmid3);
866  Point_2 pout1(m_OuterRadius, 0), pout2(0, m_OuterRadius), pout3(-m_OuterRadius, 0);
867  Circle_2 outer_circle(pout1, pout2, pout3);
868  Line_2 l_up(pnull, phightmp);
869  std::vector<CGAL::Object> res;
870  CGAL::intersection(outer_circle, l_up, std::back_inserter(res));
871  Point_2 upperright;
872  std::vector<CGAL::Object>::const_iterator iter;
873  for (iter = res.begin(); iter != res.end(); ++iter)
874  {
875  CGAL::Object obj = *iter;
876  if (const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned> *point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned>>(&obj))
877  {
878  if (CGAL::to_double(point->first.x()) > 0)
879  {
880  Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
881  upperright = pntmp;
882  }
883  }
884  else
885  {
886  std::cout << "CGAL::Object type not pair..." << std::endl;
887  exit(1);
888  }
889  }
890  Line_2 l_right(plow, upperright);
891  res.clear();
892  Point_2 midpoint;
893  CGAL::intersection(mid_circle, l_right, std::back_inserter(res));
894  for (iter = res.begin(); iter != res.end(); ++iter)
895  {
896  CGAL::Object obj = *iter;
897  if (const std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned> *point = CGAL::object_cast<std::pair<CGAL::Circular_arc_point_2<PHG4InnerHcalDetector::Circular_k>, unsigned>>(&obj))
898  {
899  if (CGAL::to_double(point->first.x()) > 0)
900  {
901  Point_2 pntmp(CGAL::to_double(point->first.x()), CGAL::to_double(point->first.y()));
902  midpoint = pntmp;
903  }
904  }
905  else
906  {
907  std::cout << "CGAL::Object type not pair..." << std::endl;
908  exit(1);
909  }
910  }
911  // length left side
912  double ll = sqrt((CGAL::to_double(midpoint.x()) - m_InnerRadius) * (CGAL::to_double(midpoint.x()) - m_InnerRadius) + CGAL::to_double(midpoint.y()) * CGAL::to_double(midpoint.y()));
913  double upside = sqrt(CGAL::to_double(midpoint.x()) * CGAL::to_double(midpoint.x()) + CGAL::to_double(midpoint.y()) * CGAL::to_double(midpoint.y()));
914  // c^2 = a^2+b^2 - 2ab*cos(gamma)
915  // gamma = acos((a^2+b^2=c^2)/2ab
916  double tiltangle = acos((ll * ll + upside * upside - m_InnerRadius * m_InnerRadius) / (2 * ll * upside));
917  tiltangle = tiltangle * rad;
918  m_TiltAngle = copysign(tiltangle, ncross);
919  m_Params->set_double_param("tilt_angle", m_TiltAngle / deg);
920  return;
921 }
922 
924 {
925  std::cout << "Inner Hcal Detector:" << std::endl;
926  if (what == "ALL" || what == "VOLUME")
927  {
928  std::cout << "Volume Envelope: " << m_VolumeEnvelope / cm3 << " cm^3" << std::endl;
929  std::cout << "Volume Steel: " << m_VolumeSteel / cm3 << " cm^3" << std::endl;
930  std::cout << "Volume Scintillator: " << m_VolumeScintillator / cm3 << " cm^3" << std::endl;
931  std::cout << "Volume Air: " << (m_VolumeEnvelope - m_VolumeSteel - m_VolumeScintillator) / cm3 << " cm^3" << std::endl;
932  }
933  return;
934 }
935 
937 {
938  auto it = m_ScintiTilePhysVolMap.find(volume);
939  if (it != m_ScintiTilePhysVolMap.end())
940  {
941  return it->second;
942  }
943  std::cout << "could not locate volume " << volume->GetName()
944  << " in Inner Hcal scintillator map" << std::endl;
945  gSystem->Exit(1);
946  // that's dumb but code checkers do not know that gSystem->Exit()
947  // terminates, so using the standard exit() makes them happy
948  exit(1);
949 }
950 
951 // This is dulplicated code, we can get rid of it when we have the code to make towergeom for real data reco.
953 {
954  PHNodeIterator iter(topNode());
955  PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
956  if (!runNode)
957  {
958  std::cout << PHWHERE << "Run Node missing, exiting." << std::endl;
959  gSystem->Exit(1);
960  exit(1);
961  }
962  PHNodeIterator runIter(runNode);
963  PHCompositeNode *RunDetNode = dynamic_cast<PHCompositeNode *>(runIter.findFirst("PHCompositeNode", m_SuperDetector));
964  if (!RunDetNode)
965  {
966  RunDetNode = new PHCompositeNode(m_SuperDetector);
967  runNode->addNode(RunDetNode);
968  }
969  m_TowerGeomNodeName = "TOWERGEOM_" + m_SuperDetector;
970  m_RawTowerGeom = findNode::getClass<RawTowerGeomContainer>(topNode(), m_TowerGeomNodeName);
971  if (!m_RawTowerGeom)
972  {
975  RunDetNode->addNode(newNode);
976  }
979  m_RawTowerGeom->set_radius(innerrad);
980  m_RawTowerGeom->set_thickness(thickness);
983  double geom_ref_radius = innerrad + thickness / 2.;
984  double phistart = m_Params->get_double_param("phistart");
985  if (!std::isfinite(phistart))
986  {
987  std::cout << PHWHERE << " phistart is not finite: " << phistart
988  << ", exiting now (this will crash anyway)" << std::endl;
989  gSystem->Exit(1);
990  }
991  for (int i = 0; i < m_Params->get_int_param(PHG4HcalDefs::n_towers); i++)
992  {
993  double phiend = phistart + 2. * M_PI / m_Params->get_int_param(PHG4HcalDefs::n_towers);
994  std::pair<double, double> range = std::make_pair(phiend, phistart);
995  phistart = phiend;
996  m_RawTowerGeom->set_phibounds(i, range);
997  }
998  double etalowbound = -m_Params->get_double_param("scinti_eta_coverage_neg");
999  for (int i = 0; i < m_Params->get_int_param("etabins"); i++)
1000  {
1001  // double etahibound = etalowbound + 2.2 / get_int_param("etabins");
1002  double etahibound = etalowbound +
1003  (m_Params->get_double_param("scinti_eta_coverage_neg") + m_Params->get_double_param("scinti_eta_coverage_pos")) / m_Params->get_int_param("etabins");
1004  std::pair<double, double> range = std::make_pair(etalowbound, etahibound);
1005  m_RawTowerGeom->set_etabounds(i, range);
1006  etalowbound = etahibound;
1007  }
1008  for (int iphi = 0; iphi < m_RawTowerGeom->get_phibins(); iphi++)
1009  {
1010  for (int ieta = 0; ieta < m_RawTowerGeom->get_etabins(); ieta++)
1011  {
1013 
1014  const double x(geom_ref_radius * cos(m_RawTowerGeom->get_phicenter(iphi)));
1015  const double y(geom_ref_radius * sin(m_RawTowerGeom->get_phicenter(iphi)));
1016  const double z(geom_ref_radius / tan(PHG4Utils::get_theta(m_RawTowerGeom->get_etacenter(ieta))));
1017 
1019  if (tg)
1020  {
1021  if (Verbosity() > 0)
1022  {
1023  std::cout << "IHCalDetector::InitRun - Tower geometry " << key << " already exists" << std::endl;
1024  }
1025 
1026  if (fabs(tg->get_center_x() - x) > 1e-4)
1027  {
1028  std::cout << "IHCalDetector::InitRun - Fatal Error - duplicated Tower geometry " << key << " with existing x = " << tg->get_center_x() << " and expected x = " << x
1029  << std::endl;
1030 
1031  return;
1032  }
1033  if (fabs(tg->get_center_y() - y) > 1e-4)
1034  {
1035  std::cout << "IHCalDetector::InitRun - Fatal Error - duplicated Tower geometry " << key << " with existing y = " << tg->get_center_y() << " and expected y = " << y
1036  << std::endl;
1037  return;
1038  }
1039  if (fabs(tg->get_center_z() - z) > 1e-4)
1040  {
1041  std::cout << "IHCalDetector::InitRun - Fatal Error - duplicated Tower geometry " << key << " with existing z= " << tg->get_center_z() << " and expected z = " << z
1042  << std::endl;
1043  return;
1044  }
1045  }
1046  else
1047  {
1048  if (Verbosity() > 0)
1049  {
1050  std::cout << "IHCalDetector::InitRun - building tower geometry " << key << "" << std::endl;
1051  }
1052 
1053  tg = new RawTowerGeomv1(key);
1054 
1055  tg->set_center_x(x);
1056  tg->set_center_y(y);
1057  tg->set_center_z(z);
1059  }
1060  }
1061  }
1062  if (Verbosity() > 0)
1063  {
1065  }
1066 }