ECCE @ EIC 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 using namespace std;
29 
30 //_______________________________________________________________________
31 PHG4HcalSubsystem::PHG4HcalSubsystem(const std::string& na, const int lyr)
32  : detector_(nullptr)
33  , steppingAction_(nullptr)
34  , radius(100)
35  , length(100)
36  , xpos(0)
37  , ypos(0)
38  , zpos(0)
39  , lengthViaRapidityCoverage(true)
40  , TrackerThickness(100)
41  , material("G4_Fe")
42  , _sciTilt(0)
43  , _sciWidth(0.7)
44  , _sciNum(256)
45  , _sciPhi0(0)
46  , active(0)
47  , absorberactive(0)
48  , layer(lyr)
49  , detector_type(na)
50  , superdetector("NONE")
51  , light_scint_model_(true)
52  , light_balance_(false)
53  , light_balance_inner_radius_(0.0 * cm)
54  , light_balance_inner_corr_(1.0)
55  , light_balance_outer_radius_(10.0 * cm)
56  , light_balance_outer_corr_(1.0)
57 {
58  // put the layer into the name so we get unique names
59  // for multiple SVX layers
60  ostringstream nam;
61  nam << na << "_" << lyr;
62  Name(nam.str().c_str());
63 }
64 
65 //_______________________________________________________________________
67 {
68  // create hit list only for active layers
69  PHNodeIterator iter(topNode);
70  PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
71  // create detector
72  detector_ = new PHG4HcalDetector(this, topNode, Name(), layer);
74  G4double detlength = length;
76  {
78  }
79  detector_->SetLength(detlength);
91  if (active)
92  {
93  ostringstream nodename;
94  if (superdetector != "NONE")
95  {
96  nodename << "G4HIT_" << superdetector;
97  }
98  else
99  {
100  nodename << "G4HIT_" << detector_type << "_" << layer;
101  }
102  PHG4HitContainer* cylinder_hits = findNode::getClass<PHG4HitContainer>(topNode, nodename.str());
103  if (!cylinder_hits)
104  {
105  dstNode->addNode(new PHIODataNode<PHObject>(cylinder_hits = new PHG4HitContainer(nodename.str()), nodename.str(), "PHObject"));
106  }
107  cylinder_hits->AddLayer(layer);
108  if (absorberactive)
109  {
110  nodename.str("");
111  if (superdetector != "NONE")
112  {
113  nodename << "G4HIT_ABSORBER_" << superdetector;
114  }
115  else
116  {
117  nodename << "G4HIT_ABSORBER_" << detector_type << "_" << layer;
118  }
119  PHG4HitContainer* cylinder_hits_2 = findNode::getClass<PHG4HitContainer>(topNode, nodename.str());
120  if (!cylinder_hits_2)
121  {
122  dstNode->addNode(new PHIODataNode<PHObject>(cylinder_hits_2 = new PHG4HitContainer(nodename.str()), nodename.str(), "PHObject"));
123  }
124  cylinder_hits_2->AddLayer(layer);
125  }
127  steppingAction_->set_zmin(zpos - detlength / 2.);
128  steppingAction_->set_zmax(zpos + detlength / 2.);
129  if (light_balance_)
130  {
136  }
137  }
138 
139  return 0;
140 }
141 
142 //_______________________________________________________________________
144 {
145  // pass top node to stepping action so that it gets
146  // relevant nodes needed internally
147  if (steppingAction_)
148  {
150  }
151  return 0;
152 }
153 
154 //_______________________________________________________________________
156 {
157  return detector_;
158 }
159 
160 //_______________________________________________________________________
162 {
163  return steppingAction_;
164 }
165 
166 void PHG4HcalSubsystem::Print(const std::string& what) const
167 {
168  detector_->Print(what);
169  return;
170 }
171 
173 {
174  if (ncross == 0)
175  {
176  cout << Name() << " invalid crossing number " << ncross
177  << " how do you think we can construct a meaningful detector with this number????" << endl;
178  exit(1);
179  }
180  // The delta phi angle between 2 adjacent slats is just 360/nslats. If
181  // there is just one crossing the tips of each slat can extend from
182  // phi-delta-phi/2 to phi+delta-phi/2. G4 rotates around the center of a
183  // slat so we are dealing in increments of just delta-phi/2. This
184  // has to be multiplied by the number of crossings we want.
185  //
186  // To find this tilt angle we have to calculate a triangle
187  // the long sides are the outer radius and the
188  // inner radius + 1/2 tracker thickness
189  // the angle between these sides is just (360/nslat)/2*(cross)
190  // the we use a/sin(alpha) = b/sin(beta)
191  // and c^2 = a^2+b^2 -2ab*cos(gamma)
192  // the solution is not unique, we have to pick 180-beta
193  // but the slat angle is 180-beta (it's on the other side of the triangle
194  // just draw the damned thing if this is confusing)
195  double sign = 1;
196  if (ncross < 0)
197  {
198  sign = -1;
199  }
200  int ncr = fabs(ncross);
201  double thick = TrackerThickness;
202  double c = radius + thick / 2.;
203  double b = radius + thick;
204 
205  double alpha = 0;
206  alpha = ((360. / _sciNum * M_PI / 180.) / 2.) * ncr;
207  double sinb = sin(alpha) * b / (sqrt(b * b + c * c - 2 * b * c * cos(alpha)));
208  double beta = asin(sinb) * 180. / M_PI; // this is already the slat angle
209  _sciTilt = beta * sign;
210  // print out triangle
211  // double tbeta = 180 - beta;
212  // double gamma = 180 - (alpha * 180. / M_PI) - tbeta;
213  // double a = b * sin(alpha) / sinb;
214  // cout << "triangle length: a: " << a
215  // << ", b: " << b << ", c: " << c << endl;
216  // cout << "triangle angle: alpha: " << alpha*180./M_PI
217  // << ", beta: " << tbeta
218  // << ", gamma: " << gamma
219  // << endl;
220  cout << Name() << ": SetTiltViaNcross(" << ncross << ") setting slat angle to: " << _sciTilt << " degrees" << endl;
221  return;
222 }