Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4CEmcTestBeamDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4CEmcTestBeamDetector.cc
2 
3 #include <g4main/PHG4Detector.h> // for PHG4Detector
4 
5 #include <phool/recoConsts.h>
6 
7 #include <Geant4/G4Box.hh>
8 #include <Geant4/G4Colour.hh>
9 #include <Geant4/G4LogicalVolume.hh>
10 #include <Geant4/G4PVPlacement.hh>
11 #include <Geant4/G4RotationMatrix.hh> // for G4RotationMatrix
12 #include <Geant4/G4String.hh> // for G4String
13 #include <Geant4/G4SystemOfUnits.hh>
14 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
15 #include <Geant4/G4Transform3D.hh> // for G4Transform3D
16 #include <Geant4/G4Tubs.hh>
17 #include <Geant4/G4Types.hh> // for G4double
18 #include <Geant4/G4VisAttributes.hh>
19 
20 #include <algorithm> // for copy
21 #include <cmath> // for cos, sin, NAN, acos, atan
22 #include <iostream> // for operator<<, ostringstream
23 #include <sstream>
24 
25 class G4Material;
26 class G4VSolid;
27 class PHCompositeNode;
28 
29 using namespace std;
30 
31 static double no_overlap = 0.0001 * cm; // added safety margin against overlaps by using same boundary between volumes
32 
34  : PHG4Detector(subsys, Node, dnam)
35  , gap(0.25 * mm)
36  , place_in_x(0 * cm)
37  , place_in_y(0 * cm)
38  , place_in_z(0 * cm)
39  , plate_x(135 * mm)
40  , plate_z(135 * mm)
41  , x_rot(0)
42  , y_rot(0)
43  , z_rot(0)
44  , alpha(NAN)
45  , inner_radius(NAN)
46  , outer_radius(NAN)
47  , tower_angular_coverage(NAN)
48  , cemc_angular_coverage(NAN)
49  , active_scinti_fraction(0.78)
50  , sandwiches_per_tower(12)
51  , // 12 tungsten/scintillator fiber snadwiches per tower
52  num_towers(7)
53  , active(0)
54  , absorberactive(0)
55  , layer(lyr)
56  , blackhole(0)
57 {
58  w_dimension[0] = plate_x;
59  w_dimension[1] = 0.5 * mm;
60  w_dimension[2] = plate_z;
61  sc_dimension[0] = plate_x;
62  sc_dimension[1] = 1 * mm;
63  sc_dimension[2] = plate_z;
64  sandwich_thickness = 2 * w_dimension[1] + sc_dimension[1]; // two tungsten plats, one scintillator
65  for (int i = 0; i < 4; i++)
66  {
67  sandwich_vol.push_back(nullptr);
68  }
69 }
70 
71 //_______________________________________________________________
72 //_______________________________________________________________
74 {
75  if (active && volume == sandwich_vol[2])
76  {
77  return 1;
78  }
79  if (absorberactive && (volume == sandwich_vol[0] || volume == sandwich_vol[1]))
80  {
81  return -1;
82  }
83  if (absorberactive && sandwich_vol[3] && volume == sandwich_vol[3])
84  {
85  return -1;
86  }
87  return 0;
88 }
89 
90 void PHG4CEmcTestBeamDetector::ConstructMe(G4LogicalVolume* logicWorld)
91 {
94  G4Material* Air = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
95  G4VSolid* cemc_tub = new G4Tubs("CEmcTub", inner_radius - 2 * no_overlap, outer_radius + 2 * no_overlap, (w_dimension[2] + 2 * no_overlap) / 2., 0, cemc_angular_coverage);
96  G4LogicalVolume* cemc_log = new G4LogicalVolume(cemc_tub, Air, G4String("CEmc"), nullptr, nullptr, nullptr);
97 
98  G4RotationMatrix cemc_rotm;
99  // put our cemc at center displacement in x
100  double radius_at_center = inner_radius + (outer_radius - inner_radius) / 2.;
101  double xcenter = radius_at_center * cos(cemc_angular_coverage / 2.);
102  double ycenter = radius_at_center * sin(cemc_angular_coverage / 2.);
103  cemc_rotm.rotateX(x_rot);
104  cemc_rotm.rotateY(y_rot);
105  cemc_rotm.rotateZ(z_rot);
106  new G4PVPlacement(G4Transform3D(cemc_rotm, G4ThreeVector(place_in_x - xcenter, place_in_y - ycenter, place_in_z)), cemc_log, "CEmc", logicWorld, false, false, OverlapCheck());
107  G4VSolid* tower_tub = new G4Tubs("TowerTub", inner_radius - no_overlap, outer_radius + no_overlap, (w_dimension[2] + no_overlap) / 2., 0, tower_angular_coverage);
108  G4LogicalVolume* tower_log = new G4LogicalVolume(tower_tub, Air, G4String("CEmcTower"), nullptr, nullptr, nullptr);
109  // G4VisAttributes* towerVisAtt = new G4VisAttributes();
110  // towerVisAtt->SetVisibility(true);
111  // towerVisAtt->SetForceSolid(true);
112  // towerVisAtt->SetColour(G4Colour::Blue());
113  // tower_log->SetVisAttributes(towerVisAtt);
114 
115  ostringstream tower_vol_name;
116  for (int i = 0; i < 7; i++)
117  {
118  tower_vol_name << "CEmcTower_" << i;
119  double phi = -i * tower_angular_coverage;
120  G4RotationMatrix* tower_rotm = new G4RotationMatrix();
121  tower_rotm->rotateZ(phi * rad);
122  new G4PVPlacement(tower_rotm, G4ThreeVector(0, 0, 0), tower_log, tower_vol_name.str(), cemc_log, false, i, OverlapCheck());
123  tower_vol_name.str("");
124  }
125  ConstructTowerVolume(tower_log);
126  return;
127 }
128 
129 // here we build up the tower from the sandwichs (tungsten + scintillator)
130 int PHG4CEmcTestBeamDetector::ConstructTowerVolume(G4LogicalVolume* tower_log)
131 {
133  G4Material* Air = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
134  G4VSolid* sandwich_box = new G4Box("Sandwich_box",
135  w_dimension[0] / 2., sandwich_thickness / 2., w_dimension[2] / 2.);
136 
137  G4LogicalVolume* sandwich_log = new G4LogicalVolume(sandwich_box,
138  Air,
139  G4String("CEmcSandwich"),
140  nullptr, nullptr, nullptr);
141  G4VisAttributes* sandwichVisAtt = new G4VisAttributes();
142  sandwichVisAtt->SetVisibility(true);
143  sandwichVisAtt->SetForceSolid(true);
144  sandwichVisAtt->SetColour(G4Colour::White());
145  sandwich_log->SetVisAttributes(sandwichVisAtt);
146  ostringstream sandwich_name;
147  for (int i = 0; i < 12; i++)
148  {
149  G4RotationMatrix* sandwich_rotm = new G4RotationMatrix();
150  double phi = -i * alpha;
151  sandwich_rotm->rotateZ(phi * rad);
152  sandwich_name << "CEmcSandwich_" << i;
153  double xshift = cos(phi) * (inner_radius + (outer_radius - inner_radius) / 2.);
154  double yshift = -sin(phi) * (inner_radius + (outer_radius - inner_radius) / 2.);
155  // we need to shift everything up by sandwich_thickness/2, calculate the shift in x
156  // if we push the sandwich up by sandwich_thickness/2
157  double xcorr = atan(phi) * sandwich_thickness / 2.;
158  double ycorr = sandwich_thickness / 2.;
159 
160  new G4PVPlacement(sandwich_rotm, G4ThreeVector(xshift + xcorr, yshift + ycorr, 0),
161  sandwich_log,
162  sandwich_name.str(),
163  tower_log, false, i, OverlapCheck());
164  sandwich_name.str("");
165  }
166  ConstructSandwichVolume(sandwich_log); // put W and scinti into sandwich
167  return 0;
168 }
169 
170 // here we put a single tungsten + scintillator sandwich together
172 {
173  vector<G4LogicalVolume*> block_logic;
174  G4Material* AbsorberMaterial = GetDetectorMaterial("G4_W");
175  G4Material* ScintiMaterial = GetDetectorMaterial("G4_POLYSTYRENE");
176 
178  {
179  cout << "invalid active scintillator fraction " << active_scinti_fraction
180  << " try between 0 and 1" << endl;
181  }
182 
183  double sc_active_thickness = sc_dimension[1] * active_scinti_fraction;
184  double sc_passive_thickness = sc_dimension[1] - sc_active_thickness;
185 
186  G4VSolid* block_w = new G4Box("Tungsten_box",
187  w_dimension[0] / 2., w_dimension[1] / 2., w_dimension[2] / 2.);
188  G4VSolid* block_sc = new G4Box("Scinti_box",
189  sc_dimension[0] / 2., sc_active_thickness / 2., sc_dimension[2] / 2.);
190  G4VSolid* block_passive_sc = nullptr;
191  block_logic.push_back(new G4LogicalVolume(block_w,
192  AbsorberMaterial,
193  "Plate_log_W",
194  nullptr, nullptr, nullptr));
195  block_logic.push_back(new G4LogicalVolume(block_sc,
196  ScintiMaterial,
197  "Plate_log_Sc",
198  nullptr, nullptr, nullptr));
199  G4VisAttributes* matVis = new G4VisAttributes();
200  G4VisAttributes* matVis1 = new G4VisAttributes();
201  matVis->SetVisibility(true);
202  matVis->SetForceSolid(true);
203  matVis->SetColour(G4Colour::Red());
204  matVis1->SetVisibility(true);
205  matVis1->SetForceSolid(true);
206  matVis1->SetColour(G4Colour::Green());
207  block_logic[0]->SetVisAttributes(matVis);
208  block_logic[1]->SetVisAttributes(matVis1);
209 
210  // here is the math to set the positions of those tiles inside the sandwich. The center of the coordinate system is in the center
211  // of the sandwich (w -> thickness of tungsten layer, sc_a -> thickness of actvie scinti, sc_p -> thickness of passive scinti)
212  // distance to the bottom of the sandwich:
213  // -(w + sc_a + sc_p + w)/2 = -w/2 -sc_a/2 -sc_p/2 -w/2
214 
215  // implement the absorber at the bottom and at the top of the sandwich
216  // moving up by w/2 for the lowest layer tungsten:
217  // -w/2 -sc_a/2 -sc_p/2 -w/2 + w/2 = -w/2 -sc_a/2 -sc_p/2 = -(w+sc)/2
218  // where sc_a+sc_p = sc = scintillator thickness (sc_dimension[1])
219 
220  sandwich_vol[0] = new G4PVPlacement(nullptr, G4ThreeVector(0, -(w_dimension[1] + sc_dimension[1]) / 2., 0),
221  block_logic[0],
222  "CEmc_W_plate_down",
223  sandwich, false, 0, OverlapCheck());
224 
225  // top of the sandwich - starting at the bottom and add the tungsten and scintillator layer and half tungsten
226  // -w/2 -sc/2 -w/2 + w + sc + w/2 = +sc/2 +w/2 = (w+sc)/2
227  sandwich_vol[1] = new G4PVPlacement(nullptr, G4ThreeVector(0, (w_dimension[1] + sc_dimension[1]) / 2., 0),
228  block_logic[0],
229  "CEmc_W_plate_up",
230  sandwich, false, 0, OverlapCheck());
231 
232  // since we split the scintillator into an active and passive part to accomodate
233  // for the scintillator fibers not occupying 100% of the volume.
234  // this is mocked up by 2 separate layers of plastic, active scintillator is the middle layer
235  // again going to the bottom of the sandwich box (sc_a -> active sc, sc_p -> passive sc)
236  // -w/2 -sc_a/2 -sc_p/2 -w/2
237  // -(w + sc_a + sc_p + w)/2, adding the w layer and half of the sc_a to get to the middle of the sc_a layer:
238  // -w/2 -sc_a/2 - sc_p/2 -w/2 + w + sc_a/2 = -sc_p/2
239  // if fraction of active sc is 1, sc_passive_thickness is zero
240  sandwich_vol[2] = new G4PVPlacement(nullptr, G4ThreeVector(0, -sc_passive_thickness / 2., 0),
241  block_logic[1],
242  "CEmc_Scinti_plate",
243  sandwich, false, 0, OverlapCheck());
244 
245  if (sc_passive_thickness > 0)
246  {
247  G4VisAttributes* matVis2 = new G4VisAttributes();
248  matVis1->SetVisibility(true);
249  matVis1->SetForceSolid(true);
250  matVis1->SetColour(G4Colour::Blue());
251  block_passive_sc = new G4Box("Passive_Scinti_box",
252  sc_dimension[0] / 2., sc_passive_thickness / 2., sc_dimension[2] / 2.);
253  block_logic.push_back(new G4LogicalVolume(block_passive_sc,
254  ScintiMaterial,
255  "Plate_log_Passive_Sc",
256  nullptr, nullptr, nullptr));
257  block_logic[2]->SetVisAttributes(matVis2);
258  // here we go again - bottome of sandwich box is
259  // -(w + sc_a + sc_p +w)/2, now add w and sc_a and half of sc_p:
260  // -w/2 -sc_a/2 - sc_p/2 -w/2 + w + sc_a + sc_p/2 = sc_a/2
261  sandwich_vol[3] = new G4PVPlacement(nullptr, G4ThreeVector(0, sc_active_thickness / 2., 0),
262  block_logic[2],
263  "CEmc_Passive_Si_plate",
264  sandwich, false, 0, OverlapCheck());
265  }
266  else
267  {
268  sandwich_vol[3] = nullptr;
269  }
270  return 0;
271 }
272 
274 {
275  // calculate the inner/outer radius of the g4tub using the test setup numbers
276  // 1mm scintillator, 0.5mm tungsten, 0.25 mm gap at the back of the detector
277  // the geometry with figure is explained in our wiki:
278  // https://www.phenix.bnl.gov/WWW/offline/wikioff/index.php/CEmc
279  // get the alpha angle via law of cos (a is the 0.25 mm gap):
280  // a^2 = b^2+c^2 - 2bc*cos(alpha)
281  // cos(alpha) = (b^2+c^2 - a^2)/2bc
282  // with b=c
283  // cos(alpha) = 1-(a^2/2b^2)
284  alpha = acos(1 - (gap * gap) / (2 * plate_x * plate_x));
285  // inner radius
287  outer_radius = sqrt((inner_radius + plate_x) * (inner_radius + plate_x) + sandwich_thickness * sandwich_thickness);
288 
289  // twice no_overlap to apply safety margin also for towers
290  // angular coverage of 1 tower is 12 sandwiches = sandwiches_per_tower*alpha
291  tower_angular_coverage = sandwiches_per_tower * alpha + 0.00005 * deg; // safety margin added to remove 142 nm overlap
292  // cemc prototype has 7 towers
294  return;
295 }