Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4ZDCDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4ZDCDetector.cc
1 #include "PHG4ZDCDetector.h"
2 
3 #include "PHG4ZDCDefs.h"
4 #include "PHG4ZDCDisplayAction.h"
5 
6 #include <phparameter/PHParameters.h>
7 
8 #include <g4gdml/PHG4GDMLUtility.hh>
9 
10 #include <g4main/PHG4Detector.h> // for PHG4Detector
11 #include <g4main/PHG4DisplayAction.h> // for PHG4DisplayAction
12 #include <g4main/PHG4Subsystem.h>
13 
14 #include <phool/recoConsts.h>
15 
16 #include <TSystem.h>
17 
18 #include <Geant4/G4Box.hh>
19 #include <Geant4/G4LogicalVolume.hh>
20 #include <Geant4/G4PVPlacement.hh>
21 #include <Geant4/G4PhysicalConstants.hh>
22 #include <Geant4/G4RotationMatrix.hh> // for G4RotationMatrix
23 #include <Geant4/G4String.hh> // for G4String
24 #include <Geant4/G4SubtractionSolid.hh>
25 #include <Geant4/G4SystemOfUnits.hh>
26 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
27 #include <Geant4/G4Transform3D.hh> // for G4Transform3D
28 #include <Geant4/G4Tubs.hh>
29 #include <Geant4/G4Types.hh> // for G4double, G4int
30 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
31 
32 #include <cassert>
33 #include <cmath>
34 #include <iostream>
35 #include <sstream>
36 
37 class G4Material;
38 class G4VSolid;
39 class PHCompositeNode;
40 
41 //_______________________________________________________________________
42 PHG4ZDCDetector::PHG4ZDCDetector(PHG4Subsystem* subsys, PHCompositeNode* Node, PHParameters* parameters, const std::string& dnam, const int detid)
43  : PHG4Detector(subsys, Node, dnam)
44  , m_DisplayAction(dynamic_cast<PHG4ZDCDisplayAction*>(subsys->GetDisplayAction()))
45  , m_Params(parameters)
46  , m_GdmlConfig(PHG4GDMLUtility::GetOrMakeConfigNode(Node))
47  , m_Angle(M_PI_4 * radian)
48  , m_TPlate(2.3 * mm)
49  , m_HPlate(400.0 * mm)
50  , m_WPlate(100.0 * mm)
51  , m_TAbsorber(5.0 * mm)
52  , m_HAbsorber(150.0 * mm)
53  , m_WAbsorber(100.0 * mm)
54  , m_DFiber(0.5 * mm)
55  , m_HFiber(400.0 * mm)
56  , m_WFiber(100.0 * mm)
57  , m_GFiber(0.0001 * mm)
58  , m_Gap(0.2 * mm)
59  , m_TSMD(10.0 * mm)
60  , m_HSMD(160.0 * mm)
61  , m_WSMD(105.0 * mm)
62  , m_RHole(63.9 * mm)
63  , m_TWin(4.7 * mm)
64  , m_RWin(228.60 * mm)
65  , m_PlaceHole(122.56 * mm)
66  , m_Pxwin(0.0 * mm)
67  , m_Pywin(0.0 * mm)
68  , m_Pzwin(915.5 * mm)
69  , m_NMod(3)
70  , m_NLay(27)
71  , m_ActiveFlag(m_Params->get_int_param("active"))
72  , m_AbsorberActiveFlag(m_Params->get_int_param("absorberactive"))
73  , m_SupportActiveFlag(m_Params->get_int_param("supportactive"))
74  , m_Layer(detid)
75  , m_SuperDetector("NONE")
76 {
77  assert(m_GdmlConfig);
78 }
79 
80 //_______________________________________________________________________
81 int PHG4ZDCDetector::IsInZDC(G4VPhysicalVolume* volume) const
82 {
83  G4LogicalVolume* mylogvol = volume->GetLogicalVolume();
84 
85  if (m_ActiveFlag)
86  {
87  if (m_ScintiLogicalVolSet.find(mylogvol) != m_ScintiLogicalVolSet.end())
88  {
89  return 1;
90  }
91  if (m_FiberLogicalVolSet.find(mylogvol) != m_FiberLogicalVolSet.end())
92  {
93  return 2;
94  }
95  }
96  if (m_AbsorberActiveFlag)
97  {
98  if (m_AbsorberLogicalVolSet.find(mylogvol) != m_AbsorberLogicalVolSet.end())
99  {
100  return -1;
101  }
102  }
103  if (m_SupportActiveFlag)
104  {
105  if (m_SupportLogicalVolSet.find(mylogvol) != m_SupportLogicalVolSet.end())
106  {
107  return -2;
108  }
109  }
110  return 0;
111 }
112 
113 //_______________________________________________________________________
114 void PHG4ZDCDetector::ConstructMe(G4LogicalVolume* logicWorld)
115 {
116  if (Verbosity() > 0)
117  {
118  std::cout << "PHG4ZDCDetector: Begin Construction" << std::endl;
119  }
120 
121  if (m_Layer != PHG4ZDCDefs::NORTH &&
122  m_Layer != PHG4ZDCDefs::SOUTH)
123  {
124  std::cout << "use either PHG4ZDCDefs::NORTH or PHG4ZDCDefs::SOUTH for ZDC Subsystem" << std::endl;
125  gSystem->Exit(1);
126  return;
127  }
128 
129  recoConsts* rc = recoConsts::instance();
130  G4Material* WorldMaterial = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
131  G4Material* Fe = GetDetectorMaterial("G4_Fe");
132  G4Material* W = GetDetectorMaterial("G4_W");
133  G4Material* PMMA = GetDetectorMaterial("G4_PLEXIGLASS");
134  G4Material* Scint = GetDetectorMaterial("G4_POLYSTYRENE"); //place holder
135 
136  G4double TGap = m_DFiber + m_Gap;
137  G4double Mod_Length = 2 * m_TPlate + m_NLay * (TGap + m_TAbsorber);
138  G4double Det_Length = m_NMod * Mod_Length + m_TSMD;
139  G4double RTT = 1. / sin(m_Angle);
140  G4double First_Pos = -RTT * Det_Length / 2;
141  G4double Room = 3.5 * mm;
142  G4double Mother_2Z = RTT * Det_Length + 2. * (m_HFiber - m_HAbsorber / 2.) * cos(m_Angle);
143  /* Create exit windows */
144  G4VSolid* ExitWindow_nocut_solid = new G4Tubs(G4String("ExitWindow_nocut_solid"),
145  0.0, m_RWin, m_TWin, 0.0, CLHEP::twopi);
146 
147  G4VSolid* Hole_solid = new G4Tubs(G4String("Hole_solid"),
148  0.0, m_RHole, 2 * m_TWin, 0.0, CLHEP::twopi);
149  G4VSolid* ExitWindow_1cut_solid = new G4SubtractionSolid("ExitWindow_1cut_solid", ExitWindow_nocut_solid, Hole_solid, nullptr, G4ThreeVector(m_PlaceHole, 0, 0));
150 
151  G4VSolid* ExitWindow_2cut_solid = new G4SubtractionSolid("ExitWindow_2cut_solid", ExitWindow_1cut_solid, Hole_solid, nullptr, G4ThreeVector(-m_PlaceHole, 0, 0));
152 
153  G4LogicalVolume* ExitWindow_log = new G4LogicalVolume(ExitWindow_2cut_solid, GetDetectorMaterial("G4_STAINLESS-STEEL"), G4String("ExitWindow_log"), nullptr, nullptr, nullptr);
154 
155  GetDisplayAction()->AddVolume(ExitWindow_log, "Window");
156  m_SupportLogicalVolSet.insert(ExitWindow_log);
157  G4RotationMatrix Window_rotm;
158  Window_rotm.rotateX(m_Params->get_double_param("rot_x") * deg);
159  Window_rotm.rotateY(m_Params->get_double_param("rot_y") * deg);
160  Window_rotm.rotateZ(m_Params->get_double_param("rot_z") * deg);
161 
162  if (m_Layer == PHG4ZDCDefs::NORTH)
163  {
164  new G4PVPlacement(G4Transform3D(Window_rotm, G4ThreeVector(m_Pxwin, m_Pywin, m_Params->get_double_param("place_z") * cm - m_Pzwin)),
165  ExitWindow_log, "Window_North", logicWorld, false, PHG4ZDCDefs::NORTH, OverlapCheck());
166  }
167 
168  else if (m_Layer == PHG4ZDCDefs::SOUTH)
169  {
170  new G4PVPlacement(G4Transform3D(Window_rotm, G4ThreeVector(m_Pxwin, m_Pywin, -m_Params->get_double_param("place_z") * cm + m_Pzwin)),
171  ExitWindow_log, "Window_South", logicWorld, false, PHG4ZDCDefs::SOUTH, OverlapCheck());
172  }
173  /* ZDC detector here */
174  /* Create the box envelope = 'world volume' for ZDC */
175 
176  G4double Mother_X = m_WAbsorber / 2. + Room;
177  G4double Mother_Y = (m_HFiber - m_HAbsorber / 2.) * sin(m_Angle) + Room;
178  G4double Mother_Z = Mother_2Z / 2. + Room;
179  G4VSolid* ZDC_envelope_solid = new G4Box(G4String("ZDC_envelope_solid"),
180  Mother_X,
181  Mother_Y,
182  Mother_Z);
183 
184  G4LogicalVolume* ZDC_envelope_log = new G4LogicalVolume(ZDC_envelope_solid, WorldMaterial, G4String("ZDC_envelope_log"), nullptr, nullptr, nullptr);
185 
186  /* Define visualization attributes */
187  GetDisplayAction()->AddVolume(ZDC_envelope_log, "Envelope");
188 
189  /* Define rotation attributes for envelope cone */
190  G4RotationMatrix ZDC_rotm;
191  ZDC_rotm.rotateX(m_Params->get_double_param("rot_x") * deg);
192  ZDC_rotm.rotateY(m_Params->get_double_param("rot_y") * deg);
193  ZDC_rotm.rotateZ(m_Params->get_double_param("rot_z") * deg);
194 
195  /* Create logical volumes for a plate to contain fibers */
196  G4VSolid* fiber_plate_solid = new G4Box(G4String("fiber_plate_solid"),
197  m_WFiber / 2.,
198  m_HFiber / 2.,
199  TGap / 2.);
200 
201  G4LogicalVolume* fiber_plate_log = new G4LogicalVolume(fiber_plate_solid, WorldMaterial, G4String("fiber_plate_log"), nullptr, nullptr, nullptr);
202  GetDisplayAction()->AddVolume(fiber_plate_log, "fiber_plate_air");
203  /* front and back plate */
204  G4VSolid* fb_plate_solid = new G4Box(G4String("fb_plate_solid"),
205  m_WPlate / 2.,
206  m_HPlate / 2.,
207  m_TPlate / 2.);
208 
209  G4LogicalVolume* fb_plate_log = new G4LogicalVolume(fb_plate_solid, Fe, G4String("fb_plate_log"), nullptr, nullptr, nullptr);
210  m_SupportLogicalVolSet.insert(fb_plate_log);
211  GetDisplayAction()->AddVolume(fb_plate_log, "FrontBackPlate");
212 
213  /* absorber */
214  G4VSolid* absorber_solid = new G4Box(G4String("absorber_solid"),
215  m_WAbsorber / 2.,
216  m_HAbsorber / 2.,
217  m_TAbsorber / 2.);
218 
219  G4LogicalVolume* absorber_log = new G4LogicalVolume(absorber_solid, W, G4String("absorber_log"), nullptr, nullptr, nullptr);
220  m_AbsorberLogicalVolSet.insert(absorber_log);
221  GetDisplayAction()->AddVolume(absorber_log, "Absorber");
222 
223  /* SMD */
224  //volume that contains scintillators
225  G4VSolid* SMD_solid = new G4Box(G4String("SMD_solid"),
226  m_WSMD / 2.,
227  m_HSMD / 2.,
228  m_TSMD / 2.);
229 
230  G4LogicalVolume* SMD_log = new G4LogicalVolume(SMD_solid, WorldMaterial, G4String("SMD_log"), nullptr, nullptr, nullptr);
231  GetDisplayAction()->AddVolume(SMD_log, "SMD");
232  // small scintillators block
233  G4double scintx = 15 * mm;
234  G4double scinty = 20 * mm;
235  G4VSolid* Scint_solid = new G4Box(G4String("Scint_solid"),
236  scintx / 2.,
237  scinty / 2.,
238  m_TSMD / 2.);
239 
240  G4LogicalVolume* Scint_log = new G4LogicalVolume(Scint_solid, Scint, G4String("Scint_log"), nullptr, nullptr, nullptr);
241  m_ScintiLogicalVolSet.insert(Scint_log);
242  GetDisplayAction()->AddVolume(Scint_log, "Scint_solid");
243 
244  //put scintillators in the SMD volume
245  double scint_XPos = -m_WSMD / 2.;
246  double scint_Xstep = scintx / 2.;
247  double scint_Ystep = scinty / 2.;
248  int Nx = m_WSMD / scintx;
249  int Ny = m_HSMD / scinty;
250 
251  for (int i = 0; i < Nx; i++)
252  {
253  scint_XPos += scint_Xstep;
254  double scint_YPos = -m_HSMD / 2.;
255  for (int j = 0; j < Ny; j++)
256  {
257  int copyno = Nx * j + i;
258  scint_YPos += scint_Ystep;
259  G4RotationMatrix SMDRotation;
260  new G4PVPlacement(G4Transform3D(SMDRotation, G4ThreeVector(scint_XPos, scint_YPos, 0.0)),
261  Scint_log,
262  "single_scint",
263  SMD_log,
264  false, copyno, OverlapCheck());
265 
266  scint_YPos += scint_Ystep;
267  }
268  scint_XPos += scint_Xstep;
269  }
270 
271  /* Create logical volumes for fibers */
272  G4VSolid* single_fiber_solid = new G4Tubs(G4String("single_fiber_solid"),
273  0.0, (m_DFiber / 2.) - m_GFiber, (m_HFiber / 2.), 0.0, CLHEP::twopi);
274 
275  G4LogicalVolume* single_fiber_log = new G4LogicalVolume(single_fiber_solid, PMMA, G4String("single_fiber_log"), nullptr, nullptr, nullptr);
276  m_FiberLogicalVolSet.insert(single_fiber_log);
277  GetDisplayAction()->AddVolume(single_fiber_log, "Fiber");
278 
279  /* Rotation Matrix for fibers */
280  G4RotationMatrix* FiberRotation = new G4RotationMatrix();
281  FiberRotation->rotateX(90. * deg);
282 
283  /* Place fibers in the fiber plate */
284  G4double fiber_XPos = -m_WFiber / 2.;
285  G4double fiber_step = m_DFiber / 2.;
286  int Nfiber = m_WFiber / m_DFiber;
287 
288  for (int i = 0; i < Nfiber; i++)
289  {
290  fiber_XPos += fiber_step;
291  int copyno = i;
292 
293  new G4PVPlacement(FiberRotation, G4ThreeVector(fiber_XPos, 0.0, 0.0),
294  single_fiber_log,
295  G4String("single_fiber_scint"),
296  fiber_plate_log,
297  false, copyno, OverlapCheck());
298  fiber_XPos += fiber_step;
299  }
300 
301  /* Rotation for plates in ZDC */
302  G4RotationMatrix* PlateRotation = new G4RotationMatrix();
303 
304  PlateRotation->rotateX(-m_Angle);
305 
306  /* construct ZDC */
307  G4double Plate_Step = m_TPlate * RTT;
308  G4double Absorber_Step = m_TAbsorber * RTT;
309  G4double Gap_Step = TGap * RTT;
310  G4double SMD_Step = m_TSMD * RTT;
311  G4double ZPos = First_Pos - Plate_Step / 2.;
312  G4double Gap_YPos = (m_HFiber - m_HAbsorber) / 2. * sin(m_Angle);
313  G4double Gap_ZPos = (m_HFiber - m_HAbsorber) / 2. * cos(m_Angle);
314 
315  G4double Plate_YPos = (m_HPlate - m_HAbsorber) / 2. * sin(m_Angle);
316  G4double Plate_ZPos = (m_HPlate - m_HAbsorber) / 2. * cos(m_Angle);
317 
318  G4double SMD_YPos = (m_HSMD - m_HAbsorber) / 2. * sin(m_Angle);
319  G4double SMD_ZPos = (m_HSMD - m_HAbsorber) / 2. * cos(m_Angle);
320 
321  /* start the loop: for every module --- front plate-absorber-fiber plate-absorber-.....-fiber plate-back plate */
322  int copyno_plate = 0;
323  for (int i = 0; i < m_NMod; i++)
324  {
325  //place the SMD in between the 1st and 2nd module
326  if (i == 1)
327  {
328  ZPos += (SMD_Step / 2.);
329  new G4PVPlacement(PlateRotation, G4ThreeVector(0.0, SMD_YPos, SMD_ZPos + ZPos),
330  SMD_log,
331  G4String("SMD"),
332  ZDC_envelope_log,
333  false, 0, OverlapCheck()); //using copy number for now, need to find a better way
334  ZPos += (SMD_Step / 2.);
335  }
336  /* place the front plate */
337  ZPos += (Plate_Step / 2.);
338  new G4PVPlacement(PlateRotation, G4ThreeVector(0.0, Plate_YPos, Plate_ZPos + ZPos),
339  fb_plate_log,
340  G4String("front_plate"),
341  ZDC_envelope_log,
342  false, copyno_plate, OverlapCheck());
343  ZPos += (Plate_Step / 2.);
344  copyno_plate++;
345  for (int j = 0; j < m_NLay; j++)
346  {
347  /* place the Absorber */
348  ZPos += (Absorber_Step / 2.);
349  new G4PVPlacement(PlateRotation, G4ThreeVector(0.0, 0.0, ZPos),
350  absorber_log,
351  G4String("single_absorber"),
352  ZDC_envelope_log,
353  false, i * 100 + j, OverlapCheck());
354  ZPos += (Absorber_Step / 2.);
355 
356  /* place the fiber plate */
357  ZPos += (Gap_Step / 2.);
358  int copyno = 27 * i + j;
359  std::string name_fiber_plate = "Fiber_Plate_" + std::to_string(copyno);
360  new G4PVPlacement(PlateRotation, G4ThreeVector(0.0, Gap_YPos, Gap_ZPos + ZPos),
361  fiber_plate_log,
362  name_fiber_plate,
363  ZDC_envelope_log,
364  false, copyno, OverlapCheck());
365  ZPos += (Gap_Step / 2.);
366  }
367  /* place the back plate */
368  ZPos += (Plate_Step / 2.);
369  new G4PVPlacement(PlateRotation, G4ThreeVector(0.0, Plate_YPos, Plate_ZPos + ZPos),
370  fb_plate_log,
371  G4String("back_plate"),
372  ZDC_envelope_log,
373  false, copyno_plate, OverlapCheck());
374  copyno_plate++;
375  ZPos += (Plate_Step / 2.);
376  }
377 
378  /* Place envelope cone in simulation */
379 
380  if (m_Layer == PHG4ZDCDefs::NORTH)
381  {
382  new G4PVPlacement(G4Transform3D(ZDC_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)),
383  ZDC_envelope_log, "ZDC_Envelope_North", logicWorld, false, PHG4ZDCDefs::NORTH, OverlapCheck());
384  }
385  else if (m_Layer == PHG4ZDCDefs::SOUTH)
386  {
387  ZDC_rotm.rotateY(180 * deg);
388  new G4PVPlacement(G4Transform3D(ZDC_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)),
389  ZDC_envelope_log, "ZDC_Envelope_South", logicWorld, false, PHG4ZDCDefs::SOUTH, OverlapCheck());
390  }
391  return;
392 }