EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4TpcDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4TpcDetector.cc
1 #include "PHG4TpcDetector.h"
2 #include "PHG4TpcDefs.h"
3 #include "PHG4TpcDisplayAction.h"
4 
5 #include <g4main/PHG4Detector.h> // for PHG4Detector
6 #include <g4main/PHG4DisplayAction.h> // for PHG4DisplayAction
7 #include <g4main/PHG4Subsystem.h>
8 
9 #include <TSystem.h>
10 
11 #include <phparameter/PHParameters.h>
12 
13 #include <phool/recoConsts.h>
14 
15 #include <Geant4/G4LogicalVolume.hh>
16 #include <Geant4/G4Material.hh>
17 #include <Geant4/G4PVPlacement.hh>
18 #include <Geant4/G4String.hh> // for G4String
19 #include <Geant4/G4SystemOfUnits.hh>
20 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
21 #include <Geant4/G4Tubs.hh>
22 #include <Geant4/G4UserLimits.hh>
23 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
24 
25 #include <cassert>
26 #include <cmath>
27 #include <iostream> // for basic_ostream::operator<<
28 #include <map> // for map
29 #include <sstream>
30 
31 class G4VSolid;
32 class PHCompositeNode;
33 
34 //_______________________________________________________________
35 PHG4TpcDetector::PHG4TpcDetector(PHG4Subsystem *subsys, PHCompositeNode *Node, PHParameters *parameters, const std::string &dnam)
36  : PHG4Detector(subsys, Node, dnam)
37  , m_DisplayAction(dynamic_cast<PHG4TpcDisplayAction *>(subsys->GetDisplayAction()))
38  , m_Params(parameters)
39  , m_ActiveFlag(m_Params->get_int_param("active"))
40  , m_AbsorberActiveFlag(m_Params->get_int_param("absorberactive"))
41  , m_InnerCageRadius(m_Params->get_double_param("gas_inner_radius") * cm - m_Params->get_double_param("cage_layer_9_thickness") * cm - m_Params->get_double_param("cage_layer_8_thickness") * cm - m_Params->get_double_param("cage_layer_7_thickness") * cm - m_Params->get_double_param("cage_layer_6_thickness") * cm - m_Params->get_double_param("cage_layer_5_thickness") * cm - m_Params->get_double_param("cage_layer_4_thickness") * cm - m_Params->get_double_param("cage_layer_3_thickness") * cm - m_Params->get_double_param("cage_layer_2_thickness") * cm - m_Params->get_double_param("cage_layer_1_thickness") * cm)
42  , m_OuterCageRadius(m_Params->get_double_param("gas_outer_radius") * cm + m_Params->get_double_param("cage_layer_9_thickness") * cm + m_Params->get_double_param("cage_layer_8_thickness") * cm + m_Params->get_double_param("cage_layer_7_thickness") * cm + m_Params->get_double_param("cage_layer_6_thickness") * cm + m_Params->get_double_param("cage_layer_5_thickness") * cm + m_Params->get_double_param("cage_layer_4_thickness") * cm + m_Params->get_double_param("cage_layer_3_thickness") * cm + m_Params->get_double_param("cage_layer_2_thickness") * cm + m_Params->get_double_param("cage_layer_1_thickness") * cm)
43 {
44 }
45 
46 //_______________________________________________________________
47 int PHG4TpcDetector::IsInTpc(G4VPhysicalVolume *volume) const
48 {
49  if (m_ActiveFlag)
50  {
51  if (m_ActiveVolumeSet.find(volume) != m_ActiveVolumeSet.end())
52  {
53  return 1;
54  }
55  }
56  if (m_AbsorberActiveFlag)
57  {
58  if (m_AbsorberVolumeSet.find(volume) != m_AbsorberVolumeSet.end())
59  {
60  return -1;
61  }
62  }
63  return 0;
64 }
65 
66 //_______________________________________________________________
67 void PHG4TpcDetector::ConstructMe(G4LogicalVolume *logicWorld)
68 {
69  // create Tpc envelope
70  // tpc consists of (from inside to gas volume, outside is reversed up to now)
71  // 1st layer cu
72  // 2nd layer FR4
73  // 3rd layer HoneyComb
74  // 4th layer cu
75  // 5th layer FR4
76  // 6th layer Kapton
77  // 7th layer cu
78  // 8th layer Kapton
79  // 9th layer cu
80 
81  double steplimits = m_Params->get_double_param("steplimits") * cm;
82  if (std::isfinite(steplimits))
83  {
84  m_G4UserLimits = new G4UserLimits(steplimits);
85  }
86 
87  G4VSolid *tpc_envelope = new G4Tubs("tpc_envelope", m_InnerCageRadius, m_OuterCageRadius, m_Params->get_double_param("tpc_length") * cm / 2., 0., 2 * M_PI);
88 
89  recoConsts *rc = recoConsts::instance();
90  G4LogicalVolume *tpc_envelope_logic = new G4LogicalVolume(tpc_envelope,
91  GetDetectorMaterial(rc->get_StringFlag("WorldMaterial")),
92  "tpc_envelope");
93  m_DisplayAction->AddVolume(tpc_envelope_logic, "TpcEnvelope");
94 
95  ConstructTpcExternalSupports(logicWorld);
96 
97  ConstructTpcCageVolume(tpc_envelope_logic);
98  ConstructTpcGasVolume(tpc_envelope_logic);
99 
100  new G4PVPlacement(0, 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),
101  tpc_envelope_logic, "tpc_envelope",
102  logicWorld, 0, false, OverlapCheck());
103 }
104 
105 int PHG4TpcDetector::ConstructTpcGasVolume(G4LogicalVolume *tpc_envelope)
106 {
107  static std::map<int, std::string> tpcgasvolname =
108  {{PHG4TpcDefs::North, "tpc_gas_north"},
109  {PHG4TpcDefs::South, "tpc_gas_south"}};
110 
111  // Window / central membrane
112  double tpc_window_thickness = m_Params->get_double_param("window_thickness") * cm;
113  double tpc_half_length = (m_Params->get_double_param("tpc_length") * cm - tpc_window_thickness) / 2.;
114 
115  //'window' (modernly called central membrane only) material is ENIG, not Copper:
116  //thickness in this recipe are just a ratio. we set the usual thickness below.
117  std::vector<double> thickness;
118  std::vector<std::string> material;
119  material.push_back("G4_Ni");
120  thickness.push_back(.240 * cm);
121  material.push_back("G4_Au");
122  thickness.push_back(.008 * cm);
123  G4Material *temp = nullptr;
124  temp=GetDetectorMaterial("ENIG", false);
125  if (temp == nullptr)
126  {
127  CreateCompositeMaterial("ENIG", material, thickness); //see new function below
128  }
129 
130 
131 
132 
133  G4VSolid *tpc_window = new G4Tubs("tpc_window", m_Params->get_double_param("gas_inner_radius") * cm, m_Params->get_double_param("gas_outer_radius") * cm, tpc_window_thickness / 2., 0., 2 * M_PI);
134  //we build our CM surface:
135  G4LogicalVolume *tpc_window_logic = new G4LogicalVolume(tpc_window,
136  GetDetectorMaterial("ENIG"),
137  "tpc_window");
138  //previously: GetDetectorMaterial(m_Params->get_string_param("window_surface1_material")),
139 
140 
141  // G4VisAttributes *visatt = new G4VisAttributes();
142  // visatt->SetVisibility(true);
143  // visatt->SetForceSolid(true);
144  // visatt->SetColor(PHG4TPCColorDefs::tpc_cu_color);
145  // tpc_window_logic->SetVisAttributes(visatt);
146 
147  m_DisplayAction->AddVolume(tpc_window_logic, "TpcWindow");
148  G4VPhysicalVolume *tpc_window_phys = new G4PVPlacement(0, G4ThreeVector(0, 0, 0),
149  tpc_window_logic, "tpc_window",
150  tpc_envelope, false, PHG4TpcDefs::Window, OverlapCheck());
151 
152  m_AbsorberVolumeSet.insert(tpc_window_phys);
153 
154  //now build the FR4 layer beneath that:
155  // Window / central membrane core
156  double tpc_window_surface1_thickness = m_Params->get_double_param("window_surface1_thickness") * cm;
157  double tpc_window_surface2_thickness = m_Params->get_double_param("window_surface2_thickness") * cm;
158  double tpc_window_surface2_core_thickness = tpc_window_thickness - 2 * tpc_window_surface1_thickness;
159  double tpc_window_core_thickness = tpc_window_surface2_core_thickness - 2 * (tpc_window_surface2_thickness);
160 
161  G4VSolid *tpc_window_surface2_core =
162  new G4Tubs("tpc_window_surface2_core", m_Params->get_double_param("gas_inner_radius") * cm, m_Params->get_double_param("gas_outer_radius") * cm,
163  tpc_window_surface2_core_thickness / 2., 0., 2 * M_PI);
164  G4LogicalVolume *tpc_window_surface2_core_logic = new G4LogicalVolume(tpc_window_surface2_core,
165  GetDetectorMaterial(m_Params->get_string_param("window_surface2_material")),
166  "tpc_window_surface2_core");
167 
168  m_DisplayAction->AddVolume(tpc_window_surface2_core_logic, "TpcWindow");
169  // visatt = new G4VisAttributes();
170  // visatt->SetVisibility(true);
171  // visatt->SetForceSolid(true);
172  // visatt->SetColor(PHG4TPCColorDefs::tpc_pcb_color);
173  // tpc_window_surface2_core_logic->SetVisAttributes(visatt);
174  G4VPhysicalVolume *tpc_window_surface2_core_phys = new G4PVPlacement(0, G4ThreeVector(0, 0, 0),
175  tpc_window_surface2_core_logic, "tpc_window_surface2_core",
176  tpc_window_logic, false, PHG4TpcDefs::WindowCore, OverlapCheck());
177 
178  m_AbsorberVolumeSet.insert(tpc_window_surface2_core_phys);
179 
180  //and now the honeycomb core:
181  G4VSolid *tpc_window_core =
182  new G4Tubs("tpc_window", m_Params->get_double_param("gas_inner_radius") * cm, m_Params->get_double_param("gas_outer_radius") * cm,
183  tpc_window_core_thickness / 2., 0., 2 * M_PI);
184  G4LogicalVolume *tpc_window_core_logic = new G4LogicalVolume(tpc_window_core,
185  GetDetectorMaterial(m_Params->get_string_param("window_core_material")),
186  "tpc_window_core");
187 
188  m_DisplayAction->AddVolume(tpc_window_core_logic, "TpcHoneyComb");
189  G4VPhysicalVolume *tpc_window_core_phys = new G4PVPlacement(0, G4ThreeVector(0, 0, 0),
190  tpc_window_core_logic, "tpc_window_core",
191  tpc_window_surface2_core_logic, false, PHG4TpcDefs::WindowCore, OverlapCheck());
192 
193  m_AbsorberVolumeSet.insert(tpc_window_core_phys);
194 
195  // Gas
196  G4VSolid *tpc_gas = new G4Tubs("tpc_gas", m_Params->get_double_param("gas_inner_radius") * cm, m_Params->get_double_param("gas_outer_radius") * cm, tpc_half_length / 2., 0., 2 * M_PI);
197 
198  G4LogicalVolume *tpc_gas_logic = new G4LogicalVolume(tpc_gas,
199  GetDetectorMaterial(m_Params->get_string_param("tpc_gas")),
200  "tpc_gas");
201 
202  tpc_gas_logic->SetUserLimits(m_G4UserLimits);
203  m_DisplayAction->AddVolume(tpc_gas_logic, "TpcGas");
204  G4VPhysicalVolume *tpc_gas_phys = new G4PVPlacement(0, G4ThreeVector(0, 0, (tpc_half_length + tpc_window_thickness) / 2.),
205  tpc_gas_logic, tpcgasvolname[PHG4TpcDefs::North],
206  tpc_envelope, false, PHG4TpcDefs::North, OverlapCheck());
207  std::cout << "north copy no: " << tpc_gas_phys->GetCopyNo() << std::endl;
208 
209  m_ActiveVolumeSet.insert(tpc_gas_phys);
210  tpc_gas_phys = new G4PVPlacement(0, G4ThreeVector(0, 0, -(tpc_half_length + tpc_window_thickness) / 2.),
211  tpc_gas_logic, tpcgasvolname[PHG4TpcDefs::South],
212  tpc_envelope, false, PHG4TpcDefs::South, OverlapCheck());
213 
214  std::cout << "south copy no: " << tpc_gas_phys->GetCopyNo() << std::endl;
215  m_ActiveVolumeSet.insert(tpc_gas_phys);
216 
217 #if G4VERSION_NUMBER >= 1033
218  const G4RegionStore *theRegionStore = G4RegionStore::GetInstance();
219  G4Region *tpcregion = theRegionStore->GetRegion("REGION_TPCGAS");
220  tpc_gas_logic->SetRegion(tpcregion);
221  tpcregion->AddRootLogicalVolume(tpc_gas_logic);
222 #endif
223  return 0;
224 }
225 
226 int PHG4TpcDetector::ConstructTpcExternalSupports(G4LogicalVolume *logicWorld)
227 {
228  //note that these elements are outside the tpc logical volume!
229 
230  // Two two-inch diam. 304 Stainless Steel solid 'hanger beams' at 32.05" from beam center
231  // at +/- 41.39 degrees left and right of vertical
232  //stainless steel: 0.695 iron, 0.190 chromium, 0.095 nickel, 0.020 manganese.
233  //if we're being pedantic, that is. But store-bought stainless is probably okay.
234  // G4Material *StainlessSteel = new G4Material("StainlessSteel", density = 8.02*g/cm3, 5);
235  //StainlessSteel->AddMaterial(matman->FindOrBuildMaterial("G4_Si"), 0.01);
236  //StainlessSteel->AddMaterial(matman->FindOrBuildMaterial("G4_Mn"), 0.02);
237  //StainlessSteel->AddMaterial(matman->FindOrBuildMaterial("G4_Cr"), 0.19);
238  //StainlessSteel->AddMaterial(matman->FindOrBuildMaterial("G4_Ni"), 0.10);
239  //StainlessSteel->AddMaterial(matman->FindOrBuildMaterial("G4_Fe"), 0.68);
240  G4Material *stainlessSteel=GetDetectorMaterial("G4_STAINLESS-STEEL");
241  double inch=2.54*cm;
242  double hangerAngle=41.39*M_PI/180.;
243  double hangerRadius=32.05*inch;
244  double hangerX=std::sin(hangerAngle)*hangerRadius;
245  double hangerY=std::cos(hangerAngle)*hangerRadius;
246  double hangerDiameter=2.*inch;
247  G4VSolid *hangerBeam = new G4Tubs("tpc_hanger_beam", 0,hangerDiameter/2.,(m_Params->get_double_param("tpc_length") * cm )/ 2., 0., 2*M_PI);
248  G4LogicalVolume *hangerBeamLogic = new G4LogicalVolume(hangerBeam,
249  stainlessSteel,
250  "tpc_hanger_beam");
251 
252  m_DisplayAction->AddVolume(hangerBeamLogic, "TpcHangerBeam");
253  G4VPhysicalVolume *tpc_hanger_beam_phys[2]={nullptr,nullptr};
254  tpc_hanger_beam_phys[0] = new G4PVPlacement(0, G4ThreeVector(hangerX,hangerY, 0),
255  hangerBeamLogic, "tpc_hanger_beam_left",
256  logicWorld, false, 0, OverlapCheck());
257  tpc_hanger_beam_phys[1] = new G4PVPlacement(0, G4ThreeVector(-hangerX,hangerY, 0),
258  hangerBeamLogic, "tpc_hanger_beam_right",
259  logicWorld, false, 1, OverlapCheck());
260  m_AbsorberVolumeSet.insert(tpc_hanger_beam_phys[0]);
261  m_AbsorberVolumeSet.insert(tpc_hanger_beam_phys[1]);
262 
263 
264 
265  //Twelve one-inch diam carbon fiber rods of thickness 1/16" at 31.04" from beam center
266  //borrowed from the INTT specification of carbon fiber
267  //note that this defines a clocking!
268  G4Material *carbonFiber=GetDetectorMaterial("CFRP_INTT");
269  double rodAngleStart=M_PI/12.;
270  double rodAngularSpacing=2*M_PI/12.;
271  double rodRadius=31.5*inch;
272  double rodWallThickness=1./8.*inch;
273  double rodDiameter=3./4.*inch;
274  G4VSolid *tieRod = new G4Tubs("tpc_tie_rod",rodDiameter/2.-rodWallThickness, rodDiameter/2.,(m_Params->get_double_param("tpc_length") * cm )/ 2., 0., 2*M_PI);
275  G4LogicalVolume *tieRodLogic = new G4LogicalVolume(tieRod,
276  carbonFiber,
277  "tpc_tie_rod");
278 
279  G4VPhysicalVolume *tpc_tie_rod_phys[12]={nullptr,nullptr,nullptr,nullptr,nullptr,nullptr,
280  nullptr,nullptr,nullptr,nullptr,nullptr,nullptr};
281 
282  std::ostringstream name;
283  for (int i=0;i<12;i++){
284  double ang=rodAngleStart+rodAngularSpacing*i;
285  name.str("");
286  name << "tpc_tie_rod_" << i;
287  tpc_tie_rod_phys[i] = new G4PVPlacement(0, G4ThreeVector(rodRadius*cos(ang),rodRadius*sin(ang), 0),
288  tieRodLogic, name.str(),
289  logicWorld, false, i, OverlapCheck());
290  m_AbsorberVolumeSet.insert(tpc_tie_rod_phys[i]);
291  }
292 
293 
294 
295  // G4VisAttributes *visatt = new G4VisAttributes();
296  // visatt->SetVisibility(true);
297  // visatt->SetForceSolid(true);
298  // visatt->SetColor(PHG4TPCColorDefs::tpc_cu_color);
299  // tpc_window_logic->SetVisAttributes(visatt);
300 
301 
302  return 0;
303 }
304 
305 int PHG4TpcDetector::ConstructTpcCageVolume(G4LogicalVolume *tpc_envelope)
306 {
307  // 8th layer cu
308  // 9th layer FR4
309  // 10th layer HoneyComb
310  // 11th layer FR4
311  // 12th layer Kapton
312  // 13th layer FR4
313  static const int nlayers = 9;
314  static const double thickness[nlayers] = {m_Params->get_double_param("cage_layer_1_thickness") * cm,
315  m_Params->get_double_param("cage_layer_2_thickness") * cm,
316  m_Params->get_double_param("cage_layer_3_thickness") * cm,
317  m_Params->get_double_param("cage_layer_4_thickness") * cm,
318  m_Params->get_double_param("cage_layer_5_thickness") * cm,
319  m_Params->get_double_param("cage_layer_6_thickness") * cm,
320  m_Params->get_double_param("cage_layer_7_thickness") * cm,
321  m_Params->get_double_param("cage_layer_8_thickness") * cm,
322  m_Params->get_double_param("cage_layer_9_thickness") * cm};
323 
324  static const std::string material[nlayers] = {m_Params->get_string_param("cage_layer_1_material"),
325  m_Params->get_string_param("cage_layer_2_material"),
326  m_Params->get_string_param("cage_layer_3_material"),
327  m_Params->get_string_param("cage_layer_4_material"),
328  m_Params->get_string_param("cage_layer_5_material"),
329  m_Params->get_string_param("cage_layer_6_material"),
330  m_Params->get_string_param("cage_layer_7_material"),
331  m_Params->get_string_param("cage_layer_8_material"),
332  m_Params->get_string_param("cage_layer_9_material")};
333 
334 
335 
336 
337 
338  double tpc_cage_radius = m_InnerCageRadius;
339  std::ostringstream name;
340  for (int i = 0; i < nlayers; i++)
341  {
342  name.str("");
343  int layerno = i + 1;
344  name << "tpc_cage_layer_" << layerno;
345  G4VSolid *tpc_cage_layer = new G4Tubs(name.str(), tpc_cage_radius, tpc_cage_radius + thickness[i], m_Params->get_double_param("tpc_length") * cm / 2., 0., 2 * M_PI);
346  G4LogicalVolume *tpc_cage_layer_logic = new G4LogicalVolume(tpc_cage_layer,
347  GetDetectorMaterial(material[i]),
348  name.str());
349  m_DisplayAction->AddTpcInnerLayer(tpc_cage_layer_logic);
350  G4VPhysicalVolume *tpc_cage_layer_phys = new G4PVPlacement(0, G4ThreeVector(0, 0, 0),
351  tpc_cage_layer_logic, name.str(),
352  tpc_envelope, false, layerno, OverlapCheck());
353  m_AbsorberVolumeSet.insert(tpc_cage_layer_phys);
354  tpc_cage_radius += thickness[i];
355  }
356  // outer cage
357 
358  tpc_cage_radius = m_OuterCageRadius;
359  for (int i = 0; i < nlayers; i++)
360  {
361  tpc_cage_radius -= thickness[i];
362  name.str("");
363  int layerno = 10 + 1 + i; // so the accompanying inner layer is layer - 10
364  name << "tpc_cage_layer_" << layerno;
365  G4VSolid *tpc_cage_layer = new G4Tubs(name.str(), tpc_cage_radius, tpc_cage_radius + thickness[i], m_Params->get_double_param("tpc_length") * cm / 2., 0., 2 * M_PI);
366  G4LogicalVolume *tpc_cage_layer_logic = new G4LogicalVolume(tpc_cage_layer,
367  GetDetectorMaterial(material[i]),
368  name.str());
369  m_DisplayAction->AddTpcOuterLayer(tpc_cage_layer_logic);
370  G4VPhysicalVolume *tpc_cage_layer_phys = new G4PVPlacement(0, G4ThreeVector(0, 0, 0),
371  tpc_cage_layer_logic, name.str(),
372  tpc_envelope, false, layerno, OverlapCheck());
373  m_AbsorberVolumeSet.insert(tpc_cage_layer_phys);
374  }
375 
376  return 0;
377 }
378 
379 
380 void PHG4TpcDetector ::CreateCompositeMaterial(
381  std::string compositeName,
382  std::vector<std::string> materialName,
383  std::vector<double> thickness)
384 {
385  //takes in a list of material names known to Geant already, and thicknesses, and creates a new material called compositeName.
386 
387  //check that desired material name doesn't already exist
388  //note that this throws a warning.
389  std::cout << __PRETTY_FUNCTION__ << " NOTICE: Checking if material " << compositeName << " exists. This will return a warning if it doesn't, but that is okay." << std::endl;
390  G4Material *tempmat = GetDetectorMaterial(compositeName, false);
391 
392  if (tempmat != nullptr)
393  {
394  std::cout << __PRETTY_FUNCTION__ << " Fatal Error: composite material " << compositeName << " already exists" << std::endl;
395  assert(!tempmat);
396  }
397 
398  //check that both arrays have the same depth
399  assert(materialName.size() == thickness.size());
400 
401  //sum up the areal density and total thickness so we can divvy it out
402  double totalArealDensity = 0, totalThickness = 0;
403  for (std::vector<double>::size_type i = 0; i < thickness.size(); i++)
404  {
405  tempmat = GetDetectorMaterial(materialName[i]);
406  if (tempmat == nullptr)
407  {
408  std::cout << __PRETTY_FUNCTION__ << " Fatal Error: component material " << materialName[i] << " does not exist." << std::endl;
409  gSystem->Exit(1);
410  exit(1);
411  }
412  totalArealDensity += tempmat->GetDensity() * thickness[i];
413  totalThickness += thickness[i];
414  }
415 
416  //register a new material with the average density of the whole:
417  double compositeDensity = totalArealDensity / totalThickness;
418  G4Material *composite = new G4Material(compositeName, compositeDensity, thickness.size());
419 
420  //now calculate the fraction due to each material, and register those
421  for (std::vector<double>::size_type i = 0; i < thickness.size(); i++)
422  {
423  tempmat = GetDetectorMaterial(materialName[i]); //don't need to check this, since we did in the previous loop.
424  composite->AddMaterial(tempmat, thickness[i] * tempmat->GetDensity() / totalArealDensity);
425  }
426 
427  //how to register our finished material?
428  return;
429 }