g4simulation/g4caloprototype/PHG4HcalPrototypeDetector.cc

0001 // This is G4 detector implementation for the hcal prototype
0002 // Created on 1/27/2014, Liang, HeXC
0003 // Updated on 3/21/2014, Liang, HeXC
0004 // Updated on 4/18/2014, Liang, HeXC, Updating detector construction (proper tilting angles!)
0005
0006 #include "PHG4HcalPrototypeDetector.h"
0007
0008 #include "PHG4HcalPrototypeDetectorMessenger.h"  // for PHG4HcalPrototypeDet...
0009
0010 #include <g4main/PHG4Detector.h>  // for PHG4Detector
0011
0012 #include <Geant4/G4Box.hh>
0013 #include <Geant4/G4Colour.hh>
0014 #include <Geant4/G4LogicalVolume.hh>
0015 #include <Geant4/G4Material.hh>
0016 #include <Geant4/G4NistManager.hh>
0017 #include <Geant4/G4PVPlacement.hh>
0018 #include <Geant4/G4RotationMatrix.hh>  // for G4RotationMatrix
0019 #include <Geant4/G4RunManager.hh>
0020 #include <Geant4/G4String.hh>         // for G4String
0021 #include <Geant4/G4SystemOfUnits.hh>  // for cm, mm, deg, rad
0022 #include <Geant4/G4ThreeVector.hh>    // for G4ThreeVector
0023 #include <Geant4/G4Transform3D.hh>    // for G4Transform3D
0024 #include <Geant4/G4Trap.hh>
0025 #include <Geant4/G4Types.hh>  // for G4double, G4int, G4bool
0026 #include <Geant4/G4UnionSolid.hh>
0027 #include <Geant4/G4VisAttributes.hh>
0028 #include <Geant4/G4ios.hh>  // for G4cout, G4endl
0029
0030 #include <cmath>  // for cos, sin, M_PI, asin
0031 #include <sstream>
0032
0033 class PHCompositeNode;
0034
0035 using namespace std;
0036
0037 // static double no_overlap = 0.00015 * cm; // added safety margin against overlaps by using same boundary between volumes
0038
0039 PHG4HcalPrototypeDetector::PHG4HcalPrototypeDetector(PHG4Subsystem* subsys, PHCompositeNode* Node, const std::string& dnam, const int lyr)
0040   : PHG4Detector(subsys, Node, dnam)
0041   , nScint360(15)
0042   ,  // This is legacy variable in case one wants to build a cylindrical calorimeter
0043   nHcal1Layers(16)
0044   , nHcal2Layers(nHcal1Layers)
0045   , hcal2ScintSizeX(2.54 * 26.1 * cm)
0046   ,  // from Don's drawing
0047   hcal2ScintSizeY(0.85 * cm)
0048   ,  // Changed from 8.0cm to 8.5cm following
0049   hcal2ScintSizeZ(2.54 * 29.53 * cm)
0050   ,  // from Don's drawing
0051   hcal1ScintSizeX(2.54 * 12.43 * cm)
0052   ,  // from Don's drawing
0053   hcal1ScintSizeY(0.85 * cm)
0054   ,  // Changed from 8.0cm to 8.5cm following
0055   hcal1ScintSizeZ(2.54 * 17.1 * cm)
0056   ,  // from Don's drawing
0057   hcal1TiltAngle(0.27)
0058   , hcal2TiltAngle(0.14)
0059   , hcal1DPhi(0.27)
0060   ,  // calculated based on Don's drawing.  (twopi/16.0);
0061   hcal2DPhi(0.288)
0062   ,  // calculated based on Don's drawing.  (twopi/16.0);
0063   hcal1RadiusIn(1855 * mm)
0064   , hcal2RadiusIn(hcal1RadiusIn + hcal1ScintSizeX)
0065   ,
0066   // The frame box for the HCAL.
0067   hcalBoxSizeX(1.1 * (hcal2ScintSizeX + hcal1ScintSizeX))
0068   , hcalBoxSizeY(2.54 * 40.0 * cm)
0069   ,  // need to get more accurate dimension for the box
0070   hcalBoxSizeZ(1.1 * hcal2ScintSizeZ)
0071   , hcalBoxRotationAngle_z(0.0 * rad)
0072   ,  // Rotation along z-axis
0073   hcalBoxRotationAngle_y(0.0 * rad)
0074   ,  // Rotation along y-axis
0075   hcal2Abs_dxa(2.54 * 27.1 * cm)
0076   ,  // these numbers are from Don't drawings
0077   hcal2Abs_dxb(hcal2Abs_dxa)
0078   , hcal2Abs_dya(2.54 * 1.099 * cm)
0079   , hcal2Abs_dyb(2.54 * 1.776 * cm)
0080   , hcal2Abs_dz(hcal2ScintSizeZ)
0081   , hcal1Abs_dxa(hcal1ScintSizeX)
0082   ,  // these numbers are from Don't drawings
0083   hcal1Abs_dxb(hcal1ScintSizeX)
0084   , hcal1Abs_dya(2.54 * 0.787 * cm)
0085   , hcal1Abs_dyb(2.54 * 1.115 * cm)
0086   , hcal1Abs_dz(hcal1ScintSizeZ)
0087   , hcalJunctionSizeX(2.54 * 1.0 * cm)
0088   , hcalJunctionSizeY(hcal2ScintSizeY)
0089   , hcalJunctionSizeZ(hcal2Abs_dz)
0090   , physiWorld(nullptr)
0091   , logicWorld(nullptr)
0092   , logicHcalBox(nullptr)
0093   , solidHcalBox(nullptr)
0094   , physiHcalBox(nullptr)
0095   , logicHcal2ScintLayer(nullptr)
0096   , logicHcal1ScintLayer(nullptr)
0097   , logicHcal2AbsLayer(nullptr)
0098   , logicHcal1AbsLayer(nullptr)
0099   , world_mat(nullptr)
0100   , steel(nullptr)
0101   , scint_mat(nullptr)
0102   , active(0)
0103   , absorberactive(0)
0104   , layer(lyr)
0105   , blackhole(0)
0106 {
0107   // create commands for interactive definition of the detector
0108   fDetectorMessenger = new PHG4HcalPrototypeDetectorMessenger(this);
0109 }
0110
0111 //_______________________________________________________________
0112 //_______________________________________________________________
0113 int PHG4HcalPrototypeDetector::IsInHcalPrototype(G4VPhysicalVolume* /*volume*/) const
0114 {
0115   return 0;  // not sure what value to return for now.
0116 }
0117
0118 void PHG4HcalPrototypeDetector::ConstructMe(G4LogicalVolume* world)
0119 {
0120   logicWorld = world;
0121
0122   DefineMaterials();
0123
0124   ConstructDetector();
0125 }
0126
0127 void PHG4HcalPrototypeDetector::DefineMaterials()
0128 {
0129   // Water is defined from NIST material database
0130   G4NistManager* nist = G4NistManager::Instance();
0131
0132   world_mat = nist->FindOrBuildMaterial("G4_AIR");
0133
0134   steel = nist->FindOrBuildMaterial("G4_Fe");  // Need to fix this *******
0135   scint_mat = nist->FindOrBuildMaterial("G4_POLYSTYRENE");
0136
0137   G4cout << *(G4Material::GetMaterialTable()) << G4endl;
0138 }
0139
0140 // We now build our detector solids, etc
0141 //
0142 G4VPhysicalVolume* PHG4HcalPrototypeDetector::ConstructDetector()
0143 {
0144   // Option to switch on/off checking of volumes overlaps
0145   //
0146   G4bool checkOverlaps = true;
0147
0148   // HCAL Frame Box for enclosing the inner and the outer hcal
0149   // which allows us to rotate the whole calorimeter easily
0150   solidHcalBox = new G4Box("HcalBox",                                              //its name
0151                            hcalBoxSizeX / 2, hcalBoxSizeY / 2, hcalBoxSizeZ / 2);  //its size
0152
0153   logicHcalBox = new G4LogicalVolume(solidHcalBox,  //its solid
0154                                      world_mat,     //its material
0155                                      "HcalBox");    //its name
0156
0157   // Place the HcalBox inside the World
0158   G4ThreeVector boxVector = G4ThreeVector(hcal1RadiusIn + hcalBoxSizeX / 2.0, 0, 0);
0159   G4RotationMatrix rotBox = G4RotationMatrix();
0160
0161   rotBox.rotateZ(hcalBoxRotationAngle_z);  // Rotate the whole calorimeter box along z-axis
0162   rotBox.rotateY(hcalBoxRotationAngle_y);  // Rotate the whole calorimeter box along z-axis
0163
0164   G4Transform3D boxTransform = G4Transform3D(rotBox, boxVector);
0165
0166   physiHcalBox = new G4PVPlacement(boxTransform,    // rotation + positioning
0167                                    logicHcalBox,    //its logical volume
0168                                    "HcalBox",       //its name
0169                                    logicWorld,      //its mother  volume
0170                                    false,           //no boolean operation
0171                                    0,               //copy number
0172                                    checkOverlaps);  //checking overlaps
0173
0174   // HCal Box visualization attributes
0175   G4VisAttributes* hcalBoxVisAtt = new G4VisAttributes(G4Colour(0.0, 0.0, 1.0));  // blue
0176   hcalBoxVisAtt->SetVisibility(true);
0177   //hcalBoxVisAtt->SetForceSolid(true);
0178   logicHcalBox->SetVisAttributes(hcalBoxVisAtt);
0179
0180   // Make outer hcal section
0181   //
0182   // Build scintillator layer box as a place holder for the scintillator sheets
0183   // Make it 0.05cm thinner than the true gap width for avoiding geometry overlaps
0184   G4Box* hcal2ScintLayer = new G4Box("hcal2ScintLayer",                                                             //its name
0185                                      hcal2ScintSizeX / 2, (hcal2ScintSizeY - 0.05 * cm) / 2, hcal2ScintSizeZ / 2);  //its size
0186
0187   logicHcal2ScintLayer =
0188       new G4LogicalVolume(hcal2ScintLayer,     // its solid name
0189                           world_mat,           // material, the same as the material of the world valume
0190                           "hcal2ScintLayer");  // its name
0191
0192   // Scintillator visualization attributes
0193   G4VisAttributes* scintVisAtt = new G4VisAttributes(G4Colour(1.0, 0.0, 0.0));  //Red
0194   scintVisAtt->SetVisibility(true);
0195   //scintVisAtt->SetForceSolid(true);
0196   logicHcal2ScintLayer->SetVisAttributes(scintVisAtt);
0197
0198   // Construct outer scintillator 1U
0199   G4double outer1UpDz = 649.8 * mm;
0200   G4double outer1UpDy1 = 2 * 0.35 * cm;
0201   G4double outer1UpDx2 = 179.3 * mm;
0202   G4double outer1UpDx4 = 113.2 * mm;
0203
0204   G4Trap* outer1USheetSolid = new G4Trap("outer1USheet",  //its name
0205                                          outer1UpDy1, outer1UpDz,
0206                                          outer1UpDx2, outer1UpDx4);  //its size
0207
0208   G4LogicalVolume* logicOuter1USheet = new G4LogicalVolume(outer1USheetSolid,
0209                                                            scint_mat,
0210                                                            "outer1USheet");
0211
0212   // Scintillator visualization attributes
0213   G4VisAttributes* scintSheetVisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 0.0));  //White
0214   scintSheetVisAtt->SetVisibility(true);
0215   scintSheetVisAtt->SetForceSolid(true);
0216
0217   logicOuter1USheet->SetVisAttributes(scintSheetVisAtt);
0218
0219   // Detector placement
0220   // Position the outer right 1U sheet
0221   G4ThreeVector threeVecOuter1U_1 = G4ThreeVector(0 * cm, 0 * cm, -0.252 * (outer1UpDx2 + outer1UpDx4));
0222   G4RotationMatrix rot1U_1 = G4RotationMatrix();
0223   rot1U_1.rotateZ(90 * deg);
0224   rot1U_1.rotateX(-90 * deg);
0225
0226   G4Transform3D transformOuter1U_1 = G4Transform3D(rot1U_1, threeVecOuter1U_1);
0227
0228   new G4PVPlacement(transformOuter1U_1,
0229                     logicOuter1USheet,
0230                     "outer1USheet",
0231                     logicHcal2ScintLayer,
0232                     false,
0233                     0,  // Copy one
0234                     checkOverlaps);
0235
0236   // Detector placement
0237   // Position the outer left 1U sheet
0238   G4ThreeVector threeVecOuter1U_2 = G4ThreeVector(0 * cm, 0 * cm, 0.252 * (outer1UpDx2 + outer1UpDx4));
0239   G4RotationMatrix rot1U_2 = G4RotationMatrix();
0240   rot1U_2.rotateZ(90 * deg);
0241   rot1U_2.rotateX(90 * deg);
0242
0243   G4Transform3D transformOuter1U_2 = G4Transform3D(rot1U_2, threeVecOuter1U_2);
0244
0245   new G4PVPlacement(transformOuter1U_2,
0246                     logicOuter1USheet,
0247                     "outer1USheet",
0248                     logicHcal2ScintLayer,
0249                     false,
0250                     1,  // Copy two
0251                     checkOverlaps);
0252
0253   // Construct outer scintillator 2U
0254   G4double outer2UpDz = 0.5 * 649.8 * mm;
0255   G4double outer2UpTheta = 8.8 * M_PI / 180.;
0256   G4double outer2UpPhi = 0.0 * M_PI / 180.;
0257   G4double outer2UpDy1 = 0.35 * cm;
0258   G4double outer2UpDy2 = 0.35 * cm;
0259   G4double outer2UpDx1 = 0.5 * 179.3 * mm, outer2UpDx2 = 0.5 * 179.3 * mm;
0260   G4double outer2UpDx3 = 0.5 * 113.2 * mm, outer2UpDx4 = 0.5 * 113.2 * mm;
0261   G4double outer2UpAlp1 = 0. * M_PI / 180., outer2UpAlp2 = 0. * M_PI / 180;
0262
0263   G4Trap* outer2USheetSolid = new G4Trap("outer2USheet",
0264                                          outer2UpDz,
0265                                          outer2UpTheta,
0266                                          outer2UpPhi,
0267                                          outer2UpDy1,
0268                                          outer2UpDx1,
0269                                          outer2UpDx2,
0270                                          outer2UpAlp1,
0271                                          outer2UpDy2,
0272                                          outer2UpDx3,
0273                                          outer2UpDx4,
0274                                          outer2UpAlp2);
0275
0276   G4LogicalVolume* logicOuter2USheet = new G4LogicalVolume(outer2USheetSolid,
0277                                                            scint_mat,
0278                                                            "outer2USheet");
0279
0280   logicOuter2USheet->SetVisAttributes(scintSheetVisAtt);
0281
0282   // Detector placement
0283   // Position the right most sheet
0284   G4ThreeVector threeVecOuter2U_1 = G4ThreeVector(0 * cm, 0 * cm, -0.755 * (outer1UpDx2 + outer1UpDx4));
0285   G4RotationMatrix rot2U_1 = G4RotationMatrix();
0286   rot2U_1.rotateY(-90 * deg);
0287
0288   G4Transform3D transformOuter2U_1 = G4Transform3D(rot2U_1, threeVecOuter2U_1);
0289
0290   new G4PVPlacement(transformOuter2U_1,
0291                     logicOuter2USheet,
0292                     "outer2USheet",
0293                     logicHcal2ScintLayer,
0294                     false,
0295                     0,  // copy one
0296                     checkOverlaps);
0297
0298   // Detector placement
0299   // Position the outer left most sheet
0300   G4ThreeVector threeVecOuter2U_2 = G4ThreeVector(0 * cm, 0 * cm, 0.755 * (outer1UpDx2 + outer1UpDx4));
0301   G4RotationMatrix rot2U_2 = G4RotationMatrix();
0302   rot2U_2.rotateY(-90 * deg);
0303   rot2U_2.rotateX(180 * deg);
0304
0305   G4Transform3D transformOuter2U_2 = G4Transform3D(rot2U_2, threeVecOuter2U_2);
0306
0307   new G4PVPlacement(transformOuter2U_2,
0308                     logicOuter2USheet,
0309                     "outer2USheet",
0310                     logicHcal2ScintLayer,
0311                     false,
0312                     1,  // copy two
0313                     checkOverlaps);
0314
0315   G4Trap* hcal2AbsLayer =
0316       new G4Trap("hcal2AbsLayer",  //its name
0317                  hcal2Abs_dz, hcal2Abs_dxa,
0318                  hcal2Abs_dyb, hcal2Abs_dya);  //its size
0319
0320   // Add extra absorber materials at the junction
0321   G4Box* hcalJunction = new G4Box("HcalJunction",                                                        //its name
0322                                   hcalJunctionSizeX / 2, hcalJunctionSizeY / 2, hcalJunctionSizeZ / 2);  //its size
0323
0324   // define the transformation for placing the junction to the inner side of the hcal2 absorber layer
0325   G4ThreeVector threeVecJunction = G4ThreeVector(-hcal2Abs_dya / 2.0 - 1.1 * hcalJunctionSizeY, hcal2Abs_dxa / 2.0 - hcalJunctionSizeX / 2, 0.0 * cm);
0326   G4RotationMatrix rotJunction = G4RotationMatrix();
0327   rotJunction.rotateZ(90 * deg);
0328   rotJunction.rotateX(0 * deg);
0329   G4Transform3D transformJunction = G4Transform3D(rotJunction, threeVecJunction);
0330
0331   // attach the junction to the absorber using G4UnionSolid
0332   G4UnionSolid* hcal2AbsLayerJunct = new G4UnionSolid("hcal2AbsLayerJunct", hcal2AbsLayer, hcalJunction, transformJunction);
0333
0334   logicHcal2AbsLayer = new G4LogicalVolume(hcal2AbsLayerJunct,  //  hcal2AbsLayer,
0335                                            steel,
0336                                            "hcal2AbsLayer");
0337
0338   // hcal2Aborber visualization attributes
0339   G4VisAttributes* hcal2AbsVisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 1.0));
0340   hcal2AbsVisAtt->SetVisibility(true);
0341   logicHcal2AbsLayer->SetVisAttributes(hcal2AbsVisAtt);
0342
0343   // place scintillator layers insides the HCal2 absorber volume
0344   //
0345
0346   G4double theta = 0.;
0347   G4double theta2 = hcal2TiltAngle;                                                                               // I am not convinced that I got the tilt angle right.  So I simply forced it to
0348   G4double rScintLayerCenter = hcal2RadiusIn + hcal2ScintSizeX / 2.0 + 2.54 * 1.0 * cm;                           // move the scintillator toward the rear end of HCAL2
0349   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
0350   G4double rAbsLayerCenter = hcal2RadiusIn + hcal2Abs_dxa / 2.0;
0351   G4double RmidAbsLayerX = 0.81 * (hcal2Abs_dxa - hcal1ScintSizeX) / 2.0 + 2 * 2.54 * cm;
0352   G4int LayerNum = 1;
0353   G4double xPadding = 0.0;
0354   G4double yPadding = 0.0;
0355   G4double tiltPadding = 0.0;
0356
0357   for (G4int iLayer = -nHcal2Layers / 2; iLayer < nHcal2Layers / 2; iLayer++)
0358   {
0359     //
0360     // Place outer absorber layer first
0361     //
0362     theta = hcal2DPhi / nHcal2Layers * iLayer;
0363     G4double xposShift = rAbsLayerCenter * (cos(theta) - 1.0);
0364     G4double yposShift = rAbsLayerCenter * sin(theta);
0365     xPadding = 0.013 * RmidAbsLayerX * LayerNum;  // adding an extra padding for placing the absorber and scintillator
0366     G4ThreeVector absTrans = G4ThreeVector(RmidAbsLayerX + xposShift - xPadding, yposShift, 0);
0367     G4RotationMatrix rotAbsLayer = G4RotationMatrix();
0368     rotAbsLayer.rotateY(90 * deg);
0369     rotAbsLayer.rotateX(90 * deg);
0370     rotAbsLayer.rotateY(-90 * deg);
0371     tiltPadding = LayerNum * 0.005 * rad;
0372     rotAbsLayer.rotateZ((iLayer + nHcal2Layers / 2) * 0.02 * rad + tiltPadding);  // Added a funge increment factor 0.01
0373
0374     G4Transform3D transformAbs = G4Transform3D(rotAbsLayer, absTrans);
0375
0376     new G4PVPlacement(transformAbs,        //rotation,position
0377                       logicHcal2AbsLayer,  //its logical volume
0378                       "hcal2AbsLayer",     //its name
0379                       logicHcalBox,        //its mother  volume
0380                       false,               //no boolean operation
0381                       LayerNum - 1,        //copy number
0382                       checkOverlaps);      //overlaps checking
0383     //
0384     // Place outer hcal scintillator layer
0385     //
0386     theta += 0.5 * hcal2DPhi / nHcal2Layers;
0387     G4cout << "M_PI: " << M_PI << "     theta: " << theta << "    TileAngle: " << theta2 << G4endl;
0388     xposShift = rScintLayerCenter * (cos(theta) - 1.0);
0389     yposShift = rScintLayerCenter * sin(theta);
0390     G4ThreeVector myTrans = G4ThreeVector(RmidScintLayerX + xposShift - xPadding, yposShift + 0.3 * cm, 0);  // hcal2RadiusIn + hcal1ScintSizeX/2.0
0391     G4RotationMatrix rotm = G4RotationMatrix();
0392     rotm.rotateZ((iLayer + nHcal2Layers / 2 + 1) * 0.020 * rad + tiltPadding);  // Added a funge increment factor 0.02
0393     G4Transform3D transform = G4Transform3D(rotm, myTrans);
0394
0395     G4cout << "  iLayer " << iLayer << G4endl;
0396
0397     new G4PVPlacement(transform,             //rotation,position
0398                       logicHcal2ScintLayer,  //its logical volume
0399                       "hcal2ScintLayer",     //its name
0400                       logicHcalBox,          //its mother volume
0401                       false,                 //no boolean operation
0402                       LayerNum - 1,          //copy number
0403                       checkOverlaps);        //checking overlaps
0404     LayerNum++;
0405   }
0406
0407   // Make INNER hcal section
0408   //
0409   // Build scintillator layer box as a place holder for the scintillator sheets
0410   // Make it 0.05cm thinner than the true gap width for avoiding geometry overlaps
0411   G4Box* hcal1ScintLayer = new G4Box("hcal1ScintLayer",                                                             //its name
0412                                      hcal1ScintSizeX / 2, (hcal1ScintSizeY - 0.05 * cm) / 2, hcal1ScintSizeZ / 2);  //its size
0413
0414   logicHcal1ScintLayer =
0415       new G4LogicalVolume(hcal1ScintLayer,     // its solid name
0416                           world_mat,           // simply use world material
0417                           "hcal1ScintLayer");  // its name
0418
0419   logicHcal1ScintLayer->SetVisAttributes(scintVisAtt);
0420
0421   // Constructing scintillator sheets
0422   // Construct inner scintillator 1U
0423   // The numbers are read off from Don's drawings
0424   G4double inner1UpDz = 316.8 * mm;
0425   G4double inner1UpDy1 = 2 * 0.35 * cm;
0426   G4double inner1UpDx2 = 108.6 * mm;
0427   G4double inner1UpDx4 = 77.4 * mm;
0428
0429   G4Trap* inner1USheetSolid = new G4Trap("inner1USheet",  //its name
0430                                          inner1UpDy1, inner1UpDz,
0431                                          inner1UpDx2, inner1UpDx4);  //its size
0432
0433   G4LogicalVolume* logicInner1USheet = new G4LogicalVolume(inner1USheetSolid,
0434                                                            scint_mat,
0435                                                            "inner1USheet");
0436   logicInner1USheet->SetVisAttributes(scintSheetVisAtt);
0437
0438   // Detector placement
0439   // Position the inner right 1U sheet
0440   G4ThreeVector threeVecInner1U_1_inner = G4ThreeVector(0 * cm, 0 * cm, -0.252 * (inner1UpDx2 + inner1UpDx4));
0441
0442   G4Transform3D transformInner1U_1 = G4Transform3D(rot1U_1, threeVecInner1U_1_inner);
0443
0444   new G4PVPlacement(transformInner1U_1,
0445                     logicInner1USheet,
0446                     "inner1USheet",
0447                     logicHcal1ScintLayer,
0448                     false,
0449                     0,  // Copy one
0450                     checkOverlaps);
0451
0452   // Detector placement
0453   // Position the inner left 1U sheet
0454   G4ThreeVector threeVecInner1U_2_inner = G4ThreeVector(0 * cm, 0 * cm, 0.252 * (inner1UpDx2 + inner1UpDx4));
0455   //G4RotationMatrix rot1U_2  = G4RotationMatrix();
0456   //rot1U_2.rotateZ(90*deg);
0457   //rot1U_2.rotateX(90*deg);
0458
0459   G4Transform3D transformInner1U_2 = G4Transform3D(rot1U_2, threeVecInner1U_2_inner);
0460
0461   new G4PVPlacement(transformInner1U_2,
0462                     logicInner1USheet,
0463                     "inner1USheet",
0464                     //logicWorld,
0465                     logicHcal1ScintLayer,
0466                     false,
0467                     1,  // Copy two
0468                     checkOverlaps);
0469
0470   // Construct inner scintillator 2U
0471   G4double inner2UpDz = 0.5 * 316.8 * mm;
0472   G4double inner2UpTheta = 8.8 * M_PI / 180.;
0473   G4double inner2UpPhi = 0.0 * M_PI / 180.;
0474   G4double inner2UpDy1 = 0.35 * cm;
0475   G4double inner2UpDy2 = 0.35 * cm;
0476   G4double inner2UpDx1 = 0.5 * 108.6 * mm, inner2UpDx2 = 0.5 * 108.6 * mm;
0477   G4double inner2UpDx3 = 0.5 * 77.4 * mm, inner2UpDx4 = 0.5 * 77.4 * mm;
0478   G4double inner2UpAlp1 = 0. * M_PI / 180., inner2UpAlp2 = 0. * M_PI / 180;
0479
0480   G4Trap* inner2USheetSolid = new G4Trap("inner2USheet",
0481                                          inner2UpDz,
0482                                          inner2UpTheta,
0483                                          inner2UpPhi,
0484                                          inner2UpDy1,
0485                                          inner2UpDx1,
0486                                          inner2UpDx2,
0487                                          inner2UpAlp1,
0488                                          inner2UpDy2,
0489                                          inner2UpDx3,
0490                                          inner2UpDx4,
0491                                          inner2UpAlp2);
0492
0493   G4LogicalVolume* logicInner2USheet = new G4LogicalVolume(inner2USheetSolid,
0494                                                            scint_mat,
0495                                                            "inner2USheet");
0496
0497   logicInner2USheet->SetVisAttributes(scintSheetVisAtt);
0498
0499   // Detector placement
0500   // Position the right most sheet
0501   G4ThreeVector threeVecInner2U_1_inner = G4ThreeVector(0 * cm, 0 * cm, -0.755 * (inner1UpDx2 + inner1UpDx4));
0502   //G4RotationMatrix rot2U_1  = G4RotationMatrix();
0503   //rot2U_1.rotateY(-90*deg);
0504
0505   G4Transform3D transformInner2U_1 = G4Transform3D(rot2U_1, threeVecInner2U_1_inner);
0506
0507   new G4PVPlacement(transformInner2U_1,
0508                     logicInner2USheet,
0509                     "inner2USheet",
0510                     //logicWorld,
0511                     logicHcal1ScintLayer,
0512                     false,
0513                     0,  // copy one
0514                     checkOverlaps);
0515
0516   // Detector placement
0517   // Position the inner left most sheet
0518   G4ThreeVector threeVecInner2U_2_inner = G4ThreeVector(0 * cm, 0 * cm, 0.755 * (inner1UpDx2 + inner1UpDx4));
0519   //G4RotationMatrix rot2U_2  = G4RotationMatrix();
0520   //rot2U_2.rotateY(-90*deg);
0521   //rot2U_2.rotateX(180*deg);
0522
0523   G4Transform3D transformInner2U_2 = G4Transform3D(rot2U_2, threeVecInner2U_2_inner);
0524
0525   new G4PVPlacement(transformInner2U_2,
0526                     logicInner2USheet,
0527                     "inner2USheet",
0528                     //logicWorld,
0529                     logicHcal1ScintLayer,
0530                     false,
0531                     1,  // copy two
0532                     checkOverlaps);
0533
0534   // Build hcal1 absorber layer in trapezoid shape
0535   //
0536   /*
0537   G4double hcal1Abs_dxa = hcal1ScintSizeX; // hcal1Abs_dxb = hcal1ScintSizeX;
0538   G4double hcal1Abs_dya = 0.787*2.54*cm, hcal1Abs_dyb = 1.115*2.54*cm;
0539   // these numbers are from Don't drawings
0540   G4double hcal1Abs_dz  = hcal1ScintSizeZ;
0541   */
0542   G4Trap* hcal1AbsLayer =
0543       new G4Trap("hcal1AbsLayer",  //its name
0544                  hcal1Abs_dz, hcal1Abs_dxa,
0545                  hcal1Abs_dyb, hcal1Abs_dya);  //its size
0546
0547   logicHcal1AbsLayer = new G4LogicalVolume(hcal1AbsLayer,
0548                                            steel,
0549                                            "hcal1AbsLayer");
0550
0551   logicHcal1AbsLayer->SetVisAttributes(hcal2AbsVisAtt);
0552
0553   // place scintillator layers insides the HCal1 absorber volume
0554   //
0555   // Calculate the title angle
0556
0557   theta2 = hcal1TiltAngle;
0558   rAbsLayerCenter = hcal1RadiusIn + hcal1ScintSizeX / 2.0;
0559   //  RmidAbsLayerX =  (hcal2Abs_dxa - hcal1ScintSizeX)/2.0;
0560   RmidAbsLayerX = hcal2Abs_dxa / 2.0;
0561   LayerNum = 1;  // These three parameters are for making fine adjustment of the placement of inner hcal
0562   yPadding = -2.54 * 0.8 * cm;
0563   for (G4int iLayer = -nHcal2Layers / 2 - 1; iLayer < nHcal2Layers / 2 - 1; iLayer++)
0564   {  // shift one layer down
0565     //
0566     // Place inner absorber layer first
0567     //
0568     theta = hcal1DPhi / nHcal1Layers * iLayer;  // add an off set 0.01 rad
0569     G4double xposShift = rAbsLayerCenter * (cos(theta) - 1.0);
0570     G4double yposShift = rAbsLayerCenter * sin(theta);
0571     xPadding = 0.005 * RmidAbsLayerX * LayerNum - 2.54 * 1.9 * cm;
0572     G4ThreeVector absTrans = G4ThreeVector(-RmidAbsLayerX * cos(theta) + xposShift - xPadding, yposShift + yPadding, 0);
0573     G4RotationMatrix rotAbsLayer = G4RotationMatrix();
0574     rotAbsLayer.rotateY(90 * deg);
0575     rotAbsLayer.rotateX(90 * deg);
0576     rotAbsLayer.rotateY(-90 * deg);
0577     tiltPadding = LayerNum * 0.005 * rad;
0578     rotAbsLayer.rotateZ((theta - theta2) * rad + tiltPadding);
0579     //rotAbsLayer.rotateZ(-(iLayer + nHcal2Layers/2)*0.02*rad + tiltPadding);
0580
0581     G4Transform3D transformAbs = G4Transform3D(rotAbsLayer, absTrans);
0582
0583     new G4PVPlacement(transformAbs,        //rotation,position
0584                       logicHcal1AbsLayer,  //its logical volume
0585                       "hcal1AbsLayer",     //its name
0586                       logicHcalBox,        //its mother  volume
0587                       false,               //no boolean operation
0588                       LayerNum - 1,        //copy number
0589                       checkOverlaps);      //overlaps checking
0590
0591     //
0592     // Place inner hcal scintillator layer
0593     //
0594
0595     theta += 0.5 * hcal1DPhi / nHcal1Layers;
0596     //G4cout << "M_PI: " << M_PI << "     theta: " << theta << "    TileAngle: " << theta2 << G4endl;
0597     //G4double phi = iLayer*dPhi;
0598     //    G4ThreeVector myTrans = G4ThreeVector(-RmidAbsLayerX*cos(theta) + xposShift,  yposShift, 0);
0599     xposShift = rAbsLayerCenter * (cos(theta) - 1.0);
0600     yposShift = rAbsLayerCenter * sin(theta);
0601     G4ThreeVector myTrans = G4ThreeVector(-RmidAbsLayerX * cos(theta) + xposShift - xPadding, yposShift + yPadding, 0);
0602     //G4cout << " x: " << Rmid*cos(theta) << "    y: " << Rmid*sin(theta) << "   1.0*rad: " << 1.8*rad << G4endl;
0603     G4RotationMatrix rotm = G4RotationMatrix();
0604     rotm.rotateZ((theta - theta2 + 0.01) * rad + tiltPadding);  // Added a funge increment factor 0.01
0605     //rotm.rotateZ(-(iLayer + nHcal2Layers/2+1)*0.02*rad + tiltPadding);        // Added a funge increment factor 0.01
0606     G4Transform3D transform = G4Transform3D(rotm, myTrans);
0607
0608     G4cout << "  iLayer " << iLayer << G4endl;
0609
0610     new G4PVPlacement(transform,             //rotation,position
0611                       logicHcal1ScintLayer,  //its logical volume
0612                       "hcal1ScintLayer",     //its name
0613                       logicHcalBox,          //its mother volume: logicHcal1Ab
0614                       false,                 //no boolean operation
0615                       LayerNum - 1,          //copy number
0616                       checkOverlaps);        //checking overlaps
0617     LayerNum++;
0618   }
0619
0620   return physiWorld;
0621 }
0622
0623 void PHG4HcalPrototypeDetector::CalculateGeometry()
0624 {
0625   return;
0626 }
0627
0628 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0629
0630 void PHG4HcalPrototypeDetector::SetMaterial(G4String /*materialChoice*/)
0631 {
0632 }
0633
0634 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0635
0636 void PHG4HcalPrototypeDetector::UpdateGeometry()
0637 {
0638   G4RunManager::GetRunManager()->PhysicsHasBeenModified();
0639   G4RunManager::GetRunManager()->DefineWorldVolume(ConstructDetector());
0640 }
0641
0642 // The following code is copied from the sPHENIX G4 simulation
0643 void PHG4HcalPrototypeDetector::SetTiltViaNcross(const int ncross)
0644 {
0645   G4double sign = 1;
0646   if (ncross < 0)
0647   {
0648     sign = -1;
0649   }
0650   G4int ncr = fabs(ncross);
0651
0652   // Determine title angle for the outer hcal
0653   //
0654   G4double cSide = hcal2RadiusIn + hcal2ScintSizeX / 2;
0655   G4double bSide = hcal2RadiusIn + hcal2ScintSizeX;
0656
0657   G4double alpha = 0;
0658   if (ncr > 1)
0659   {
0660     alpha = (360. / nHcal2Layers * M_PI / 180.) * (ncr - 1) / 2.0;
0661   }
0662   else
0663   {
0664     alpha = (360. / nHcal2Layers * M_PI / 180.) / 2.;
0665   }
0666
0667   G4double sinbSide = sin(alpha) * bSide / (sqrt(bSide * bSide + cSide * cSide - 2 * bSide * cSide * cos(alpha)));
0668   G4double beta = asin(sinbSide);  // This is the slat angle
0669
0670   hcal2TiltAngle = beta * sign;
0671
0672   // Determine title angle for the inner hcal
0673   //
0674   cSide = hcal1RadiusIn + hcal1ScintSizeX / 2;
0675   bSide = hcal1RadiusIn + hcal1ScintSizeX;
0676
0677   sinbSide = sin(alpha) * bSide / (sqrt(bSide * bSide + cSide * cSide - 2.0 * bSide * cSide * cos(alpha)));
0678   beta = asin(sinbSide);  // This is the slat angle
0679
0680   hcal1TiltAngle = beta * sign;
0681
0682   G4cout << " alpha : " << alpha << G4endl;
0683   G4cout << " SetTitlViaNCross(" << ncross << ") setting the outer hcal slat tilt angle to : " << hcal2TiltAngle << " radian" << G4endl;
0684   G4cout << " SetTitlViaNCross(" << ncross << ") setting the inner hcal slat tilt angle to : " << hcal1TiltAngle << " radian" << G4endl;
0685   return;
0686 }
0687
0688 // Detector construction messengers for setting plate angles
0689 //
0690 void PHG4HcalPrototypeDetector::SetOuterHcalDPhi(G4double dphi)
0691 {
0692   hcal2DPhi = dphi;
0693   G4cout << "In SetOuterHcalDPhi: " << hcal2DPhi << " is set!!! " << G4endl;
0694 }
0695
0696 void PHG4HcalPrototypeDetector::SetOuterPlateTiltAngle(G4double dtheta)
0697 {
0698   hcal2TiltAngle = dtheta;
0699   G4cout << "In SetOuterPlateTiltAngle: " << hcal2TiltAngle << " is set!!! " << G4endl;
0700 }
0701
0702 void PHG4HcalPrototypeDetector::SetInnerHcalDPhi(G4double dphi)
0703 {
0704   hcal1DPhi = dphi;
0705   G4cout << "In SetInnerHcalDPhi: " << hcal1DPhi << " is set!!! " << G4endl;
0706 }
0707
0708 void PHG4HcalPrototypeDetector::SetInnerPlateTiltAngle(G4double dtheta)
0709 {
0710   hcal1TiltAngle = dtheta;
0711   G4cout << "In SetInnerPlateTiltAngle: " << hcal1TiltAngle << " is set!!! " << G4endl;
0712 }