g4simulation/g4detectors/PHG4CEmcTestBeamDetector.cc

0001 #include "PHG4CEmcTestBeamDetector.h"
0002
0003 #include <g4main/PHG4Detector.h>  // for PHG4Detector
0004
0005 #include <phool/recoConsts.h>
0006
0007 #include <Geant4/G4Box.hh>
0008 #include <Geant4/G4Colour.hh>
0009 #include <Geant4/G4LogicalVolume.hh>
0010 #include <Geant4/G4PVPlacement.hh>
0011 #include <Geant4/G4RotationMatrix.hh>  // for G4RotationMatrix
0012 #include <Geant4/G4String.hh>          // for G4String
0013 #include <Geant4/G4SystemOfUnits.hh>
0014 #include <Geant4/G4ThreeVector.hh>  // for G4ThreeVector
0015 #include <Geant4/G4Transform3D.hh>  // for G4Transform3D
0016 #include <Geant4/G4Tubs.hh>
0017 #include <Geant4/G4Types.hh>  // for G4double
0018 #include <Geant4/G4VisAttributes.hh>
0019
0020 #include <algorithm>  // for copy
0021 #include <cmath>      // for cos, sin, NAN, acos, atan
0022 #include <iostream>   // for operator<<, ostringstream
0023 #include <sstream>
0024
0025 class G4Material;
0026 class G4VSolid;
0027 class PHCompositeNode;
0028
0029 using namespace std;
0030
0031 static double no_overlap = 0.0001 * cm;  // added safety margin against overlaps by using same boundary between volumes
0032
0033 PHG4CEmcTestBeamDetector::PHG4CEmcTestBeamDetector(PHG4Subsystem* subsys, PHCompositeNode* Node, const std::string& dnam, const int lyr)
0034   : PHG4Detector(subsys, Node, dnam)
0035   , gap(0.25 * mm)
0036   , place_in_x(0 * cm)
0037   , place_in_y(0 * cm)
0038   , place_in_z(0 * cm)
0039   , plate_x(135 * mm)
0040   , plate_z(135 * mm)
0041   , x_rot(0)
0042   , y_rot(0)
0043   , z_rot(0)
0044   , alpha(NAN)
0045   , inner_radius(NAN)
0046   , outer_radius(NAN)
0047   , tower_angular_coverage(NAN)
0048   , cemc_angular_coverage(NAN)
0049   , active_scinti_fraction(0.78)
0050   , sandwiches_per_tower(12)
0051   ,  // 12 tungsten/scintillator fiber snadwiches per tower
0052   num_towers(7)
0053   , active(0)
0054   , absorberactive(0)
0055   , layer(lyr)
0056   , blackhole(0)
0057 {
0058   w_dimension[0] = plate_x;
0059   w_dimension[1] = 0.5 * mm;
0060   w_dimension[2] = plate_z;
0061   sc_dimension[0] = plate_x;
0062   sc_dimension[1] = 1 * mm;
0063   sc_dimension[2] = plate_z;
0064   sandwich_thickness = 2 * w_dimension[1] + sc_dimension[1];  // two tungsten plats, one scintillator
0065   for (int i = 0; i < 4; i++)
0066   {
0067     sandwich_vol.push_back(nullptr);
0068   }
0069 }
0070
0071 //_______________________________________________________________
0072 //_______________________________________________________________
0073 int PHG4CEmcTestBeamDetector::IsInCEmcTestBeam(G4VPhysicalVolume* volume) const
0074 {
0075   if (active && volume == sandwich_vol[2])
0076   {
0077     return 1;
0078   }
0079   if (absorberactive && (volume == sandwich_vol[0] || volume == sandwich_vol[1]))
0080   {
0081     return -1;
0082   }
0083   if (absorberactive && sandwich_vol[3] && volume == sandwich_vol[3])
0084   {
0085     return -1;
0086   }
0087   return 0;
0088 }
0089
0090 void PHG4CEmcTestBeamDetector::ConstructMe(G4LogicalVolume* logicWorld)
0091 {
0092   CalculateGeometry();
0093   recoConsts* rc = recoConsts::instance();
0094   G4Material* Air = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
0095   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);
0096   G4LogicalVolume* cemc_log = new G4LogicalVolume(cemc_tub, Air, G4String("CEmc"), nullptr, nullptr, nullptr);
0097
0098   G4RotationMatrix cemc_rotm;
0099   // put our cemc at center displacement in x
0100   double radius_at_center = inner_radius + (outer_radius - inner_radius) / 2.;
0101   double xcenter = radius_at_center * cos(cemc_angular_coverage / 2.);
0102   double ycenter = radius_at_center * sin(cemc_angular_coverage / 2.);
0103   cemc_rotm.rotateX(x_rot);
0104   cemc_rotm.rotateY(y_rot);
0105   cemc_rotm.rotateZ(z_rot);
0106   new G4PVPlacement(G4Transform3D(cemc_rotm, G4ThreeVector(place_in_x - xcenter, place_in_y - ycenter, place_in_z)), cemc_log, "CEmc", logicWorld, false, false, OverlapCheck());
0107   G4VSolid* tower_tub = new G4Tubs("TowerTub", inner_radius - no_overlap, outer_radius + no_overlap, (w_dimension[2] + no_overlap) / 2., 0, tower_angular_coverage);
0108   G4LogicalVolume* tower_log = new G4LogicalVolume(tower_tub, Air, G4String("CEmcTower"), nullptr, nullptr, nullptr);
0109   //  G4VisAttributes* towerVisAtt = new G4VisAttributes();
0110   // towerVisAtt->SetVisibility(true);
0111   // towerVisAtt->SetForceSolid(true);
0112   // towerVisAtt->SetColour(G4Colour::Blue());
0113   //  tower_log->SetVisAttributes(towerVisAtt);
0114
0115   ostringstream tower_vol_name;
0116   for (int i = 0; i < 7; i++)
0117   {
0118     tower_vol_name << "CEmcTower_" << i;
0119     double phi = -i * tower_angular_coverage;
0120     G4RotationMatrix* tower_rotm = new G4RotationMatrix();
0121     tower_rotm->rotateZ(phi * rad);
0122     new G4PVPlacement(tower_rotm, G4ThreeVector(0, 0, 0), tower_log, tower_vol_name.str(), cemc_log, false, i, OverlapCheck());
0123     tower_vol_name.str("");
0124   }
0125   ConstructTowerVolume(tower_log);
0126   return;
0127 }
0128
0129 // here we build up the tower from the sandwichs (tungsten + scintillator)
0130 int PHG4CEmcTestBeamDetector::ConstructTowerVolume(G4LogicalVolume* tower_log)
0131 {
0132   recoConsts* rc = recoConsts::instance();
0133   G4Material* Air = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
0134   G4VSolid* sandwich_box = new G4Box("Sandwich_box",
0135                                      w_dimension[0] / 2., sandwich_thickness / 2., w_dimension[2] / 2.);
0136
0137   G4LogicalVolume* sandwich_log = new G4LogicalVolume(sandwich_box,
0138                                                       Air,
0139                                                       G4String("CEmcSandwich"),
0140                                                       nullptr, nullptr, nullptr);
0141   G4VisAttributes* sandwichVisAtt = new G4VisAttributes();
0142   sandwichVisAtt->SetVisibility(true);
0143   sandwichVisAtt->SetForceSolid(true);
0144   sandwichVisAtt->SetColour(G4Colour::White());
0145   sandwich_log->SetVisAttributes(sandwichVisAtt);
0146   ostringstream sandwich_name;
0147   for (int i = 0; i < 12; i++)
0148   {
0149     G4RotationMatrix* sandwich_rotm = new G4RotationMatrix();
0150     double phi = -i * alpha;
0151     sandwich_rotm->rotateZ(phi * rad);
0152     sandwich_name << "CEmcSandwich_" << i;
0153     double xshift = cos(phi) * (inner_radius + (outer_radius - inner_radius) / 2.);
0154     double yshift = -sin(phi) * (inner_radius + (outer_radius - inner_radius) / 2.);
0155     // we need to shift everything up by sandwich_thickness/2, calculate the shift in x
0156     // if we push the sandwich up by sandwich_thickness/2
0157     double xcorr = atan(phi) * sandwich_thickness / 2.;
0158     double ycorr = sandwich_thickness / 2.;
0159
0160     new G4PVPlacement(sandwich_rotm, G4ThreeVector(xshift + xcorr, yshift + ycorr, 0),
0161                       sandwich_log,
0162                       sandwich_name.str(),
0163                       tower_log, false, i, OverlapCheck());
0164     sandwich_name.str("");
0165   }
0166   ConstructSandwichVolume(sandwich_log);  // put W and scinti into sandwich
0167   return 0;
0168 }
0169
0170 // here we put a single tungsten + scintillator sandwich together
0171 int PHG4CEmcTestBeamDetector::ConstructSandwichVolume(G4LogicalVolume* sandwich)
0172 {
0173   vector<G4LogicalVolume*> block_logic;
0174   G4Material* AbsorberMaterial = GetDetectorMaterial("G4_W");
0175   G4Material* ScintiMaterial = GetDetectorMaterial("G4_POLYSTYRENE");
0176
0177   if (active_scinti_fraction > 1 || active_scinti_fraction < 0)
0178   {
0179     cout << "invalid active scintillator fraction " << active_scinti_fraction
0180          << " try between 0 and 1" << endl;
0181   }
0182
0183   double sc_active_thickness = sc_dimension[1] * active_scinti_fraction;
0184   double sc_passive_thickness = sc_dimension[1] - sc_active_thickness;
0185
0186   G4VSolid* block_w = new G4Box("Tungsten_box",
0187                                 w_dimension[0] / 2., w_dimension[1] / 2., w_dimension[2] / 2.);
0188   G4VSolid* block_sc = new G4Box("Scinti_box",
0189                                  sc_dimension[0] / 2., sc_active_thickness / 2., sc_dimension[2] / 2.);
0190   G4VSolid* block_passive_sc = nullptr;
0191   block_logic.push_back(new G4LogicalVolume(block_w,
0192                                             AbsorberMaterial,
0193                                             "Plate_log_W",
0194                                             nullptr, nullptr, nullptr));
0195   block_logic.push_back(new G4LogicalVolume(block_sc,
0196                                             ScintiMaterial,
0197                                             "Plate_log_Sc",
0198                                             nullptr, nullptr, nullptr));
0199   G4VisAttributes* matVis = new G4VisAttributes();
0200   G4VisAttributes* matVis1 = new G4VisAttributes();
0201   matVis->SetVisibility(true);
0202   matVis->SetForceSolid(true);
0203   matVis->SetColour(G4Colour::Red());
0204   matVis1->SetVisibility(true);
0205   matVis1->SetForceSolid(true);
0206   matVis1->SetColour(G4Colour::Green());
0207   block_logic[0]->SetVisAttributes(matVis);
0208   block_logic[1]->SetVisAttributes(matVis1);
0209
0210   // here is the math to set the positions of those tiles inside the sandwich. The center of the coordinate system is in the center
0211   // of the sandwich (w -> thickness of tungsten layer, sc_a -> thickness of actvie scinti, sc_p -> thickness of passive scinti)
0212   // distance to the bottom of the sandwich:
0213   // -(w + sc_a + sc_p + w)/2 = -w/2 -sc_a/2 -sc_p/2 -w/2
0214
0215   // implement the absorber at the bottom and at the top of the sandwich
0216   // moving up by w/2 for the lowest layer tungsten:
0217   // -w/2 -sc_a/2 -sc_p/2 -w/2 + w/2 = -w/2 -sc_a/2 -sc_p/2 = -(w+sc)/2
0218   // where sc_a+sc_p = sc = scintillator thickness (sc_dimension[1])
0219
0220   sandwich_vol[0] = new G4PVPlacement(nullptr, G4ThreeVector(0, -(w_dimension[1] + sc_dimension[1]) / 2., 0),
0221                                       block_logic[0],
0222                                       "CEmc_W_plate_down",
0223                                       sandwich, false, 0, OverlapCheck());
0224
0225   // top of the sandwich - starting at the bottom and add the tungsten and scintillator layer and half tungsten
0226   // -w/2 -sc/2 -w/2 + w + sc + w/2 = +sc/2 +w/2 = (w+sc)/2
0227   sandwich_vol[1] = new G4PVPlacement(nullptr, G4ThreeVector(0, (w_dimension[1] + sc_dimension[1]) / 2., 0),
0228                                       block_logic[0],
0229                                       "CEmc_W_plate_up",
0230                                       sandwich, false, 0, OverlapCheck());
0231
0232   // since we split the scintillator into an active and passive part to accomodate
0233   // for the scintillator fibers not occupying 100% of the volume.
0234   // this is mocked up by 2 separate layers of plastic,  active scintillator is the middle layer
0235   // again going to the bottom of the sandwich box (sc_a -> active sc, sc_p -> passive sc)
0236   //  -w/2 -sc_a/2 -sc_p/2 -w/2
0237   // -(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:
0238   // -w/2 -sc_a/2 - sc_p/2 -w/2 + w + sc_a/2 = -sc_p/2
0239   // if fraction of active sc is 1, sc_passive_thickness is zero
0240   sandwich_vol[2] = new G4PVPlacement(nullptr, G4ThreeVector(0, -sc_passive_thickness / 2., 0),
0241                                       block_logic[1],
0242                                       "CEmc_Scinti_plate",
0243                                       sandwich, false, 0, OverlapCheck());
0244
0245   if (sc_passive_thickness > 0)
0246   {
0247     G4VisAttributes* matVis2 = new G4VisAttributes();
0248     matVis1->SetVisibility(true);
0249     matVis1->SetForceSolid(true);
0250     matVis1->SetColour(G4Colour::Blue());
0251     block_passive_sc = new G4Box("Passive_Scinti_box",
0252                                  sc_dimension[0] / 2., sc_passive_thickness / 2., sc_dimension[2] / 2.);
0253     block_logic.push_back(new G4LogicalVolume(block_passive_sc,
0254                                               ScintiMaterial,
0255                                               "Plate_log_Passive_Sc",
0256                                               nullptr, nullptr, nullptr));
0257     block_logic[2]->SetVisAttributes(matVis2);
0258     // here we go again - bottome of sandwich box is
0259     //  -(w + sc_a + sc_p +w)/2, now add w and sc_a and half of sc_p:
0260     //  -w/2 -sc_a/2 - sc_p/2 -w/2 + w + sc_a + sc_p/2 = sc_a/2
0261     sandwich_vol[3] = new G4PVPlacement(nullptr, G4ThreeVector(0, sc_active_thickness / 2., 0),
0262                                         block_logic[2],
0263                                         "CEmc_Passive_Si_plate",
0264                                         sandwich, false, 0, OverlapCheck());
0265   }
0266   else
0267   {
0268     sandwich_vol[3] = nullptr;
0269   }
0270   return 0;
0271 }
0272
0273 void PHG4CEmcTestBeamDetector::CalculateGeometry()
0274 {
0275   // calculate the inner/outer radius of the g4tub using the test setup numbers
0276   // 1mm scintillator, 0.5mm tungsten, 0.25 mm gap at the back of the detector
0277   // the geometry with figure is explained in our wiki:
0278   // https://www.phenix.bnl.gov/WWW/offline/wikioff/index.php/CEmc
0279   // get the alpha angle via law of cos (a is the 0.25 mm gap):
0280   // a^2 = b^2+c^2 - 2bc*cos(alpha)
0281   // cos(alpha) =  (b^2+c^2 - a^2)/2bc
0282   // with b=c
0283   // cos(alpha) = 1-(a^2/2b^2)
0284   alpha = acos(1 - (gap * gap) / (2 * plate_x * plate_x));
0285   // inner radius
0286   inner_radius = sandwich_thickness / tan(alpha);
0287   outer_radius = sqrt((inner_radius + plate_x) * (inner_radius + plate_x) + sandwich_thickness * sandwich_thickness);
0288
0289   // twice no_overlap to apply safety margin also for towers
0290   // angular coverage of 1 tower is 12 sandwiches = sandwiches_per_tower*alpha
0291   tower_angular_coverage = sandwiches_per_tower * alpha + 0.00005 * deg;  // safety margin added to remove 142 nm overlap
0292   // cemc prototype has 7 towers
0293   cemc_angular_coverage = num_towers * tower_angular_coverage;
0294   return;
0295 }