Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4HcalPrototypeDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4HcalPrototypeDetector.cc
1 // This is G4 detector implementation for the hcal prototype
2 // Created on 1/27/2014, Liang, HeXC
3 // Updated on 3/21/2014, Liang, HeXC
4 // Updated on 4/18/2014, Liang, HeXC, Updating detector construction (proper tilting angles!)
5 
6 #include "PHG4HcalPrototypeDetector.h"
7 
8 #include "PHG4HcalPrototypeDetectorMessenger.h" // for PHG4HcalPrototypeDet...
9 
10 #include <g4main/PHG4Detector.h> // for PHG4Detector
11 
12 #include <Geant4/G4Box.hh>
13 #include <Geant4/G4Colour.hh>
14 #include <Geant4/G4LogicalVolume.hh>
15 #include <Geant4/G4Material.hh>
16 #include <Geant4/G4NistManager.hh>
17 #include <Geant4/G4PVPlacement.hh>
18 #include <Geant4/G4RotationMatrix.hh> // for G4RotationMatrix
19 #include <Geant4/G4RunManager.hh>
20 #include <Geant4/G4String.hh> // for G4String
21 #include <Geant4/G4SystemOfUnits.hh> // for cm, mm, deg, rad
22 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
23 #include <Geant4/G4Transform3D.hh> // for G4Transform3D
24 #include <Geant4/G4Trap.hh>
25 #include <Geant4/G4Types.hh> // for G4double, G4int, G4bool
26 #include <Geant4/G4UnionSolid.hh>
27 #include <Geant4/G4VisAttributes.hh>
28 #include <Geant4/G4ios.hh> // for G4cout, G4endl
29 
30 #include <cmath> // for cos, sin, M_PI, asin
31 #include <sstream>
32 
33 class PHCompositeNode;
34 
35 using namespace std;
36 
37 // static double no_overlap = 0.00015 * cm; // added safety margin against overlaps by using same boundary between volumes
38 
39 PHG4HcalPrototypeDetector::PHG4HcalPrototypeDetector(PHG4Subsystem* subsys, PHCompositeNode* Node, const std::string& dnam, const int lyr)
40  : PHG4Detector(subsys, Node, dnam)
41  , nScint360(15)
42  , // This is legacy variable in case one wants to build a cylindrical calorimeter
43  nHcal1Layers(16)
44  , nHcal2Layers(nHcal1Layers)
45  , hcal2ScintSizeX(2.54 * 26.1 * cm)
46  , // from Don's drawing
47  hcal2ScintSizeY(0.85 * cm)
48  , // Changed from 8.0cm to 8.5cm following
49  hcal2ScintSizeZ(2.54 * 29.53 * cm)
50  , // from Don's drawing
51  hcal1ScintSizeX(2.54 * 12.43 * cm)
52  , // from Don's drawing
53  hcal1ScintSizeY(0.85 * cm)
54  , // Changed from 8.0cm to 8.5cm following
55  hcal1ScintSizeZ(2.54 * 17.1 * cm)
56  , // from Don's drawing
57  hcal1TiltAngle(0.27)
58  , hcal2TiltAngle(0.14)
59  , hcal1DPhi(0.27)
60  , // calculated based on Don's drawing. (twopi/16.0);
61  hcal2DPhi(0.288)
62  , // calculated based on Don's drawing. (twopi/16.0);
63  hcal1RadiusIn(1855 * mm)
64  , hcal2RadiusIn(hcal1RadiusIn + hcal1ScintSizeX)
65  ,
66  // The frame box for the HCAL.
67  hcalBoxSizeX(1.1 * (hcal2ScintSizeX + hcal1ScintSizeX))
68  , hcalBoxSizeY(2.54 * 40.0 * cm)
69  , // need to get more accurate dimension for the box
70  hcalBoxSizeZ(1.1 * hcal2ScintSizeZ)
71  , hcalBoxRotationAngle_z(0.0 * rad)
72  , // Rotation along z-axis
73  hcalBoxRotationAngle_y(0.0 * rad)
74  , // Rotation along y-axis
75  hcal2Abs_dxa(2.54 * 27.1 * cm)
76  , // these numbers are from Don't drawings
77  hcal2Abs_dxb(hcal2Abs_dxa)
78  , hcal2Abs_dya(2.54 * 1.099 * cm)
79  , hcal2Abs_dyb(2.54 * 1.776 * cm)
80  , hcal2Abs_dz(hcal2ScintSizeZ)
81  , hcal1Abs_dxa(hcal1ScintSizeX)
82  , // these numbers are from Don't drawings
83  hcal1Abs_dxb(hcal1ScintSizeX)
84  , hcal1Abs_dya(2.54 * 0.787 * cm)
85  , hcal1Abs_dyb(2.54 * 1.115 * cm)
86  , hcal1Abs_dz(hcal1ScintSizeZ)
87  , hcalJunctionSizeX(2.54 * 1.0 * cm)
88  , hcalJunctionSizeY(hcal2ScintSizeY)
89  , hcalJunctionSizeZ(hcal2Abs_dz)
90  , physiWorld(nullptr)
91  , logicWorld(nullptr)
92  , logicHcalBox(nullptr)
93  , solidHcalBox(nullptr)
94  , physiHcalBox(nullptr)
95  , logicHcal2ScintLayer(nullptr)
96  , logicHcal1ScintLayer(nullptr)
97  , logicHcal2AbsLayer(nullptr)
98  , logicHcal1AbsLayer(nullptr)
99  , world_mat(nullptr)
100  , steel(nullptr)
101  , scint_mat(nullptr)
102  , active(0)
103  , absorberactive(0)
104  , layer(lyr)
105  , blackhole(0)
106 {
107  // create commands for interactive definition of the detector
108  fDetectorMessenger = new PHG4HcalPrototypeDetectorMessenger(this);
109 }
110 
111 //_______________________________________________________________
112 //_______________________________________________________________
113 int PHG4HcalPrototypeDetector::IsInHcalPrototype(G4VPhysicalVolume* /*volume*/) const
114 {
115  return 0; // not sure what value to return for now.
116 }
117 
118 void PHG4HcalPrototypeDetector::ConstructMe(G4LogicalVolume* world)
119 {
120  logicWorld = world;
121 
122  DefineMaterials();
123 
124  ConstructDetector();
125 }
126 
127 void PHG4HcalPrototypeDetector::DefineMaterials()
128 {
129  // Water is defined from NIST material database
130  G4NistManager* nist = G4NistManager::Instance();
131 
132  world_mat = nist->FindOrBuildMaterial("G4_AIR");
133 
134  steel = nist->FindOrBuildMaterial("G4_Fe"); // Need to fix this *******
135  scint_mat = nist->FindOrBuildMaterial("G4_POLYSTYRENE");
136 
137  G4cout << *(G4Material::GetMaterialTable()) << G4endl;
138 }
139 
140 // We now build our detector solids, etc
141 //
142 G4VPhysicalVolume* PHG4HcalPrototypeDetector::ConstructDetector()
143 {
144  // Option to switch on/off checking of volumes overlaps
145  //
146  G4bool checkOverlaps = true;
147 
148  // HCAL Frame Box for enclosing the inner and the outer hcal
149  // which allows us to rotate the whole calorimeter easily
150  solidHcalBox = new G4Box("HcalBox", //its name
151  hcalBoxSizeX / 2, hcalBoxSizeY / 2, hcalBoxSizeZ / 2); //its size
152 
153  logicHcalBox = new G4LogicalVolume(solidHcalBox, //its solid
154  world_mat, //its material
155  "HcalBox"); //its name
156 
157  // Place the HcalBox inside the World
158  G4ThreeVector boxVector = G4ThreeVector(hcal1RadiusIn + hcalBoxSizeX / 2.0, 0, 0);
159  G4RotationMatrix rotBox = G4RotationMatrix();
160 
161  rotBox.rotateZ(hcalBoxRotationAngle_z); // Rotate the whole calorimeter box along z-axis
162  rotBox.rotateY(hcalBoxRotationAngle_y); // Rotate the whole calorimeter box along z-axis
163 
164  G4Transform3D boxTransform = G4Transform3D(rotBox, boxVector);
165 
166  physiHcalBox = new G4PVPlacement(boxTransform, // rotation + positioning
167  logicHcalBox, //its logical volume
168  "HcalBox", //its name
169  logicWorld, //its mother volume
170  false, //no boolean operation
171  0, //copy number
172  checkOverlaps); //checking overlaps
173 
174  // HCal Box visualization attributes
175  G4VisAttributes* hcalBoxVisAtt = new G4VisAttributes(G4Colour(0.0, 0.0, 1.0)); // blue
176  hcalBoxVisAtt->SetVisibility(true);
177  //hcalBoxVisAtt->SetForceSolid(true);
178  logicHcalBox->SetVisAttributes(hcalBoxVisAtt);
179 
180  // Make outer hcal section
181  //
182  // Build scintillator layer box as a place holder for the scintillator sheets
183  // Make it 0.05cm thinner than the true gap width for avoiding geometry overlaps
184  G4Box* hcal2ScintLayer = new G4Box("hcal2ScintLayer", //its name
185  hcal2ScintSizeX / 2, (hcal2ScintSizeY - 0.05 * cm) / 2, hcal2ScintSizeZ / 2); //its size
186 
187  logicHcal2ScintLayer =
188  new G4LogicalVolume(hcal2ScintLayer, // its solid name
189  world_mat, // material, the same as the material of the world valume
190  "hcal2ScintLayer"); // its name
191 
192  // Scintillator visualization attributes
193  G4VisAttributes* scintVisAtt = new G4VisAttributes(G4Colour(1.0, 0.0, 0.0)); //Red
194  scintVisAtt->SetVisibility(true);
195  //scintVisAtt->SetForceSolid(true);
196  logicHcal2ScintLayer->SetVisAttributes(scintVisAtt);
197 
198  // Construct outer scintillator 1U
199  G4double outer1UpDz = 649.8 * mm;
200  G4double outer1UpDy1 = 2 * 0.35 * cm;
201  G4double outer1UpDx2 = 179.3 * mm;
202  G4double outer1UpDx4 = 113.2 * mm;
203 
204  G4Trap* outer1USheetSolid = new G4Trap("outer1USheet", //its name
205  outer1UpDy1, outer1UpDz,
206  outer1UpDx2, outer1UpDx4); //its size
207 
208  G4LogicalVolume* logicOuter1USheet = new G4LogicalVolume(outer1USheetSolid,
209  scint_mat,
210  "outer1USheet");
211 
212  // Scintillator visualization attributes
213  G4VisAttributes* scintSheetVisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 0.0)); //White
214  scintSheetVisAtt->SetVisibility(true);
215  scintSheetVisAtt->SetForceSolid(true);
216 
217  logicOuter1USheet->SetVisAttributes(scintSheetVisAtt);
218 
219  // Detector placement
220  // Position the outer right 1U sheet
221  G4ThreeVector threeVecOuter1U_1 = G4ThreeVector(0 * cm, 0 * cm, -0.252 * (outer1UpDx2 + outer1UpDx4));
222  G4RotationMatrix rot1U_1 = G4RotationMatrix();
223  rot1U_1.rotateZ(90 * deg);
224  rot1U_1.rotateX(-90 * deg);
225 
226  G4Transform3D transformOuter1U_1 = G4Transform3D(rot1U_1, threeVecOuter1U_1);
227 
228  new G4PVPlacement(transformOuter1U_1,
229  logicOuter1USheet,
230  "outer1USheet",
231  logicHcal2ScintLayer,
232  false,
233  0, // Copy one
234  checkOverlaps);
235 
236  // Detector placement
237  // Position the outer left 1U sheet
238  G4ThreeVector threeVecOuter1U_2 = G4ThreeVector(0 * cm, 0 * cm, 0.252 * (outer1UpDx2 + outer1UpDx4));
239  G4RotationMatrix rot1U_2 = G4RotationMatrix();
240  rot1U_2.rotateZ(90 * deg);
241  rot1U_2.rotateX(90 * deg);
242 
243  G4Transform3D transformOuter1U_2 = G4Transform3D(rot1U_2, threeVecOuter1U_2);
244 
245  new G4PVPlacement(transformOuter1U_2,
246  logicOuter1USheet,
247  "outer1USheet",
248  logicHcal2ScintLayer,
249  false,
250  1, // Copy two
251  checkOverlaps);
252 
253  // Construct outer scintillator 2U
254  G4double outer2UpDz = 0.5 * 649.8 * mm;
255  G4double outer2UpTheta = 8.8 * M_PI / 180.;
256  G4double outer2UpPhi = 0.0 * M_PI / 180.;
257  G4double outer2UpDy1 = 0.35 * cm;
258  G4double outer2UpDy2 = 0.35 * cm;
259  G4double outer2UpDx1 = 0.5 * 179.3 * mm, outer2UpDx2 = 0.5 * 179.3 * mm;
260  G4double outer2UpDx3 = 0.5 * 113.2 * mm, outer2UpDx4 = 0.5 * 113.2 * mm;
261  G4double outer2UpAlp1 = 0. * M_PI / 180., outer2UpAlp2 = 0. * M_PI / 180;
262 
263  G4Trap* outer2USheetSolid = new G4Trap("outer2USheet",
264  outer2UpDz,
265  outer2UpTheta,
266  outer2UpPhi,
267  outer2UpDy1,
268  outer2UpDx1,
269  outer2UpDx2,
270  outer2UpAlp1,
271  outer2UpDy2,
272  outer2UpDx3,
273  outer2UpDx4,
274  outer2UpAlp2);
275 
276  G4LogicalVolume* logicOuter2USheet = new G4LogicalVolume(outer2USheetSolid,
277  scint_mat,
278  "outer2USheet");
279 
280  logicOuter2USheet->SetVisAttributes(scintSheetVisAtt);
281 
282  // Detector placement
283  // Position the right most sheet
284  G4ThreeVector threeVecOuter2U_1 = G4ThreeVector(0 * cm, 0 * cm, -0.755 * (outer1UpDx2 + outer1UpDx4));
285  G4RotationMatrix rot2U_1 = G4RotationMatrix();
286  rot2U_1.rotateY(-90 * deg);
287 
288  G4Transform3D transformOuter2U_1 = G4Transform3D(rot2U_1, threeVecOuter2U_1);
289 
290  new G4PVPlacement(transformOuter2U_1,
291  logicOuter2USheet,
292  "outer2USheet",
293  logicHcal2ScintLayer,
294  false,
295  0, // copy one
296  checkOverlaps);
297 
298  // Detector placement
299  // Position the outer left most sheet
300  G4ThreeVector threeVecOuter2U_2 = G4ThreeVector(0 * cm, 0 * cm, 0.755 * (outer1UpDx2 + outer1UpDx4));
301  G4RotationMatrix rot2U_2 = G4RotationMatrix();
302  rot2U_2.rotateY(-90 * deg);
303  rot2U_2.rotateX(180 * deg);
304 
305  G4Transform3D transformOuter2U_2 = G4Transform3D(rot2U_2, threeVecOuter2U_2);
306 
307  new G4PVPlacement(transformOuter2U_2,
308  logicOuter2USheet,
309  "outer2USheet",
310  logicHcal2ScintLayer,
311  false,
312  1, // copy two
313  checkOverlaps);
314 
315  G4Trap* hcal2AbsLayer =
316  new G4Trap("hcal2AbsLayer", //its name
317  hcal2Abs_dz, hcal2Abs_dxa,
318  hcal2Abs_dyb, hcal2Abs_dya); //its size
319 
320  // Add extra absorber materials at the junction
321  G4Box* hcalJunction = new G4Box("HcalJunction", //its name
322  hcalJunctionSizeX / 2, hcalJunctionSizeY / 2, hcalJunctionSizeZ / 2); //its size
323 
324  // define the transformation for placing the junction to the inner side of the hcal2 absorber layer
325  G4ThreeVector threeVecJunction = G4ThreeVector(-hcal2Abs_dya / 2.0 - 1.1 * hcalJunctionSizeY, hcal2Abs_dxa / 2.0 - hcalJunctionSizeX / 2, 0.0 * cm);
326  G4RotationMatrix rotJunction = G4RotationMatrix();
327  rotJunction.rotateZ(90 * deg);
328  rotJunction.rotateX(0 * deg);
329  G4Transform3D transformJunction = G4Transform3D(rotJunction, threeVecJunction);
330 
331  // attach the junction to the absorber using G4UnionSolid
332  G4UnionSolid* hcal2AbsLayerJunct = new G4UnionSolid("hcal2AbsLayerJunct", hcal2AbsLayer, hcalJunction, transformJunction);
333 
334  logicHcal2AbsLayer = new G4LogicalVolume(hcal2AbsLayerJunct, // hcal2AbsLayer,
335  steel,
336  "hcal2AbsLayer");
337 
338  // hcal2Aborber visualization attributes
339  G4VisAttributes* hcal2AbsVisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 1.0));
340  hcal2AbsVisAtt->SetVisibility(true);
341  logicHcal2AbsLayer->SetVisAttributes(hcal2AbsVisAtt);
342 
343  // place scintillator layers insides the HCal2 absorber volume
344  //
345 
346  G4double theta = 0.;
347  G4double theta2 = hcal2TiltAngle; // I am not convinced that I got the tilt angle right. So I simply forced it to
348  G4double rScintLayerCenter = hcal2RadiusIn + hcal2ScintSizeX / 2.0 + 2.54 * 1.0 * cm; // move the scintillator toward the rear end of HCAL2
349  G4double RmidScintLayerX = 0.81 * (hcal2ScintSizeX - hcal1ScintSizeX) / 2.0 + 2.54 * 1.0 * cm + 2 * 2.54 * cm; // move the scintillator toward the rear end of HCAL2
350  G4double rAbsLayerCenter = hcal2RadiusIn + hcal2Abs_dxa / 2.0;
351  G4double RmidAbsLayerX = 0.81 * (hcal2Abs_dxa - hcal1ScintSizeX) / 2.0 + 2 * 2.54 * cm;
352  G4int LayerNum = 1;
353  G4double xPadding = 0.0;
354  G4double yPadding = 0.0;
355  G4double tiltPadding = 0.0;
356 
357  for (G4int iLayer = -nHcal2Layers / 2; iLayer < nHcal2Layers / 2; iLayer++)
358  {
359  //
360  // Place outer absorber layer first
361  //
362  theta = hcal2DPhi / nHcal2Layers * iLayer;
363  G4double xposShift = rAbsLayerCenter * (cos(theta) - 1.0);
364  G4double yposShift = rAbsLayerCenter * sin(theta);
365  xPadding = 0.013 * RmidAbsLayerX * LayerNum; // adding an extra padding for placing the absorber and scintillator
366  G4ThreeVector absTrans = G4ThreeVector(RmidAbsLayerX + xposShift - xPadding, yposShift, 0);
367  G4RotationMatrix rotAbsLayer = G4RotationMatrix();
368  rotAbsLayer.rotateY(90 * deg);
369  rotAbsLayer.rotateX(90 * deg);
370  rotAbsLayer.rotateY(-90 * deg);
371  tiltPadding = LayerNum * 0.005 * rad;
372  rotAbsLayer.rotateZ((iLayer + nHcal2Layers / 2) * 0.02 * rad + tiltPadding); // Added a funge increment factor 0.01
373 
374  G4Transform3D transformAbs = G4Transform3D(rotAbsLayer, absTrans);
375 
376  new G4PVPlacement(transformAbs, //rotation,position
377  logicHcal2AbsLayer, //its logical volume
378  "hcal2AbsLayer", //its name
379  logicHcalBox, //its mother volume
380  false, //no boolean operation
381  LayerNum - 1, //copy number
382  checkOverlaps); //overlaps checking
383  //
384  // Place outer hcal scintillator layer
385  //
386  theta += 0.5 * hcal2DPhi / nHcal2Layers;
387  G4cout << "M_PI: " << M_PI << " theta: " << theta << " TileAngle: " << theta2 << G4endl;
388  xposShift = rScintLayerCenter * (cos(theta) - 1.0);
389  yposShift = rScintLayerCenter * sin(theta);
390  G4ThreeVector myTrans = G4ThreeVector(RmidScintLayerX + xposShift - xPadding, yposShift + 0.3 * cm, 0); // hcal2RadiusIn + hcal1ScintSizeX/2.0
391  G4RotationMatrix rotm = G4RotationMatrix();
392  rotm.rotateZ((iLayer + nHcal2Layers / 2 + 1) * 0.020 * rad + tiltPadding); // Added a funge increment factor 0.02
393  G4Transform3D transform = G4Transform3D(rotm, myTrans);
394 
395  G4cout << " iLayer " << iLayer << G4endl;
396 
397  new G4PVPlacement(transform, //rotation,position
398  logicHcal2ScintLayer, //its logical volume
399  "hcal2ScintLayer", //its name
400  logicHcalBox, //its mother volume
401  false, //no boolean operation
402  LayerNum - 1, //copy number
403  checkOverlaps); //checking overlaps
404  LayerNum++;
405  }
406 
407  // Make INNER hcal section
408  //
409  // Build scintillator layer box as a place holder for the scintillator sheets
410  // Make it 0.05cm thinner than the true gap width for avoiding geometry overlaps
411  G4Box* hcal1ScintLayer = new G4Box("hcal1ScintLayer", //its name
412  hcal1ScintSizeX / 2, (hcal1ScintSizeY - 0.05 * cm) / 2, hcal1ScintSizeZ / 2); //its size
413 
414  logicHcal1ScintLayer =
415  new G4LogicalVolume(hcal1ScintLayer, // its solid name
416  world_mat, // simply use world material
417  "hcal1ScintLayer"); // its name
418 
419  logicHcal1ScintLayer->SetVisAttributes(scintVisAtt);
420 
421  // Constructing scintillator sheets
422  // Construct inner scintillator 1U
423  // The numbers are read off from Don's drawings
424  G4double inner1UpDz = 316.8 * mm;
425  G4double inner1UpDy1 = 2 * 0.35 * cm;
426  G4double inner1UpDx2 = 108.6 * mm;
427  G4double inner1UpDx4 = 77.4 * mm;
428 
429  G4Trap* inner1USheetSolid = new G4Trap("inner1USheet", //its name
430  inner1UpDy1, inner1UpDz,
431  inner1UpDx2, inner1UpDx4); //its size
432 
433  G4LogicalVolume* logicInner1USheet = new G4LogicalVolume(inner1USheetSolid,
434  scint_mat,
435  "inner1USheet");
436  logicInner1USheet->SetVisAttributes(scintSheetVisAtt);
437 
438  // Detector placement
439  // Position the inner right 1U sheet
440  G4ThreeVector threeVecInner1U_1_inner = G4ThreeVector(0 * cm, 0 * cm, -0.252 * (inner1UpDx2 + inner1UpDx4));
441 
442  G4Transform3D transformInner1U_1 = G4Transform3D(rot1U_1, threeVecInner1U_1_inner);
443 
444  new G4PVPlacement(transformInner1U_1,
445  logicInner1USheet,
446  "inner1USheet",
447  logicHcal1ScintLayer,
448  false,
449  0, // Copy one
450  checkOverlaps);
451 
452  // Detector placement
453  // Position the inner left 1U sheet
454  G4ThreeVector threeVecInner1U_2_inner = G4ThreeVector(0 * cm, 0 * cm, 0.252 * (inner1UpDx2 + inner1UpDx4));
455  //G4RotationMatrix rot1U_2 = G4RotationMatrix();
456  //rot1U_2.rotateZ(90*deg);
457  //rot1U_2.rotateX(90*deg);
458 
459  G4Transform3D transformInner1U_2 = G4Transform3D(rot1U_2, threeVecInner1U_2_inner);
460 
461  new G4PVPlacement(transformInner1U_2,
462  logicInner1USheet,
463  "inner1USheet",
464  //logicWorld,
465  logicHcal1ScintLayer,
466  false,
467  1, // Copy two
468  checkOverlaps);
469 
470  // Construct inner scintillator 2U
471  G4double inner2UpDz = 0.5 * 316.8 * mm;
472  G4double inner2UpTheta = 8.8 * M_PI / 180.;
473  G4double inner2UpPhi = 0.0 * M_PI / 180.;
474  G4double inner2UpDy1 = 0.35 * cm;
475  G4double inner2UpDy2 = 0.35 * cm;
476  G4double inner2UpDx1 = 0.5 * 108.6 * mm, inner2UpDx2 = 0.5 * 108.6 * mm;
477  G4double inner2UpDx3 = 0.5 * 77.4 * mm, inner2UpDx4 = 0.5 * 77.4 * mm;
478  G4double inner2UpAlp1 = 0. * M_PI / 180., inner2UpAlp2 = 0. * M_PI / 180;
479 
480  G4Trap* inner2USheetSolid = new G4Trap("inner2USheet",
481  inner2UpDz,
482  inner2UpTheta,
483  inner2UpPhi,
484  inner2UpDy1,
485  inner2UpDx1,
486  inner2UpDx2,
487  inner2UpAlp1,
488  inner2UpDy2,
489  inner2UpDx3,
490  inner2UpDx4,
491  inner2UpAlp2);
492 
493  G4LogicalVolume* logicInner2USheet = new G4LogicalVolume(inner2USheetSolid,
494  scint_mat,
495  "inner2USheet");
496 
497  logicInner2USheet->SetVisAttributes(scintSheetVisAtt);
498 
499  // Detector placement
500  // Position the right most sheet
501  G4ThreeVector threeVecInner2U_1_inner = G4ThreeVector(0 * cm, 0 * cm, -0.755 * (inner1UpDx2 + inner1UpDx4));
502  //G4RotationMatrix rot2U_1 = G4RotationMatrix();
503  //rot2U_1.rotateY(-90*deg);
504 
505  G4Transform3D transformInner2U_1 = G4Transform3D(rot2U_1, threeVecInner2U_1_inner);
506 
507  new G4PVPlacement(transformInner2U_1,
508  logicInner2USheet,
509  "inner2USheet",
510  //logicWorld,
511  logicHcal1ScintLayer,
512  false,
513  0, // copy one
514  checkOverlaps);
515 
516  // Detector placement
517  // Position the inner left most sheet
518  G4ThreeVector threeVecInner2U_2_inner = G4ThreeVector(0 * cm, 0 * cm, 0.755 * (inner1UpDx2 + inner1UpDx4));
519  //G4RotationMatrix rot2U_2 = G4RotationMatrix();
520  //rot2U_2.rotateY(-90*deg);
521  //rot2U_2.rotateX(180*deg);
522 
523  G4Transform3D transformInner2U_2 = G4Transform3D(rot2U_2, threeVecInner2U_2_inner);
524 
525  new G4PVPlacement(transformInner2U_2,
526  logicInner2USheet,
527  "inner2USheet",
528  //logicWorld,
529  logicHcal1ScintLayer,
530  false,
531  1, // copy two
532  checkOverlaps);
533 
534  // Build hcal1 absorber layer in trapezoid shape
535  //
536  /*
537  G4double hcal1Abs_dxa = hcal1ScintSizeX; // hcal1Abs_dxb = hcal1ScintSizeX;
538  G4double hcal1Abs_dya = 0.787*2.54*cm, hcal1Abs_dyb = 1.115*2.54*cm;
539  // these numbers are from Don't drawings
540  G4double hcal1Abs_dz = hcal1ScintSizeZ;
541  */
542  G4Trap* hcal1AbsLayer =
543  new G4Trap("hcal1AbsLayer", //its name
544  hcal1Abs_dz, hcal1Abs_dxa,
545  hcal1Abs_dyb, hcal1Abs_dya); //its size
546 
547  logicHcal1AbsLayer = new G4LogicalVolume(hcal1AbsLayer,
548  steel,
549  "hcal1AbsLayer");
550 
551  logicHcal1AbsLayer->SetVisAttributes(hcal2AbsVisAtt);
552 
553  // place scintillator layers insides the HCal1 absorber volume
554  //
555  // Calculate the title angle
556 
557  theta2 = hcal1TiltAngle;
558  rAbsLayerCenter = hcal1RadiusIn + hcal1ScintSizeX / 2.0;
559  // RmidAbsLayerX = (hcal2Abs_dxa - hcal1ScintSizeX)/2.0;
560  RmidAbsLayerX = hcal2Abs_dxa / 2.0;
561  LayerNum = 1; // These three parameters are for making fine adjustment of the placement of inner hcal
562  yPadding = -2.54 * 0.8 * cm;
563  for (G4int iLayer = -nHcal2Layers / 2 - 1; iLayer < nHcal2Layers / 2 - 1; iLayer++)
564  { // shift one layer down
565  //
566  // Place inner absorber layer first
567  //
568  theta = hcal1DPhi / nHcal1Layers * iLayer; // add an off set 0.01 rad
569  G4double xposShift = rAbsLayerCenter * (cos(theta) - 1.0);
570  G4double yposShift = rAbsLayerCenter * sin(theta);
571  xPadding = 0.005 * RmidAbsLayerX * LayerNum - 2.54 * 1.9 * cm;
572  G4ThreeVector absTrans = G4ThreeVector(-RmidAbsLayerX * cos(theta) + xposShift - xPadding, yposShift + yPadding, 0);
573  G4RotationMatrix rotAbsLayer = G4RotationMatrix();
574  rotAbsLayer.rotateY(90 * deg);
575  rotAbsLayer.rotateX(90 * deg);
576  rotAbsLayer.rotateY(-90 * deg);
577  tiltPadding = LayerNum * 0.005 * rad;
578  rotAbsLayer.rotateZ((theta - theta2) * rad + tiltPadding);
579  //rotAbsLayer.rotateZ(-(iLayer + nHcal2Layers/2)*0.02*rad + tiltPadding);
580 
581  G4Transform3D transformAbs = G4Transform3D(rotAbsLayer, absTrans);
582 
583  new G4PVPlacement(transformAbs, //rotation,position
584  logicHcal1AbsLayer, //its logical volume
585  "hcal1AbsLayer", //its name
586  logicHcalBox, //its mother volume
587  false, //no boolean operation
588  LayerNum - 1, //copy number
589  checkOverlaps); //overlaps checking
590 
591  //
592  // Place inner hcal scintillator layer
593  //
594 
595  theta += 0.5 * hcal1DPhi / nHcal1Layers;
596  //G4cout << "M_PI: " << M_PI << " theta: " << theta << " TileAngle: " << theta2 << G4endl;
597  //G4double phi = iLayer*dPhi;
598  // G4ThreeVector myTrans = G4ThreeVector(-RmidAbsLayerX*cos(theta) + xposShift, yposShift, 0);
599  xposShift = rAbsLayerCenter * (cos(theta) - 1.0);
600  yposShift = rAbsLayerCenter * sin(theta);
601  G4ThreeVector myTrans = G4ThreeVector(-RmidAbsLayerX * cos(theta) + xposShift - xPadding, yposShift + yPadding, 0);
602  //G4cout << " x: " << Rmid*cos(theta) << " y: " << Rmid*sin(theta) << " 1.0*rad: " << 1.8*rad << G4endl;
603  G4RotationMatrix rotm = G4RotationMatrix();
604  rotm.rotateZ((theta - theta2 + 0.01) * rad + tiltPadding); // Added a funge increment factor 0.01
605  //rotm.rotateZ(-(iLayer + nHcal2Layers/2+1)*0.02*rad + tiltPadding); // Added a funge increment factor 0.01
606  G4Transform3D transform = G4Transform3D(rotm, myTrans);
607 
608  G4cout << " iLayer " << iLayer << G4endl;
609 
610  new G4PVPlacement(transform, //rotation,position
611  logicHcal1ScintLayer, //its logical volume
612  "hcal1ScintLayer", //its name
613  logicHcalBox, //its mother volume: logicHcal1Ab
614  false, //no boolean operation
615  LayerNum - 1, //copy number
616  checkOverlaps); //checking overlaps
617  LayerNum++;
618  }
619 
620  return physiWorld;
621 }
622 
623 void PHG4HcalPrototypeDetector::CalculateGeometry()
624 {
625  return;
626 }
627 
628 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
629 
630 void PHG4HcalPrototypeDetector::SetMaterial(G4String /*materialChoice*/)
631 {
632 }
633 
634 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
635 
636 void PHG4HcalPrototypeDetector::UpdateGeometry()
637 {
638  G4RunManager::GetRunManager()->PhysicsHasBeenModified();
639  G4RunManager::GetRunManager()->DefineWorldVolume(ConstructDetector());
640 }
641 
642 // The following code is copied from the sPHENIX G4 simulation
643 void PHG4HcalPrototypeDetector::SetTiltViaNcross(const int ncross)
644 {
645  G4double sign = 1;
646  if (ncross < 0)
647  {
648  sign = -1;
649  }
650  G4int ncr = fabs(ncross);
651 
652  // Determine title angle for the outer hcal
653  //
654  G4double cSide = hcal2RadiusIn + hcal2ScintSizeX / 2;
655  G4double bSide = hcal2RadiusIn + hcal2ScintSizeX;
656 
657  G4double alpha = 0;
658  if (ncr > 1)
659  {
660  alpha = (360. / nHcal2Layers * M_PI / 180.) * (ncr - 1) / 2.0;
661  }
662  else
663  {
664  alpha = (360. / nHcal2Layers * M_PI / 180.) / 2.;
665  }
666 
667  G4double sinbSide = sin(alpha) * bSide / (sqrt(bSide * bSide + cSide * cSide - 2 * bSide * cSide * cos(alpha)));
668  G4double beta = asin(sinbSide); // This is the slat angle
669 
670  hcal2TiltAngle = beta * sign;
671 
672  // Determine title angle for the inner hcal
673  //
674  cSide = hcal1RadiusIn + hcal1ScintSizeX / 2;
675  bSide = hcal1RadiusIn + hcal1ScintSizeX;
676 
677  sinbSide = sin(alpha) * bSide / (sqrt(bSide * bSide + cSide * cSide - 2.0 * bSide * cSide * cos(alpha)));
678  beta = asin(sinbSide); // This is the slat angle
679 
680  hcal1TiltAngle = beta * sign;
681 
682  G4cout << " alpha : " << alpha << G4endl;
683  G4cout << " SetTitlViaNCross(" << ncross << ") setting the outer hcal slat tilt angle to : " << hcal2TiltAngle << " radian" << G4endl;
684  G4cout << " SetTitlViaNCross(" << ncross << ") setting the inner hcal slat tilt angle to : " << hcal1TiltAngle << " radian" << G4endl;
685  return;
686 }
687 
688 // Detector construction messengers for setting plate angles
689 //
690 void PHG4HcalPrototypeDetector::SetOuterHcalDPhi(G4double dphi)
691 {
692  hcal2DPhi = dphi;
693  G4cout << "In SetOuterHcalDPhi: " << hcal2DPhi << " is set!!! " << G4endl;
694 }
695 
696 void PHG4HcalPrototypeDetector::SetOuterPlateTiltAngle(G4double dtheta)
697 {
698  hcal2TiltAngle = dtheta;
699  G4cout << "In SetOuterPlateTiltAngle: " << hcal2TiltAngle << " is set!!! " << G4endl;
700 }
701 
702 void PHG4HcalPrototypeDetector::SetInnerHcalDPhi(G4double dphi)
703 {
704  hcal1DPhi = dphi;
705  G4cout << "In SetInnerHcalDPhi: " << hcal1DPhi << " is set!!! " << G4endl;
706 }
707 
708 void PHG4HcalPrototypeDetector::SetInnerPlateTiltAngle(G4double dtheta)
709 {
710  hcal1TiltAngle = dtheta;
711  G4cout << "In SetInnerPlateTiltAngle: " << hcal1TiltAngle << " is set!!! " << G4endl;
712 }