ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4CrystalCalorimeterDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4CrystalCalorimeterDetector.cc
3 
4 #include <phparameter/PHParameters.h>
5 
6 #include <g4main/PHG4Detector.h> // for PHG4Detector
7 #include <g4main/PHG4DisplayAction.h> // for PHG4DisplayAction
8 #include <g4main/PHG4Subsystem.h>
9 
10 #include <phool/recoConsts.h>
11 
12 #include <Geant4/G4Box.hh>
13 #include <Geant4/G4Cons.hh>
14 #include <Geant4/G4Element.hh> // for G4Element
15 #include <Geant4/G4LogicalVolume.hh>
16 #include <Geant4/G4Material.hh>
17 #include <Geant4/G4PVPlacement.hh>
18 #include <Geant4/G4RotationMatrix.hh> // for G4RotationMatrix
19 #include <Geant4/G4String.hh> // for G4String
20 #include <Geant4/G4SubtractionSolid.hh>
21 #include <Geant4/G4SystemOfUnits.hh>
22 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
23 #include <Geant4/G4Transform3D.hh> // for G4Transform3D
24 #include <Geant4/G4Types.hh> // for G4double
25 
26 #include <TSystem.h>
27 
28 #include <cmath>
29 #include <cstdlib>
30 #include <iostream>
31 #include <sstream>
32 #include <utility> // for pair, make_pair
33 
34 class G4VSolid;
35 class PHCompositeNode;
36 
37 using namespace std;
38 
39 //_______________________________________________________________________
41  : PHG4Detector(subsys, Node, dnam)
42  , m_SuperDetector("NONE")
43  , m_Params(parameters)
44  , m_DisplayAction(dynamic_cast<PHG4CrystalCalorimeterDisplayAction*>(subsys->GetDisplayAction()))
45  , _towerlogicnameprefix("CrystalCalorimeterTower")
46  , m_IsActive(m_Params->get_int_param("active"))
47  , m_AbsorberActive(m_Params->get_int_param("absorberactive"))
48 {
49 }
50 
51 //_______________________________________________________________________
53 {
54  if (m_IsActive)
55  {
56  if (m_ActiveVolumeSet.find(volume) != m_ActiveVolumeSet.end())
57  {
58  return GetCaloType();
59  }
60  }
61  if (m_AbsorberActive)
62  {
63  if (m_PassiveVolumeSet.find(volume) != m_PassiveVolumeSet.end())
64  {
65  return -1;
66  }
67  }
68  return 0;
69 }
70 
71 //_______________________________________________________________________
73 {
74  if (Verbosity() > 0)
75  {
76  cout << "PHG4CrystalCalorimeterDetector: Begin Construction" << endl;
77  }
78 
79  if (m_Params->get_string_param("mappingtower").empty())
80  {
81  cout << "ERROR in PHG4CrystalCalorimeterDetector: No tower mapping file specified. Abort detector construction." << endl;
82  cout << "Please run set_string_param(\"mappingtower\", std::string filename ) first." << endl;
83  exit(1);
84  }
85 
86  /* Read parameters for detector construction and mapping from file */
88 
89  /* Create the cone envelope = 'world volume' for the crystal calorimeter */
92 
93  G4VSolid* eemc_envelope_solid = new G4Cons("eemc_envelope_solid",
94  m_Params->get_double_param("rMin1") * cm, m_Params->get_double_param("rMax1") * cm,
95  m_Params->get_double_param("rMin2") * cm, m_Params->get_double_param("rMax2") * cm,
96  m_Params->get_double_param("dz") * cm / 2.,
97  0, 2 * M_PI);
98 
99  G4LogicalVolume* eemc_envelope_log = new G4LogicalVolume(eemc_envelope_solid, WorldMaterial, G4String("eemc_envelope"), 0, 0, 0);
100 
101  GetDisplayAction()->AddVolume(eemc_envelope_log, "Envelope");
102  /* Define rotation attributes for envelope cone */
103  G4RotationMatrix eemc_rotm;
104  eemc_rotm.rotateX(m_Params->get_double_param("rot_x") * deg);
105  eemc_rotm.rotateY(m_Params->get_double_param("rot_y") * deg);
106  eemc_rotm.rotateZ(m_Params->get_double_param("rot_z") * deg);
107 
108  /* Place envelope cone in simulation */
109  // ostringstream name_envelope;
110  // name_envelope.str("");
111  string name_envelope = _towerlogicnameprefix + "_envelope";
112 
113  new G4PVPlacement(G4Transform3D(eemc_rotm, G4ThreeVector(m_Params->get_double_param("place_x") * cm, m_Params->get_double_param("place_y") * cm, m_Params->get_double_param("place_z") * cm)),
114  eemc_envelope_log, name_envelope, logicWorld, 0, false, OverlapCheck());
115 
116  /* Construct single calorimeter tower */
117  G4LogicalVolume* singletower = ConstructTower();
118 
119  /* Place calorimeter tower within envelope */
120  PlaceTower(eemc_envelope_log, singletower);
121 
122  return;
123 }
124 
125 //_______________________________________________________________________
128 {
129  if (Verbosity() > 0)
130  {
131  cout << "PHG4CrystalCalorimeterDetector: Build logical volume for single tower..." << endl;
132  }
133 
134  G4double carbon_thickness = 0.009 * cm;
135  G4double airgap_crystal_carbon = 0.012 * cm;
136 
137  /* dimesnions of full tower */
138  G4double tower_dx = m_Params->get_double_param("crystal_dx") * cm + 2 * (carbon_thickness + airgap_crystal_carbon);
139  G4double tower_dy = m_Params->get_double_param("crystal_dy") * cm + 2 * (carbon_thickness + airgap_crystal_carbon);
140  G4double tower_dz = m_Params->get_double_param("crystal_dz") * cm;
141 
142  /* create logical volume for single tower */
145 
146  G4VSolid* single_tower_solid = new G4Box(G4String("single_tower_solid"),
147  tower_dx / 2.0,
148  tower_dy / 2.0,
149  tower_dz / 2.0);
150 
151  G4LogicalVolume* single_tower_logic = new G4LogicalVolume(single_tower_solid,
152  WorldMaterial,
153  "single_tower_logic",
154  0, 0, 0);
155 
156  /* create geometry volume for crystal inside single_tower */
157  G4VSolid* solid_crystal = new G4Box(G4String("single_crystal_solid"),
158  m_Params->get_double_param("crystal_dx") * cm / 2.0,
159  m_Params->get_double_param("crystal_dy") * cm / 2.0,
160  m_Params->get_double_param("crystal_dz") * cm / 2.0);
161 
162  /* create geometry volume for frame (carbon fiber shell) inside single_tower */
163  G4VSolid* Carbon_hunk_solid = new G4Box(G4String("Carbon_hunk_solid"),
164  tower_dx / 2.0,
165  tower_dy / 2.0,
166  ((tower_dz / 2.0) - 1 * mm));
167 
168  G4double lead_dx = tower_dx - 2.0 * carbon_thickness;
169  G4double lead_dy = tower_dy - 2.0 * carbon_thickness;
170 
171  G4VSolid* lead_solid = new G4Box(G4String("lead_solid"),
172  lead_dx / 2.0,
173  lead_dy / 2.0,
174  tower_dz / 2.0);
175 
176  G4SubtractionSolid* Carbon_shell_solid = new G4SubtractionSolid(G4String("Carbon_Shell_solid"),
177  Carbon_hunk_solid,
178  lead_solid,
179  0,
180  G4ThreeVector(0.00 * mm, 0.00 * mm, 0.00 * mm));
181 
182  /* create logical volumes for crystal inside single_tower */
183  G4Material* material_crystal = GetDetectorMaterial(m_Params->get_string_param("material"));
184 
185  G4LogicalVolume* logic_crystal = new G4LogicalVolume(solid_crystal,
186  material_crystal,
187  "single_crystal_logic",
188  0, 0, 0);
189 
190  GetDisplayAction()->AddVolume(logic_crystal, "Crystal");
191 
192  /* create logical volumes for structural frame */
193  //Carbon Fiber
194  /*
195  G4double a = 12.01 * g / mole;
196  G4Element* elC = new G4Element("Carbon", "C", 6., a);
197  G4double density_carbon_fiber = 0.144 * g / cm3;
198  G4Material* CarbonFiber = new G4Material("CarbonFiber", density_carbon_fiber, 1);
199  CarbonFiber->AddElement(elC, 1);
200 */
201  G4Material* material_shell = GetCarbonFiber();
202 
203  G4LogicalVolume* logic_shell = new G4LogicalVolume(Carbon_shell_solid,
204  material_shell,
205  "single_absorber_logic",
206  0, 0, 0);
207 
208  GetDisplayAction()->AddVolume(logic_shell, "CarbonShell");
209 
210  /* Place structural frame in logical tower volume */
211  // ostringstream name_shell;
212  // name_shell.str("");
213  // name_shell << _towerlogicnameprefix << "_single_absorber";
214  string name_shell = _towerlogicnameprefix + "_single_absorber";
215  G4VPhysicalVolume* physvol = new G4PVPlacement(0, G4ThreeVector(0, 0, 0),
216  logic_shell,
217  name_shell,
218  single_tower_logic,
219  0, 0, OverlapCheck());
220  m_PassiveVolumeSet.insert(physvol);
221  /* Place crystal in logical tower volume */
222  string name_crystal = _towerlogicnameprefix + "_single_crystal";
223 
224  physvol = new G4PVPlacement(0, G4ThreeVector(0, 0, 0),
225  logic_crystal,
226  name_crystal,
227  single_tower_logic,
228  0, 0, OverlapCheck());
229  m_ActiveVolumeSet.insert(physvol);
230  if (Verbosity() > 0)
231  {
232  cout << "PHG4CrystalCalorimeterDetector: Building logical volume for single tower done." << endl;
233  }
234 
235  return single_tower_logic;
236 }
237 
239 {
240  /* Loop over all tower positions in vector and place tower */
241  typedef std::map<std::string, towerposition>::iterator it_type;
242 
243  for (it_type iterator = _map_tower.begin(); iterator != _map_tower.end(); ++iterator)
244  {
245  if (Verbosity() > 0)
246  {
247  cout << "PHG4CrystalCalorimeterDetector: Place tower " << iterator->first
248  << " idx_j = " << iterator->second.idx_j << ", idx_k = " << iterator->second.idx_k
249  << " at x = " << iterator->second.x << " , y = " << iterator->second.y << " , z = " << iterator->second.z << endl;
250  }
251  int copyno = (iterator->second.idx_j << 16) + iterator->second.idx_k;
252  new G4PVPlacement(0, G4ThreeVector(iterator->second.x, iterator->second.y, iterator->second.z),
253  singletower,
254  iterator->first,
255  eemcenvelope,
256  0, copyno, OverlapCheck());
257  }
258 
259  return 0;
260 }
261 
263 {
264  /* Open the datafile, if it won't open return an error */
265  ifstream istream_mapping(m_Params->get_string_param("mappingtower"));
266  if (!istream_mapping.is_open())
267  {
268  cout << "ERROR in PHG4CrystalCalorimeterDetector: Failed to open mapping file " << m_Params->get_string_param("mappingtower") << endl;
269  gSystem->Exit(1);
270  }
271 
272  /* loop over lines in file */
273  string line_mapping;
274  while (getline(istream_mapping, line_mapping))
275  {
276  /* Skip lines starting with / including a '#' */
277  if (line_mapping.find("#") != string::npos)
278  {
279  if (Verbosity() > 0)
280  {
281  cout << "PHG4CrystalCalorimeterDetector: SKIPPING line in mapping file: " << line_mapping << endl;
282  }
283  continue;
284  }
285 
286  istringstream iss(line_mapping);
287 
288  /* If line starts with keyword Tower, add to tower positions */
289  if (line_mapping.find("Tower ") != string::npos)
290  {
291  unsigned idx_j, idx_k, idx_l;
292  G4double pos_x, pos_y, pos_z;
293  G4double size_x, size_y, size_z;
294  G4double rot_x, rot_y, rot_z;
295  G4double dummy;
296  string dummys;
297 
298  /* read string- break if error */
299  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))
300  {
301  cout << "ERROR in PHG4CrystalCalorimeterDetector: Failed to read line in mapping file " << m_Params->get_string_param("mappingtower") << endl;
302  gSystem->Exit(1);
303  }
304 
305  /* Construct unique name for tower */
306  /* Mapping file uses cm, this class uses mm for length */
307  ostringstream towername;
308  towername.str("");
309  towername << _towerlogicnameprefix << "_j_" << idx_j << "_k_" << idx_k;
310 
311  /* Add Geant4 units */
312  pos_x = pos_x * cm;
313  pos_y = pos_y * cm;
314  pos_z = pos_z * cm;
315 
316  /* insert tower into tower map */
317  towerposition tower_new;
318  tower_new.x = pos_x;
319  tower_new.y = pos_y;
320  tower_new.z = pos_z;
321  tower_new.idx_j = idx_j;
322  tower_new.idx_k = idx_k;
323  _map_tower.insert(make_pair(towername.str(), tower_new));
324  }
325  else
326  {
327  /* If this line is not a comment and not a tower, save parameter as string / value. */
328  string parname;
329  G4double parval;
330 
331  /* read string- break if error */
332  if (!(iss >> parname >> parval))
333  {
334  cout << "ERROR in PHG4CrystalCalorimeterDetector: Failed to read line in mapping file " << m_Params->get_string_param("mappingtower") << endl;
335  gSystem->Exit(1);
336  }
337 
338  _map_global_parameter.insert(make_pair(parname, parval));
339  }
340  }
341 
342  /* Update member variables for global parameters based on parsed parameter file */
343  std::map<string, G4double>::iterator parit;
344 
345  parit = _map_global_parameter.find("Gcrystal_dx");
346  if (parit != _map_global_parameter.end())
347  {
348  m_Params->set_double_param("crystal_dx", parit->second); // in cm
349  }
350 
351  parit = _map_global_parameter.find("Gcrystal_dy");
352  if (parit != _map_global_parameter.end())
353  {
354  m_Params->set_double_param("crystal_dy", parit->second); // in cm
355  }
356 
357  parit = _map_global_parameter.find("Gcrystal_dz");
358  if (parit != _map_global_parameter.end())
359  {
360  m_Params->set_double_param("crystal_dz", parit->second); // in cm
361  }
362 
363  parit = _map_global_parameter.find("Gr1_inner");
364  if (parit != _map_global_parameter.end())
365  {
366  m_Params->set_double_param("rMin1", parit->second);
367  }
368 
369  parit = _map_global_parameter.find("Gr1_outer");
370  if (parit != _map_global_parameter.end())
371  {
372  m_Params->set_double_param("rMax1", parit->second);
373  }
374 
375  parit = _map_global_parameter.find("Gr2_inner");
376  if (parit != _map_global_parameter.end())
377  {
378  m_Params->set_double_param("rMin2", parit->second);
379  }
380 
381  parit = _map_global_parameter.find("Gr2_outer");
382  if (parit != _map_global_parameter.end())
383  {
384  m_Params->set_double_param("rMax2", parit->second);
385  }
386 
387  parit = _map_global_parameter.find("Gdz");
388  if (parit != _map_global_parameter.end())
389  {
390  m_Params->set_double_param("dz", parit->second);
391  }
392 
393  parit = _map_global_parameter.find("Gx0");
394  if (parit != _map_global_parameter.end())
395  {
396  m_Params->set_double_param("place_x", parit->second);
397  }
398 
399  parit = _map_global_parameter.find("Gy0");
400  if (parit != _map_global_parameter.end())
401  {
402  m_Params->set_double_param("place_y", parit->second);
403  }
404 
405  parit = _map_global_parameter.find("Gz0");
406  if (parit != _map_global_parameter.end())
407  {
408  m_Params->set_double_param("place_z", parit->second);
409  }
410 
411  parit = _map_global_parameter.find("Grot_x");
412  if (parit != _map_global_parameter.end())
413  {
414  m_Params->set_double_param("rot_x", parit->second * rad / deg);
415  }
416 
417  parit = _map_global_parameter.find("Grot_y");
418  if (parit != _map_global_parameter.end())
419  {
420  m_Params->set_double_param("rot_y", parit->second * rad / deg);
421  }
422  parit = _map_global_parameter.find("Grot_z");
423  if (parit != _map_global_parameter.end())
424  {
425  m_Params->set_double_param("rot_z", parit->second * rad / deg);
426  }
427  return 0;
428 }
429 
431 {
432  static string matname = "CrystalCarbonFiber";
433  G4Material* carbonfiber = GetDetectorMaterial(matname, false); // false suppresses warning that material does not exist
434  if (!carbonfiber)
435  {
436  G4double density_carbon_fiber = 1.44 * g / cm3;
437  carbonfiber = new G4Material(matname, density_carbon_fiber, 1);
438  carbonfiber->AddElement(G4Element::GetElement("C"), 1);
439  }
440  return carbonfiber;
441 }