ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4EICForwardEcalDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4EICForwardEcalDetector.cc
2 
3 #include "PHG4ForwardEcalDetector.h" // for PHG4ForwardEcalDetector
5 
6 #include <phparameter/PHParameters.h>
7 
8 #include <Geant4/G4Box.hh>
9 #include <Geant4/G4Cons.hh>
10 #include <Geant4/G4LogicalVolume.hh>
11 #include <Geant4/G4Material.hh>
12 #include <Geant4/G4PVPlacement.hh>
13 #include <Geant4/G4RotationMatrix.hh> // for G4RotationMatrix
14 #include <Geant4/G4String.hh> // for G4String
15 #include <Geant4/G4SystemOfUnits.hh> // for cm, mm
16 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
17 #include <Geant4/G4Transform3D.hh> // for G4Transform3D
18 #include <Geant4/G4Types.hh> // for G4double, G4int
19 
20 #include <TSystem.h>
21 
22 #include <cmath> // for M_PI
23 #include <fstream>
24 #include <iostream>
25 #include <sstream>
26 #include <utility> // for pair, make_pair
27 
28 class G4VSolid;
29 class PHCompositeNode;
30 
31 using namespace std;
32 
33 //_______________________________________________________________________
35  : PHG4ForwardEcalDetector(subsys, Node, parameters, dnam)
36  , _tower_dx(30 * mm)
37  , _tower_dy(30 * mm)
38  , _tower_dz(170.0 * mm)
39  , _materialScintillator("G4_POLYSTYRENE")
40  , _materialAbsorber("G4_Pb")
41 {
42 }
43 
44 //_______________________________________________________________________
46 {
47  if (Verbosity() > 0)
48  {
49  cout << "PHG4EICForwardEcalDetector: Begin Construction" << endl;
50  }
51 
52  /* Read parameters for detector construction and mappign from file */
54 
55  /* Create the cone envelope = 'world volume' for the crystal calorimeter */
56  G4Material* Air = GetDetectorMaterial("G4_AIR");
57 
58  G4VSolid* ecal_envelope_solid = new G4Cons("hEcal_envelope_solid",
59  GetRMin(0), GetRMax(0),
60  GetRMin(1), GetRMax(1),
61  GetdZ() / 2.,
62  0, 2 * M_PI);
63 
64  G4LogicalVolume* ecal_envelope_log = new G4LogicalVolume(ecal_envelope_solid, Air, G4String("hEcal_envelope"), 0, 0, 0);
65 
66  /* Define visualization attributes for envelope cone */
67  GetDisplayAction()->AddVolume(ecal_envelope_log, "Envelope");
68 
69  /* Define rotation attributes for envelope cone */
70  G4RotationMatrix ecal_rotm;
71  ecal_rotm.rotateX(GetXRot());
72  ecal_rotm.rotateY(GetYRot());
73  ecal_rotm.rotateZ(GetZRot());
74 
75  /* Place envelope cone in simulation */
76  ostringstream name_envelope;
77  name_envelope.str("");
78  name_envelope << TowerLogicNamePrefix() << "_envelope";
79 
81  ecal_envelope_log, name_envelope.str().c_str(), logicWorld, 0, false, OverlapCheck());
82 
83  /* Construct single calorimeter tower */
84  G4LogicalVolume* singletower = ConstructTower();
85 
86  /* Place calorimeter tower within envelope */
87  PlaceTower(ecal_envelope_log, singletower);
88 
89  return;
90 }
91 
92 //_______________________________________________________________________
95 {
96  if (Verbosity() > 0)
97  {
98  cout << "PHG4EICForwardEcalDetector: Build logical volume for single tower..." << endl;
99  }
100 
101  /* create logical volume for single tower */
102  G4Material* material_air = GetDetectorMaterial("G4_AIR");
103 
104  G4VSolid* single_tower_solid = new G4Box(G4String("single_tower_solid"),
105  _tower_dx / 2.0,
106  _tower_dy / 2.0,
107  _tower_dz / 2.0);
108 
109  G4LogicalVolume* single_tower_logic = new G4LogicalVolume(single_tower_solid,
110  material_air,
111  "single_tower_logic",
112  0, 0, 0);
113 
114  /* create geometry volumes for scintillator and absorber plates to place inside single_tower */
115  G4int nlayers = 60;
116  G4double thickness_layer = _tower_dz / (float) nlayers;
117  G4double thickness_absorber = thickness_layer / 3.0 * 2.0; // 2/3rd absorber
118  G4double thickness_scintillator = thickness_layer / 3.0 * 1.0; // 1/3rd scintillator
119  // PHENIX EMCal JGL 12/27/2015
120  // G4int nlayers = 66;
121  // G4double thickness_layer = _tower_dz/(float)nlayers;
122  // G4double thickness_absorber = thickness_layer*0.23; // 27% absorber by length
123  // G4double thickness_scintillator = thickness_layer*0.73; // 73% scintillator by length
124 
125  G4VSolid* solid_absorber = new G4Box(G4String("single_plate_absorber_solid"),
126  _tower_dx / 2.0,
127  _tower_dy / 2.0,
128  thickness_absorber / 2.0);
129 
130  G4VSolid* solid_scintillator = new G4Box(G4String("single_plate_scintillator"),
131  _tower_dx / 2.0,
132  _tower_dy / 2.0,
133  thickness_scintillator / 2.0);
134 
135  /* create logical volumes for scintillator and absorber plates to place inside single_tower */
136  G4Material* material_scintillator = GetDetectorMaterial(_materialScintillator.c_str());
137  G4Material* material_absorber = GetDetectorMaterial(_materialAbsorber.c_str());
138 
139  G4LogicalVolume* logic_absorber = new G4LogicalVolume(solid_absorber,
140  material_absorber,
141  "single_plate_absorber_logic",
142  0, 0, 0);
143 
144  G4LogicalVolume* logic_scint = new G4LogicalVolume(solid_scintillator,
145  material_scintillator,
146  "hEcal_scintillator_plate_logic",
147  0, 0, 0);
148  AbsorberLogicalVolSetInsert(logic_absorber);
149  ScintiLogicalVolSetInsert(logic_scint);
150  GetDisplayAction()->AddVolume(logic_absorber, "Absorber");
151  GetDisplayAction()->AddVolume(logic_scint, "Scintillator");
152 
153  /* place physical volumes for absorber and scintillator plates */
154  G4double xpos_i = 0;
155  G4double ypos_i = 0;
156  G4double zpos_i = (-1 * _tower_dz / 2.0) + thickness_absorber / 2.0;
157 
158  ostringstream name_absorber;
159  name_absorber.str("");
160  name_absorber << TowerLogicNamePrefix() << "_single_plate_absorber";
161 
162  ostringstream name_scintillator;
163  name_scintillator.str("");
164  name_scintillator << TowerLogicNamePrefix() << "_single_plate_scintillator";
165 
166  for (int i = 1; i <= nlayers; i++)
167  {
168  new G4PVPlacement(0, G4ThreeVector(xpos_i, ypos_i, zpos_i),
169  logic_absorber,
170  name_absorber.str().c_str(),
171  single_tower_logic,
172  0, 0, OverlapCheck());
173 
174  zpos_i += (thickness_absorber / 2. + thickness_scintillator / 2.);
175 
176  new G4PVPlacement(0, G4ThreeVector(xpos_i, ypos_i, zpos_i),
177  logic_scint,
178  name_scintillator.str().c_str(),
179  single_tower_logic,
180  0, 0, OverlapCheck());
181 
182  zpos_i += (thickness_absorber / 2. + thickness_scintillator / 2.);
183  }
184 
185  GetDisplayAction()->AddVolume(single_tower_logic, "ScintillatorSingleTower");
186 
187  if (Verbosity() > 0)
188  {
189  cout << "PHG4EICForwardEcalDetector: Building logical volume for single tower done." << endl;
190  }
191 
192  return single_tower_logic;
193 }
194 
196 {
197  /* Loop over all tower positions in vector and place tower */
198  typedef std::map<std::string, towerposition>::iterator it_type;
199 
200  for (it_type iterator = _map_tower.begin(); iterator != _map_tower.end(); ++iterator)
201  {
202  if (Verbosity() > 0)
203  {
204  cout << "PHG4EICForwardEcalDetector: Place tower " << iterator->first
205  << " at x = " << iterator->second.x << " , y = " << iterator->second.y << " , z = " << iterator->second.z << endl;
206  }
207 
208  new G4PVPlacement(0, G4ThreeVector(iterator->second.x, iterator->second.y, iterator->second.z),
209  singletower,
210  iterator->first.c_str(),
211  ecalenvelope,
212  0, 0, OverlapCheck());
213  }
214 
215  return 0;
216 }
217 
219 {
220  /* Open the datafile, if it won't open return an error */
221  ifstream istream_mapping;
222  istream_mapping.open(GetParams()->get_string_param("mapping_file"));
223  if (!istream_mapping.is_open())
224  {
225  cout << "ERROR in PHG4EICForwardEcalDetector: Failed to open mapping file " << GetParams()->get_string_param("mapping_file") << endl;
226  gSystem->Exit(1);
227  }
228 
229  /* loop over lines in file */
230  string line_mapping;
231  while (getline(istream_mapping, line_mapping))
232  {
233  /* Skip lines starting with / including a '#' */
234  if (line_mapping.find("#") != string::npos)
235  {
236  if (Verbosity() > 0)
237  {
238  cout << "PHG4EICForwardEcalDetector: SKIPPING line in mapping file: " << line_mapping << endl;
239  }
240  continue;
241  }
242 
243  istringstream iss(line_mapping);
244 
245  /* If line starts with keyword Tower, add to tower positions */
246  if (line_mapping.find("Tower ") != string::npos)
247  {
248  unsigned idx_j, idx_k, idx_l;
249  G4double pos_x, pos_y, pos_z;
250  G4double size_x, size_y, size_z;
251  G4double rot_x, rot_y, rot_z;
252  G4double dummy;
253  string dummys;
254 
255  /* read string- break if error */
256  if (!(iss >> dummys >> dummy >> idx_j >> idx_k >> idx_l >> pos_x >> pos_y >> pos_z >> size_x >> size_y >> size_z >> rot_x >> rot_y >> rot_z))
257  {
258  cerr << "ERROR in PHG4EICForwardEcalDetector: Failed to read line in mapping file " << GetParams()->get_string_param("mapping_file") << endl;
259  gSystem->Exit(1);
260  }
261 
262  /* Construct unique name for tower */
263  /* Mapping file uses cm, this class uses mm for length */
264  ostringstream towername;
265  towername.str("");
266  towername << TowerLogicNamePrefix() << "_j_" << idx_j << "_k_" << idx_k;
267 
268  /* Add Geant4 units */
269  pos_x = pos_x * cm;
270  pos_y = pos_y * cm;
271  pos_z = pos_z * cm;
272 
273  /* insert tower into tower map */
274  towerposition tower_new;
275  tower_new.x = pos_x;
276  tower_new.y = pos_y;
277  tower_new.z = pos_z;
278  _map_tower.insert(make_pair(towername.str(), tower_new));
279  }
280  else
281  {
282  /* If this line is not a comment and not a tower, save parameter as string / value. */
283  string parname;
284  G4double parval;
285 
286  /* read string- break if error */
287  if (!(iss >> parname >> parval))
288  {
289  cerr << "ERROR in PHG4EICForwardEcalDetector: Failed to read line in mapping file " << GetParams()->get_string_param("mapping_file") << endl;
290  gSystem->Exit(1);
291  }
292  InsertParam(parname, parval);
293  }
294  }
295 
296  /* Update member variables for global parameters based on parsed parameter file */
297  std::map<string, double>::const_iterator parit;
298 
299  parit = FindIter("Gtower_dx");
300  if (parit != EndIter())
301  _tower_dx = parit->second * cm;
302 
303  parit = FindIter("Gtower_dy");
304  if (parit != EndIter())
305  _tower_dy = parit->second * cm;
306 
307  parit = FindIter("Gtower_dz");
308  if (parit != EndIter())
309  _tower_dz = parit->second * cm;
310 
311  ostringstream rad;
312  for (int i = 0; i < 2; i++)
313  {
314  int index = i + 1;
315  rad.str("");
316  rad << "Gr" << index << "_inner";
317  parit = FindIter(rad.str());
318  if (parit != EndIter())
319  {
320  SetRMin(i, parit->second * cm);
321  }
322  rad.str("");
323  rad << "Gr" << index << "_outer";
324  parit = FindIter(rad.str());
325  if (parit != EndIter())
326  {
327  SetRMax(i, parit->second * cm);
328  }
329  }
330  parit = FindIter("Gdz");
331  if (parit != EndIter())
332  {
333  SetdZ(parit->second * cm);
334  }
335  parit = FindIter("Gx0");
336  if (parit != EndIter())
337  {
338  SetPlaceX(parit->second * cm);
339  }
340 
341  parit = FindIter("Gy0");
342  if (parit != EndIter())
343  {
344  SetPlaceY(parit->second * cm);
345  }
346  parit = FindIter("Gz0");
347  if (parit != EndIter())
348  {
349  SetPlaceZ(parit->second * cm);
350  }
351  parit = FindIter("Grot_x");
352  if (parit != EndIter())
353  SetXRot(parit->second);
354 
355  parit = FindIter("Grot_y");
356  if (parit != EndIter())
357  SetYRot(parit->second);
358 
359  parit = FindIter("Grot_z");
360  if (parit != EndIter())
361  SetZRot(parit->second);
362 
363  return 0;
364 }