Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4HcalSubsystem.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4HcalSubsystem.cc
1 #include "PHG4HcalSubsystem.h"
2 
3 #include "PHG4HcalDetector.h"
5 
7 #include <g4main/PHG4Utils.h>
8 
10 #include <phool/PHIODataNode.h> // for PHIODataNode
11 #include <phool/PHNode.h> // for PHNode
12 #include <phool/PHNodeIterator.h> // for PHNodeIterator
13 #include <phool/PHObject.h> // for PHObject
14 #include <phool/getClass.h>
15 
16 #include <Geant4/G4String.hh> // for G4String
17 #include <Geant4/G4SystemOfUnits.hh> // for cm
18 #include <Geant4/G4Types.hh> // for G4double
19 
20 #include <cmath> // for asin, cos, sin, sqrt, M_PI
21 #include <cstdlib> // for NULL, exit
22 #include <iostream> // for operator<<, basic_ostream
23 #include <sstream>
24 
25 class PHG4Detector;
26 class PHG4SteppingAction;
27 
28 //_______________________________________________________________________
30  : detector_(nullptr)
31  , steppingAction_(nullptr)
32  , radius(100)
33  , length(100)
34  , xpos(0)
35  , ypos(0)
36  , zpos(0)
37  , lengthViaRapidityCoverage(true)
38  , TrackerThickness(100)
39  , material("G4_Fe")
40  , _sciTilt(0)
41  , _sciWidth(0.7)
42  , _sciNum(256)
43  , _sciPhi0(0)
44  , active(0)
45  , absorberactive(0)
46  , layer(lyr)
47  , detector_type(na)
48  , superdetector("NONE")
49  , light_scint_model_(true)
50  , light_balance_(false)
51  , light_balance_inner_radius_(0.0 * cm)
52  , light_balance_inner_corr_(1.0)
53  , light_balance_outer_radius_(10.0 * cm)
54  , light_balance_outer_corr_(1.0)
55 {
56  // put the layer into the name so we get unique names
57  // for multiple SVX layers
58  // std::ostringstream nam;
59  // nam << na << "_" << lyr;
60  Name(na + "_" + std::to_string(lyr));
61 }
62 
63 //_______________________________________________________________________
65 {
66  // create hit list only for active layers
67  PHNodeIterator iter(topNode);
68  PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
69  // create detector
70  detector_ = new PHG4HcalDetector(this, topNode, Name(), layer);
72  G4double detlength = length;
74  {
76  }
77  detector_->SetLength(detlength);
89  if (active)
90  {
91  std::string nodename;
92  if (superdetector != "NONE")
93  {
94  nodename = "G4HIT_" + superdetector;
95  }
96  else
97  {
98  nodename = "G4HIT_" + detector_type + "_" + std::to_string(layer);
99  }
100  PHG4HitContainer* cylinder_hits = findNode::getClass<PHG4HitContainer>(topNode, nodename);
101  if (!cylinder_hits)
102  {
103  dstNode->addNode(new PHIODataNode<PHObject>(cylinder_hits = new PHG4HitContainer(nodename), nodename, "PHObject"));
104  }
105  cylinder_hits->AddLayer(layer);
106  if (absorberactive)
107  {
108  if (superdetector != "NONE")
109  {
110  nodename = "G4HIT_ABSORBER_" + superdetector;
111  }
112  else
113  {
114  nodename = "G4HIT_ABSORBER_" + detector_type + "_" + std::to_string(layer);
115  }
116  PHG4HitContainer* cylinder_hits_2 = findNode::getClass<PHG4HitContainer>(topNode, nodename);
117  if (!cylinder_hits_2)
118  {
119  dstNode->addNode(new PHIODataNode<PHObject>(cylinder_hits_2 = new PHG4HitContainer(nodename), nodename, "PHObject"));
120  }
121  cylinder_hits_2->AddLayer(layer);
122  }
124  steppingAction_->set_zmin(zpos - detlength / 2.);
125  steppingAction_->set_zmax(zpos + detlength / 2.);
126  if (light_balance_)
127  {
133  }
134  }
135 
136  return 0;
137 }
138 
139 //_______________________________________________________________________
141 {
142  // pass top node to stepping action so that it gets
143  // relevant nodes needed internally
144  if (steppingAction_)
145  {
147  }
148  return 0;
149 }
150 
151 //_______________________________________________________________________
153 {
154  return detector_;
155 }
156 
157 //_______________________________________________________________________
159 {
160  return steppingAction_;
161 }
162 
163 void PHG4HcalSubsystem::Print(const std::string& what) const
164 {
165  detector_->Print(what);
166  return;
167 }
168 
170 {
171  if (ncross == 0)
172  {
173  std::cout << Name() << " invalid crossing number " << ncross
174  << " how do you think we can construct a meaningful detector with this number????" << std::endl;
175  exit(1);
176  }
177  // The delta phi angle between 2 adjacent slats is just 360/nslats. If
178  // there is just one crossing the tips of each slat can extend from
179  // phi-delta-phi/2 to phi+delta-phi/2. G4 rotates around the center of a
180  // slat so we are dealing in increments of just delta-phi/2. This
181  // has to be multiplied by the number of crossings we want.
182  //
183  // To find this tilt angle we have to calculate a triangle
184  // the long sides are the outer radius and the
185  // inner radius + 1/2 tracker thickness
186  // the angle between these sides is just (360/nslat)/2*(cross)
187  // the we use a/sin(alpha) = b/sin(beta)
188  // and c^2 = a^2+b^2 -2ab*cos(gamma)
189  // the solution is not unique, we have to pick 180-beta
190  // but the slat angle is 180-beta (it's on the other side of the triangle
191  // just draw the damned thing if this is confusing)
192  double sign = 1;
193  if (ncross < 0)
194  {
195  sign = -1;
196  }
197  int ncr = std::abs(ncross);
198  double thick = TrackerThickness;
199  double c = radius + thick / 2.;
200  double b = radius + thick;
201 
202  double alpha = 0;
203  alpha = ((360. / _sciNum * M_PI / 180.) / 2.) * ncr;
204  double sinb = sin(alpha) * b / (sqrt(b * b + c * c - 2 * b * c * cos(alpha)));
205  double beta = asin(sinb) * 180. / M_PI; // this is already the slat angle
206  _sciTilt = beta * sign;
207  // print out triangle
208  // double tbeta = 180 - beta;
209  // double gamma = 180 - (alpha * 180. / M_PI) - tbeta;
210  // double a = b * sin(alpha) / sinb;
211  // std::cout << "triangle length: a: " << a
212  // << ", b: " << b << ", c: " << c << std::endl;
213  // std::cout << "triangle angle: alpha: " << alpha*180./M_PI
214  // << ", beta: " << tbeta
215  // << ", gamma: " << gamma
216  // << std::endl;
217  std::cout << Name() << ": SetTiltViaNcross(" << ncross << ") setting slat angle to: " << _sciTilt << " degrees" << std::endl;
218  return;
219 }