Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4HcalDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4HcalDetector.cc
1 #include "PHG4HcalDetector.h"
2 #include "PHG4CylinderGeomContainer.h"
3 #include "PHG4CylinderGeomv3.h"
4 
5 #include <g4main/PHG4Detector.h> // for PHG4Detector
6 #include <g4main/PHG4Utils.h>
7 
8 #include <phool/PHCompositeNode.h>
9 #include <phool/PHIODataNode.h>
10 #include <phool/PHNode.h> // for PHNode
11 #include <phool/PHNodeIterator.h> // for PHNodeIterator
12 #include <phool/PHObject.h> // for PHObject
13 #include <phool/getClass.h>
14 
15 #include <Geant4/G4Colour.hh> // for G4Colour
16 #include <Geant4/G4Cons.hh>
17 #include <Geant4/G4ExtrudedSolid.hh>
18 #include <Geant4/G4LogicalVolume.hh>
19 #include <Geant4/G4Material.hh>
20 #include <Geant4/G4PVPlacement.hh>
21 #include <Geant4/G4PhysicalConstants.hh>
22 #include <Geant4/G4RotationMatrix.hh> // for G4RotationMatrix
23 #include <Geant4/G4String.hh> // for G4String
24 #include <Geant4/G4SubtractionSolid.hh>
25 #include <Geant4/G4SystemOfUnits.hh> // for cm, deg
26 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
27 #include <Geant4/G4Transform3D.hh> // for G4Transform3D
28 #include <Geant4/G4Tubs.hh>
29 #include <Geant4/G4TwoVector.hh>
30 #include <Geant4/G4VisAttributes.hh>
31 
32 #include <algorithm> // for max
33 #include <cmath> // for sin, cos, sqrt, M_PI, asin
34 #include <cstdlib> // for exit
35 #include <iostream> // for operator<<, basic_ostream
36 #include <sstream>
37 #include <utility> // for pair
38 #include <vector> // for vector
39 
40 class PHG4CylinderGeom;
41 
42 // uncomment if you want to make a graphics display where the slats are visible
43 // it makes them stick out of the hcal for visibility
44 // NEVER EVER RUN REAL SIMS WITH THIS
45 //#define DISPLAY
46 
47 int PHG4HcalDetector::INACTIVE = -100;
48 //_______________________________________________________________
49 //note this inactive thickness is ~1.5% of a radiation length
50 PHG4HcalDetector::PHG4HcalDetector(PHG4Subsystem* subsys, PHCompositeNode* Node, const std::string& dnam, const int lyr)
51  : PHG4Detector(subsys, Node, dnam)
52  , TrackerMaterial(nullptr)
53  , TrackerThickness(100 * cm)
54  , cylinder_logic(nullptr)
55  , cylinder_physi(nullptr)
56  , radius(100 * cm)
57  , length(100 * cm)
58  , xpos(0 * cm)
59  , ypos(0 * cm)
60  , zpos(0 * cm)
61  , _sciTilt(0)
62  , _sciWidth(0.6 * cm)
63  , _sciNum(100)
64  , _sciPhi0(0)
65  , _region(nullptr)
66  , active(0)
67  , absorberactive(0)
68  , layer(lyr)
69 {
70 }
71 
72 //_______________________________________________________________
73 //_______________________________________________________________
74 int PHG4HcalDetector::IsInCylinderActive(const G4VPhysicalVolume* volume)
75 {
76  // std::cout << "checking detector" << std::endl;
77  if (active && box_vol.find(volume) != box_vol.end())
78  {
79  return box_vol.find(volume)->second;
80  }
81  if (absorberactive && volume == cylinder_physi)
82  {
83  return -1;
84  }
85  return INACTIVE;
86 }
87 
88 //_______________________________________________________________
89 void PHG4HcalDetector::ConstructMe(G4LogicalVolume* logicWorld)
90 {
91  TrackerMaterial = GetDetectorMaterial(material);
92 
93  G4Tubs* _cylinder_solid = new G4Tubs(G4String(GetName()),
94  radius,
95  radius + TrackerThickness,
96  length / 2.0, 0, twopi);
97  double innerlength = PHG4Utils::GetLengthForRapidityCoverage(radius) * 2;
98  double deltalen = (length - innerlength) / 2.; // length difference on one side
99  double cone_size_multiplier = 1.01; // 1 % larger
100  double cone_thickness = TrackerThickness * cone_size_multiplier;
101  double inner_cone_radius = radius - ((cone_thickness - TrackerThickness) / 2.);
102  double cone_length = deltalen * cone_size_multiplier;
103  G4Cons* cone2 = new G4Cons("conehead2",
104  inner_cone_radius, inner_cone_radius,
105  inner_cone_radius, inner_cone_radius + cone_thickness,
106  cone_length / 2.0, 0, twopi);
107  G4Cons* cone1 = new G4Cons("conehead",
108  inner_cone_radius, inner_cone_radius + cone_thickness,
109  inner_cone_radius, inner_cone_radius,
110  cone_length / 2.0, 0, twopi);
111  double delta_len = cone_length - deltalen;
112  G4ThreeVector zTransneg(0, 0, -(length - cone_length + delta_len) / 2.0);
113  G4ThreeVector zTranspos(0, 0, (length - cone_length + delta_len) / 2.0);
114  G4SubtractionSolid* subtraction_tmp =
115  new G4SubtractionSolid("Cylinder-Cone", _cylinder_solid, cone1, nullptr, zTransneg);
116  G4SubtractionSolid* subtraction =
117  new G4SubtractionSolid("Cylinder-Cone-Cone", subtraction_tmp, cone2, nullptr, zTranspos);
118  cylinder_logic = new G4LogicalVolume(subtraction,
119  TrackerMaterial,
120  G4String(GetName()),
121  nullptr, nullptr, nullptr);
122  G4VisAttributes* VisAtt = new G4VisAttributes();
123  VisAtt->SetColour(G4Colour::Grey());
124  VisAtt->SetVisibility(true);
125  VisAtt->SetForceSolid(true);
126  cylinder_logic->SetVisAttributes(VisAtt);
127 
128  cylinder_physi = new G4PVPlacement(nullptr, G4ThreeVector(xpos, ypos, zpos),
129  cylinder_logic,
130  G4String(GetName()),
131  logicWorld, false, false, OverlapCheck());
132  // Figure out corners of scintillator inside the containing G4Tubs.
133  // Work our way around the scintillator cross section in a counter
134  // clockwise fashion: ABCD
135  double r1 = radius;
136  double r2 = radius + TrackerThickness;
137 
138  // The coordinates of the inner corners of the scintillator
139  double x4 = r1;
140  double y4 = _sciWidth / 2.0;
141 
142  double x1 = r1;
143  double y1 = -y4;
144 
145  double a = _sciTilt * M_PI / 180.0;
146 
147  // The parametric equation for the line from A along the side of the
148  // scintillator is (x,y) = (x1,y1) + u * (cos(a), sin(a))
149  double A = 1.0;
150  double B = 2 * (x1 * cos(a) + y1 * sin(a));
151  double C = x1 * x1 + y1 * y1 - r2 * r2;
152  double D = B * B - 4 * A * C;
153 
154  // The only sensible solution, given our definitions, is u > 0.
155  double u = (-B + sqrt(D)) / 2 * A;
156 
157  // Now we can determine one of the outer corners
158  double x2 = x1 + u * cos(a);
159  double y2 = y1 + u * sin(a);
160 
161  // Similar procedure for (x3,y3) as for (x2,y2)
162  A = 1.0;
163  B = 2 * (x4 * cos(a) + y4 * sin(a));
164  C = x4 * x4 + y4 * y4 - r2 * r2;
165  D = B * B - 4 * A * C;
166  u = (-B + sqrt(D)) / 2 * A;
167 
168  double x3 = x4 + u * cos(a);
169  double y3 = y4 + u * sin(a);
170 
171  // Now we've got a four-sided "z-section".
172  G4TwoVector v1(x1, y1);
173  G4TwoVector v2(x2, y2);
174  G4TwoVector v3(x3, y3);
175  G4TwoVector v4(x4, y4);
176 
177  std::vector<G4TwoVector> vertexes;
178  vertexes.push_back(v1);
179  vertexes.push_back(v2);
180  vertexes.push_back(v3);
181  vertexes.push_back(v4);
182 
183  G4TwoVector zero(0, 0);
184  // if you want to make displays where the structure of the hcal is visible
185  // add 20 cm to the length of the scintillators
186 #ifdef DISPLAY
187  double blength = length + 20;
188 #else
189  double blength = length;
190 #endif
191  G4ExtrudedSolid* _box_solid = new G4ExtrudedSolid("_BOX",
192  vertexes,
193  blength / 2.0,
194  zero, 1.0,
195  zero, 1.0);
196 
197  // double boxlen_half = GetLength(_sciTilt * M_PI / 180.);
198  // G4Box* _box_solid = new G4Box("_BOX", boxlen_half, _sciWidth / 2.0, length / 2.0);
199 
200  G4Material* boxmat = GetDetectorMaterial("G4_POLYSTYRENE");
201  G4SubtractionSolid* subtractionbox_tmp =
202  new G4SubtractionSolid("Box-Cone", _box_solid, cone1, nullptr, zTransneg);
203  G4SubtractionSolid* subtractionbox =
204  new G4SubtractionSolid("Box-Cone-Cone", subtractionbox_tmp, cone2, nullptr, zTranspos);
205  G4LogicalVolume* box_logic = new G4LogicalVolume(subtractionbox,
206  boxmat, G4String("BOX"),
207  nullptr, nullptr, nullptr);
208  VisAtt = new G4VisAttributes();
209  PHG4Utils::SetColour(VisAtt, "G4_POLYSTYRENE");
210  VisAtt->SetVisibility(true);
211  VisAtt->SetForceSolid(true);
212  box_logic->SetVisAttributes(VisAtt);
213 
214  double phi_increment = 360. / _sciNum;
215  std::ostringstream slatname;
216  for (int i = 0; i < _sciNum; i++)
217  {
218  double phi = (i + _sciPhi0) * phi_increment;
219  G4ThreeVector myTrans = G4ThreeVector(0, 0, 0);
220  G4RotationMatrix Rot(0, 0, 0);
221  Rot.rotateZ(phi * deg);
222  slatname.str("");
223  slatname << "SLAT_" << i;
224  G4VPhysicalVolume* box_vol_tmp = new G4PVPlacement(G4Transform3D(Rot, G4ThreeVector(myTrans)),
225  box_logic,
226  G4String(slatname.str()),
227  cylinder_logic, false, false, OverlapCheck());
228  box_vol[box_vol_tmp] = i;
229  }
230  if (active)
231  {
232  std::ostringstream geonode;
233  if (superdetector != "NONE")
234  {
235  geonode << "CYLINDERGEOM_" << superdetector;
236  }
237  else
238  {
239  geonode << "CYLINDERGEOM_" << detector_type << "_" << layer;
240  }
241  PHG4CylinderGeomContainer* geo = findNode::getClass<PHG4CylinderGeomContainer>(topNode(), geonode.str());
242  if (!geo)
243  {
244  geo = new PHG4CylinderGeomContainer();
245  PHNodeIterator iter(topNode());
246  PHCompositeNode* runNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "RUN"));
247  PHIODataNode<PHObject>* newNode = new PHIODataNode<PHObject>(geo, geonode.str(), "PHObject");
248  runNode->addNode(newNode);
249  }
250  // here in the detector class we have internal units, convert to cm
251  // before putting into the geom object
252  PHG4CylinderGeom* mygeom = new PHG4CylinderGeomv3(radius / cm, (zpos - length / 2.) / cm, (zpos + length / 2.) / cm, TrackerThickness / cm, _sciNum, _sciTilt * M_PI / 180.0, _sciPhi0 * M_PI / 180.0);
253  geo->AddLayerGeom(layer, mygeom);
254  // geo->identify();
255  }
256 }
257 
258 double
259 PHG4HcalDetector::GetLength(const double phi) const
260 {
261  double c = radius + TrackerThickness / 2.;
262  double b = radius;
263  double singamma = sin(phi) * c / b;
264  double gamma = M_PI - asin(singamma);
265  double alpha = M_PI - gamma - phi;
266  double a = c * sin(alpha) / singamma;
267  return a;
268 }
269 
270 void PHG4HcalDetector::Print(const std::string& /*what*/) const
271 {
272  std::cout << "radius: " << radius << std::endl;
273  return;
274 }