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