g4simulation/g4detectors/PHG4HcalDetector.cc

0001 #include "PHG4HcalDetector.h"
0002 #include "PHG4CylinderGeomContainer.h"
0003 #include "PHG4CylinderGeomv3.h"
0004
0005 #include <g4main/PHG4Detector.h>  // for PHG4Detector
0006 #include <g4main/PHG4Utils.h>
0007
0008 #include <phool/PHCompositeNode.h>
0009 #include <phool/PHIODataNode.h>
0010 #include <phool/PHNode.h>          // for PHNode
0011 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0012 #include <phool/PHObject.h>        // for PHObject
0013 #include <phool/getClass.h>
0014
0015 #include <Geant4/G4Colour.hh>  // for G4Colour
0016 #include <Geant4/G4Cons.hh>
0017 #include <Geant4/G4ExtrudedSolid.hh>
0018 #include <Geant4/G4LogicalVolume.hh>
0019 #include <Geant4/G4Material.hh>
0020 #include <Geant4/G4PVPlacement.hh>
0021 #include <Geant4/G4PhysicalConstants.hh>
0022 #include <Geant4/G4RotationMatrix.hh>  // for G4RotationMatrix
0023 #include <Geant4/G4String.hh>          // for G4String
0024 #include <Geant4/G4SubtractionSolid.hh>
0025 #include <Geant4/G4SystemOfUnits.hh>  // for cm, deg
0026 #include <Geant4/G4ThreeVector.hh>    // for G4ThreeVector
0027 #include <Geant4/G4Transform3D.hh>    // for G4Transform3D
0028 #include <Geant4/G4Tubs.hh>
0029 #include <Geant4/G4TwoVector.hh>
0030 #include <Geant4/G4VisAttributes.hh>
0031
0032 #include <algorithm>  // for max
0033 #include <cmath>      // for sin, cos, sqrt, M_PI, asin
0034 #include <cstdlib>    // for exit
0035 #include <iostream>   // for operator<<, basic_ostream
0036 #include <sstream>
0037 #include <utility>  // for pair
0038 #include <vector>   // for vector
0039
0040 class PHG4CylinderGeom;
0041
0042 // uncomment if you want to make a graphics display where the slats are visible
0043 // it makes them stick out of the hcal for visibility
0044 // NEVER EVER RUN REAL SIMS WITH THIS
0045 //#define DISPLAY
0046
0047 int PHG4HcalDetector::INACTIVE = -100;
0048 //_______________________________________________________________
0049 // note this inactive thickness is ~1.5% of a radiation length
0050 PHG4HcalDetector::PHG4HcalDetector(PHG4Subsystem* subsys, PHCompositeNode* Node, const std::string& dnam, const int lyr)
0051   : PHG4Detector(subsys, Node, dnam)
0052   , TrackerMaterial(nullptr)
0053   , TrackerThickness(100 * cm)
0054   , cylinder_logic(nullptr)
0055   , cylinder_physi(nullptr)
0056   , radius(100 * cm)
0057   , length(100 * cm)
0058   , xpos(0 * cm)
0059   , ypos(0 * cm)
0060   , zpos(0 * cm)
0061   , _sciTilt(0)
0062   , _sciWidth(0.6 * cm)
0063   , _sciNum(100)
0064   , _sciPhi0(0)
0065   , _region(nullptr)
0066   , active(0)
0067   , absorberactive(0)
0068   , layer(lyr)
0069 {
0070 }
0071
0072 //_______________________________________________________________
0073 //_______________________________________________________________
0074 int PHG4HcalDetector::IsInCylinderActive(const G4VPhysicalVolume* volume)
0075 {
0076   //  std::cout << "checking detector" << std::endl;
0077   if (active && box_vol.find(volume) != box_vol.end())
0078   {
0079     return box_vol.find(volume)->second;
0080   }
0081   if (absorberactive && volume == cylinder_physi)
0082   {
0083     return -1;
0084   }
0085   return INACTIVE;
0086 }
0087
0088 //_______________________________________________________________
0089 void PHG4HcalDetector::ConstructMe(G4LogicalVolume* logicWorld)
0090 {
0091   TrackerMaterial = GetDetectorMaterial(material);
0092
0093   G4Tubs* _cylinder_solid = new G4Tubs(G4String(GetName()),
0094                                        radius,
0095                                        radius + TrackerThickness,
0096                                        length / 2.0, 0, twopi);
0097   double innerlength = PHG4Utils::GetLengthForRapidityCoverage(radius) * 2;
0098   double deltalen = (length - innerlength) / 2.;  // length difference on one side
0099   double cone_size_multiplier = 1.01;             // 1 % larger
0100   double cone_thickness = TrackerThickness * cone_size_multiplier;
0101   double inner_cone_radius = radius - ((cone_thickness - TrackerThickness) / 2.);
0102   double cone_length = deltalen * cone_size_multiplier;
0103   G4Cons* cone2 = new G4Cons("conehead2",
0104                              inner_cone_radius, inner_cone_radius,
0105                              inner_cone_radius, inner_cone_radius + cone_thickness,
0106                              cone_length / 2.0, 0, twopi);
0107   G4Cons* cone1 = new G4Cons("conehead",
0108                              inner_cone_radius, inner_cone_radius + cone_thickness,
0109                              inner_cone_radius, inner_cone_radius,
0110                              cone_length / 2.0, 0, twopi);
0111   double delta_len = cone_length - deltalen;
0112   G4ThreeVector zTransneg(0, 0, -(length - cone_length + delta_len) / 2.0);
0113   G4ThreeVector zTranspos(0, 0, (length - cone_length + delta_len) / 2.0);
0114   G4SubtractionSolid* subtraction_tmp =
0115       new G4SubtractionSolid("Cylinder-Cone", _cylinder_solid, cone1, nullptr, zTransneg);
0116   G4SubtractionSolid* subtraction =
0117       new G4SubtractionSolid("Cylinder-Cone-Cone", subtraction_tmp, cone2, nullptr, zTranspos);
0118   cylinder_logic = new G4LogicalVolume(subtraction,
0119                                        TrackerMaterial,
0120                                        G4String(GetName()),
0121                                        nullptr, nullptr, nullptr);
0122   G4VisAttributes* VisAtt = new G4VisAttributes();
0123   VisAtt->SetColour(G4Colour::Grey());
0124   VisAtt->SetVisibility(true);
0125   VisAtt->SetForceSolid(true);
0126   cylinder_logic->SetVisAttributes(VisAtt);
0127
0128   cylinder_physi = new G4PVPlacement(nullptr, G4ThreeVector(xpos, ypos, zpos),
0129                                      cylinder_logic,
0130                                      G4String(GetName()),
0131                                      logicWorld, false, false, OverlapCheck());
0132   // Figure out corners of scintillator inside the containing G4Tubs.
0133   // Work our way around the scintillator cross section in a counter
0134   // clockwise fashion: ABCD
0135   double r1 = radius;
0136   double r2 = radius + TrackerThickness;
0137
0138   // The coordinates of the inner corners of the scintillator
0139   double x4 = r1;
0140   double y4 = _sciWidth / 2.0;
0141
0142   double x1 = r1;
0143   double y1 = -y4;
0144
0145   double a = _sciTilt * M_PI / 180.0;
0146
0147   // The parametric equation for the line from A along the side of the
0148   // scintillator is (x,y) = (x1,y1) + u * (cos(a), sin(a))
0149   double A = 1.0;
0150   double B = 2 * (x1 * cos(a) + y1 * sin(a));
0151   double C = x1 * x1 + y1 * y1 - r2 * r2;
0152   double D = B * B - 4 * A * C;
0153
0154   // The only sensible solution, given our definitions, is u > 0.
0155   double u = (-B + sqrt(D)) / 2 * A;
0156
0157   // Now we can determine one of the outer corners
0158   double x2 = x1 + u * cos(a);
0159   double y2 = y1 + u * sin(a);
0160
0161   // Similar procedure for (x3,y3) as for (x2,y2)
0162   A = 1.0;
0163   B = 2 * (x4 * cos(a) + y4 * sin(a));
0164   C = x4 * x4 + y4 * y4 - r2 * r2;
0165   D = B * B - 4 * A * C;
0166   u = (-B + sqrt(D)) / 2 * A;
0167
0168   double x3 = x4 + u * cos(a);
0169   double y3 = y4 + u * sin(a);
0170
0171   // Now we've got a four-sided "z-section".
0172   G4TwoVector v1(x1, y1);
0173   G4TwoVector v2(x2, y2);
0174   G4TwoVector v3(x3, y3);
0175   G4TwoVector v4(x4, y4);
0176
0177   std::vector<G4TwoVector> vertexes;
0178   vertexes.push_back(v1);
0179   vertexes.push_back(v2);
0180   vertexes.push_back(v3);
0181   vertexes.push_back(v4);
0182
0183   G4TwoVector zero(0, 0);
0184   // if you want to make displays where the structure of the hcal is visible
0185   // add 20 cm to the length of the scintillators
0186 #ifdef DISPLAY
0187   double blength = length + 20;
0188 #else
0189   double blength = length;
0190 #endif
0191   G4ExtrudedSolid* _box_solid = new G4ExtrudedSolid("_BOX",
0192                                                     vertexes,
0193                                                     blength / 2.0,
0194                                                     zero, 1.0,
0195                                                     zero, 1.0);
0196
0197   //  double boxlen_half = GetLength(_sciTilt * M_PI / 180.);
0198   //  G4Box* _box_solid = new G4Box("_BOX", boxlen_half, _sciWidth / 2.0, length / 2.0);
0199
0200   G4Material* boxmat = GetDetectorMaterial("G4_POLYSTYRENE");
0201   G4SubtractionSolid* subtractionbox_tmp =
0202       new G4SubtractionSolid("Box-Cone", _box_solid, cone1, nullptr, zTransneg);
0203   G4SubtractionSolid* subtractionbox =
0204       new G4SubtractionSolid("Box-Cone-Cone", subtractionbox_tmp, cone2, nullptr, zTranspos);
0205   G4LogicalVolume* box_logic = new G4LogicalVolume(subtractionbox,
0206                                                    boxmat, G4String("BOX"),
0207                                                    nullptr, nullptr, nullptr);
0208   VisAtt = new G4VisAttributes();
0209   PHG4Utils::SetColour(VisAtt, "G4_POLYSTYRENE");
0210   VisAtt->SetVisibility(true);
0211   VisAtt->SetForceSolid(true);
0212   box_logic->SetVisAttributes(VisAtt);
0213
0214   double phi_increment = 360. / _sciNum;
0215   std::ostringstream slatname;
0216   for (int i = 0; i < _sciNum; i++)
0217   {
0218     double phi = (i + _sciPhi0) * phi_increment;
0219     G4ThreeVector myTrans = G4ThreeVector(0, 0, 0);
0220     G4RotationMatrix Rot(0, 0, 0);
0221     Rot.rotateZ(phi * deg);
0222     slatname.str("");
0223     slatname << "SLAT_" << i;
0224     G4VPhysicalVolume* box_vol_tmp = new G4PVPlacement(G4Transform3D(Rot, G4ThreeVector(myTrans)),
0225                                                        box_logic,
0226                                                        G4String(slatname.str()),
0227                                                        cylinder_logic, false, false, OverlapCheck());
0228     box_vol[box_vol_tmp] = i;
0229   }
0230   if (active)
0231   {
0232     std::ostringstream geonode;
0233     if (superdetector != "NONE")
0234     {
0235       geonode << "CYLINDERGEOM_" << superdetector;
0236     }
0237     else
0238     {
0239       geonode << "CYLINDERGEOM_" << detector_type << "_" << layer;
0240     }
0241     PHG4CylinderGeomContainer* geo = findNode::getClass<PHG4CylinderGeomContainer>(topNode(), geonode.str());
0242     if (!geo)
0243     {
0244       geo = new PHG4CylinderGeomContainer();
0245       PHNodeIterator iter(topNode());
0246       PHCompositeNode* runNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "RUN"));
0247       PHIODataNode<PHObject>* newNode = new PHIODataNode<PHObject>(geo, geonode.str(), "PHObject");
0248       runNode->addNode(newNode);
0249     }
0250     // here in the detector class we have internal units, convert to cm
0251     // before putting into the geom object
0252     PHG4CylinderGeom* mygeom = new PHG4CylinderGeomv3(radius / cm, (zpos - length / 2.) / cm, (zpos + length / 2.) / cm, TrackerThickness / cm, _sciNum, _sciTilt * M_PI / 180.0, _sciPhi0 * M_PI / 180.0);
0253     geo->AddLayerGeom(layer, mygeom);
0254     //    geo->identify();
0255   }
0256 }
0257
0258 double
0259 PHG4HcalDetector::GetLength(const double phi) const
0260 {
0261   double c = radius + TrackerThickness / 2.;
0262   double b = radius;
0263   double singamma = sin(phi) * c / b;
0264   double gamma = M_PI - asin(singamma);
0265   double alpha = M_PI - gamma - phi;
0266   double a = c * sin(alpha) / singamma;
0267   return a;
0268 }
0269
0270 void PHG4HcalDetector::Print(const std::string& /*what*/) const
0271 {
0272   std::cout << "radius: " << radius << std::endl;
0273   return;
0274 }