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