g4simulation/g4detectors/PHG4ZDCDetector.cc

0001 #include "PHG4ZDCDetector.h"
0002
0003 #include "PHG4ZDCDefs.h"
0004 #include "PHG4ZDCDisplayAction.h"
0005
0006 #include <phparameter/PHParameters.h>
0007
0008 #include <g4gdml/PHG4GDMLUtility.hh>
0009
0010 #include <g4main/PHG4Detector.h>       // for PHG4Detector
0011 #include <g4main/PHG4DisplayAction.h>  // for PHG4DisplayAction
0012 #include <g4main/PHG4Subsystem.h>
0013
0014 #include <phool/recoConsts.h>
0015
0016 #include <TSystem.h>
0017
0018 #include <Geant4/G4Box.hh>
0019 #include <Geant4/G4LogicalVolume.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>
0026 #include <Geant4/G4ThreeVector.hh>  // for G4ThreeVector
0027 #include <Geant4/G4Transform3D.hh>  // for G4Transform3D
0028 #include <Geant4/G4Tubs.hh>
0029 #include <Geant4/G4Types.hh>            // for G4double, G4int
0030 #include <Geant4/G4VPhysicalVolume.hh>  // for G4VPhysicalVolume
0031
0032 #include <cassert>
0033 #include <cmath>
0034 #include <iostream>
0035 #include <sstream>
0036
0037 class G4Material;
0038 class G4VSolid;
0039 class PHCompositeNode;
0040
0041 //_______________________________________________________________________
0042 PHG4ZDCDetector::PHG4ZDCDetector(PHG4Subsystem* subsys, PHCompositeNode* Node, PHParameters* parameters, const std::string& dnam, const int detid)
0043   : PHG4Detector(subsys, Node, dnam)
0044   , m_DisplayAction(dynamic_cast<PHG4ZDCDisplayAction*>(subsys->GetDisplayAction()))
0045   , m_Params(parameters)
0046   , m_GdmlConfig(PHG4GDMLUtility::GetOrMakeConfigNode(Node))
0047   , m_Angle(M_PI_4 * radian)
0048   , m_TPlate(2.3 * mm)
0049   , m_HPlate(400.0 * mm)
0050   , m_WPlate(100.0 * mm)
0051   , m_TAbsorber(5.0 * mm)
0052   , m_HAbsorber(150.0 * mm)
0053   , m_WAbsorber(100.0 * mm)
0054   , m_DFiber(0.5 * mm)
0055   , m_HFiber(400.0 * mm)
0056   , m_WFiber(100.0 * mm)
0057   , m_GFiber(0.0001 * mm)
0058   , m_Gap(0.2 * mm)
0059   , m_TSMD(10.0 * mm)
0060   , m_HSMD(160.0 * mm)
0061   , m_WSMD(105.0 * mm)
0062   , m_RHole(63.9 * mm)
0063   , m_TWin(4.7 * mm)
0064   , m_RWin(228.60 * mm)
0065   , m_PlaceHole(122.56 * mm)
0066   , m_Pxwin(0.0 * mm)
0067   , m_Pywin(0.0 * mm)
0068   , m_Pzwin(915.5 * mm)
0069   , m_NMod(3)
0070   , m_NLay(27)
0071   , m_ActiveFlag(m_Params->get_int_param("active"))
0072   , m_AbsorberActiveFlag(m_Params->get_int_param("absorberactive"))
0073   , m_SupportActiveFlag(m_Params->get_int_param("supportactive"))
0074   , m_Layer(detid)
0075   , m_SuperDetector("NONE")
0076 {
0077   assert(m_GdmlConfig);
0078 }
0079
0080 //_______________________________________________________________________
0081 int PHG4ZDCDetector::IsInZDC(G4VPhysicalVolume* volume) const
0082 {
0083   G4LogicalVolume* mylogvol = volume->GetLogicalVolume();
0084
0085   if (m_ActiveFlag)
0086   {
0087     if (m_ScintiLogicalVolSet.find(mylogvol) != m_ScintiLogicalVolSet.end())
0088     {
0089       return 1;
0090     }
0091     if (m_FiberLogicalVolSet.find(mylogvol) != m_FiberLogicalVolSet.end())
0092     {
0093       return 2;
0094     }
0095   }
0096   if (m_AbsorberActiveFlag)
0097   {
0098     if (m_AbsorberLogicalVolSet.find(mylogvol) != m_AbsorberLogicalVolSet.end())
0099     {
0100       return -1;
0101     }
0102   }
0103   if (m_SupportActiveFlag)
0104   {
0105     if (m_SupportLogicalVolSet.find(mylogvol) != m_SupportLogicalVolSet.end())
0106     {
0107       return -2;
0108     }
0109   }
0110   return 0;
0111 }
0112
0113 //_______________________________________________________________________
0114 void PHG4ZDCDetector::ConstructMe(G4LogicalVolume* logicWorld)
0115 {
0116   if (Verbosity() > 0)
0117   {
0118     std::cout << "PHG4ZDCDetector: Begin Construction" << std::endl;
0119   }
0120
0121   if (m_Layer != PHG4ZDCDefs::NORTH &&
0122       m_Layer != PHG4ZDCDefs::SOUTH)
0123   {
0124     std::cout << "use either PHG4ZDCDefs::NORTH or PHG4ZDCDefs::SOUTH for ZDC Subsystem" << std::endl;
0125     gSystem->Exit(1);
0126     return;
0127   }
0128
0129   recoConsts* rc = recoConsts::instance();
0130   G4Material* WorldMaterial = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
0131   G4Material* Fe = GetDetectorMaterial("G4_Fe");
0132   G4Material* W = GetDetectorMaterial("G4_W");
0133   G4Material* PMMA = GetDetectorMaterial("G4_PLEXIGLASS");
0134   G4Material* Scint = GetDetectorMaterial("G4_POLYSTYRENE");  // place holder
0135
0136   G4double TGap = m_DFiber + m_Gap;
0137   G4double Mod_Length = 2 * m_TPlate + m_NLay * (TGap + m_TAbsorber);
0138   G4double Det_Length = m_NMod * Mod_Length + m_TSMD;
0139   G4double RTT = 1. / sin(m_Angle);
0140   G4double First_Pos = -RTT * Det_Length / 2;
0141   G4double Room = 3.5 * mm;
0142   G4double Mother_2Z = RTT * Det_Length + 2. * (m_HFiber - m_HAbsorber / 2.) * cos(m_Angle);
0143   /* Create exit windows */
0144   G4VSolid* ExitWindow_nocut_solid = new G4Tubs(G4String("ExitWindow_nocut_solid"),
0145                                                 0.0, m_RWin, m_TWin, 0.0, CLHEP::twopi);
0146
0147   G4VSolid* Hole_solid = new G4Tubs(G4String("Hole_solid"),
0148                                     0.0, m_RHole, 2 * m_TWin, 0.0, CLHEP::twopi);
0149   G4VSolid* ExitWindow_1cut_solid = new G4SubtractionSolid("ExitWindow_1cut_solid", ExitWindow_nocut_solid, Hole_solid, nullptr, G4ThreeVector(m_PlaceHole, 0, 0));
0150
0151   G4VSolid* ExitWindow_2cut_solid = new G4SubtractionSolid("ExitWindow_2cut_solid", ExitWindow_1cut_solid, Hole_solid, nullptr, G4ThreeVector(-m_PlaceHole, 0, 0));
0152
0153   G4LogicalVolume* ExitWindow_log = new G4LogicalVolume(ExitWindow_2cut_solid, GetDetectorMaterial("G4_STAINLESS-STEEL"), G4String("ExitWindow_log"), nullptr, nullptr, nullptr);
0154
0155   GetDisplayAction()->AddVolume(ExitWindow_log, "Window");
0156   m_SupportLogicalVolSet.insert(ExitWindow_log);
0157   G4RotationMatrix Window_rotm;
0158   Window_rotm.rotateX(m_Params->get_double_param("rot_x") * deg);
0159   Window_rotm.rotateY(m_Params->get_double_param("rot_y") * deg);
0160   Window_rotm.rotateZ(m_Params->get_double_param("rot_z") * deg);
0161
0162   if (m_Layer == PHG4ZDCDefs::NORTH)
0163   {
0164     new G4PVPlacement(G4Transform3D(Window_rotm, G4ThreeVector(m_Pxwin, m_Pywin, m_Params->get_double_param("place_z") * cm - m_Pzwin)),
0165                       ExitWindow_log, "Window_North", logicWorld, false, PHG4ZDCDefs::NORTH, OverlapCheck());
0166   }
0167
0168   else if (m_Layer == PHG4ZDCDefs::SOUTH)
0169   {
0170     new G4PVPlacement(G4Transform3D(Window_rotm, G4ThreeVector(m_Pxwin, m_Pywin, -m_Params->get_double_param("place_z") * cm + m_Pzwin)),
0171                       ExitWindow_log, "Window_South", logicWorld, false, PHG4ZDCDefs::SOUTH, OverlapCheck());
0172   }
0173   /* ZDC detector here */
0174   /* Create the box envelope = 'world volume' for ZDC */
0175
0176   G4double Mother_X = m_WAbsorber / 2. + Room;
0177   G4double Mother_Y = (m_HFiber - m_HAbsorber / 2.) * sin(m_Angle) + Room;
0178   G4double Mother_Z = Mother_2Z / 2. + Room;
0179   G4VSolid* ZDC_envelope_solid = new G4Box(G4String("ZDC_envelope_solid"),
0180                                            Mother_X,
0181                                            Mother_Y,
0182                                            Mother_Z);
0183
0184   G4LogicalVolume* ZDC_envelope_log = new G4LogicalVolume(ZDC_envelope_solid, WorldMaterial, G4String("ZDC_envelope_log"), nullptr, nullptr, nullptr);
0185
0186   /* Define visualization attributes */
0187   GetDisplayAction()->AddVolume(ZDC_envelope_log, "Envelope");
0188
0189   /* Define rotation attributes for envelope cone */
0190   G4RotationMatrix ZDC_rotm;
0191   ZDC_rotm.rotateX(m_Params->get_double_param("rot_x") * deg);
0192   ZDC_rotm.rotateY(m_Params->get_double_param("rot_y") * deg);
0193   ZDC_rotm.rotateZ(m_Params->get_double_param("rot_z") * deg);
0194
0195   /* Create logical volumes for a plate to contain fibers */
0196   G4VSolid* fiber_plate_solid = new G4Box(G4String("fiber_plate_solid"),
0197                                           m_WFiber / 2.,
0198                                           m_HFiber / 2.,
0199                                           TGap / 2.);
0200
0201   G4LogicalVolume* fiber_plate_log = new G4LogicalVolume(fiber_plate_solid, WorldMaterial, G4String("fiber_plate_log"), nullptr, nullptr, nullptr);
0202   GetDisplayAction()->AddVolume(fiber_plate_log, "fiber_plate_air");
0203   /*  front and back plate */
0204   G4VSolid* fb_plate_solid = new G4Box(G4String("fb_plate_solid"),
0205                                        m_WPlate / 2.,
0206                                        m_HPlate / 2.,
0207                                        m_TPlate / 2.);
0208
0209   G4LogicalVolume* fb_plate_log = new G4LogicalVolume(fb_plate_solid, Fe, G4String("fb_plate_log"), nullptr, nullptr, nullptr);
0210   m_SupportLogicalVolSet.insert(fb_plate_log);
0211   GetDisplayAction()->AddVolume(fb_plate_log, "FrontBackPlate");
0212
0213   /*  absorber */
0214   G4VSolid* absorber_solid = new G4Box(G4String("absorber_solid"),
0215                                        m_WAbsorber / 2.,
0216                                        m_HAbsorber / 2.,
0217                                        m_TAbsorber / 2.);
0218
0219   G4LogicalVolume* absorber_log = new G4LogicalVolume(absorber_solid, W, G4String("absorber_log"), nullptr, nullptr, nullptr);
0220   m_AbsorberLogicalVolSet.insert(absorber_log);
0221   GetDisplayAction()->AddVolume(absorber_log, "Absorber");
0222
0223   /* SMD */
0224   // volume that contains scintillators
0225   G4VSolid* SMD_solid = new G4Box(G4String("SMD_solid"),
0226                                   m_WSMD / 2.,
0227                                   m_HSMD / 2.,
0228                                   m_TSMD / 2.);
0229
0230   G4LogicalVolume* SMD_log = new G4LogicalVolume(SMD_solid, WorldMaterial, G4String("SMD_log"), nullptr, nullptr, nullptr);
0231   GetDisplayAction()->AddVolume(SMD_log, "SMD");
0232   // small scintillators block
0233   G4double scintx = 15 * mm;
0234   G4double scinty = 20 * mm;
0235   G4VSolid* Scint_solid = new G4Box(G4String("Scint_solid"),
0236                                     scintx / 2.,
0237                                     scinty / 2.,
0238                                     m_TSMD / 2.);
0239
0240   G4LogicalVolume* Scint_log = new G4LogicalVolume(Scint_solid, Scint, G4String("Scint_log"), nullptr, nullptr, nullptr);
0241   m_ScintiLogicalVolSet.insert(Scint_log);
0242   GetDisplayAction()->AddVolume(Scint_log, "Scint_solid");
0243
0244   // put scintillators in the SMD volume
0245   double scint_XPos = -m_WSMD / 2.;
0246   double scint_Xstep = scintx / 2.;
0247   double scint_Ystep = scinty / 2.;
0248   int Nx = m_WSMD / scintx;
0249   int Ny = m_HSMD / scinty;
0250
0251   for (int i = 0; i < Nx; i++)
0252   {
0253     scint_XPos += scint_Xstep;
0254     double scint_YPos = -m_HSMD / 2.;
0255     for (int j = 0; j < Ny; j++)
0256     {
0257       int copyno = Nx * j + i;
0258       scint_YPos += scint_Ystep;
0259       G4RotationMatrix SMDRotation;
0260       new G4PVPlacement(G4Transform3D(SMDRotation, G4ThreeVector(scint_XPos, scint_YPos, 0.0)),
0261                         Scint_log,
0262                         "single_scint",
0263                         SMD_log,
0264                         false, copyno, OverlapCheck());
0265
0266       scint_YPos += scint_Ystep;
0267     }
0268     scint_XPos += scint_Xstep;
0269   }
0270
0271   /* Create logical volumes for fibers */
0272   G4VSolid* single_fiber_solid = new G4Tubs(G4String("single_fiber_solid"),
0273                                             0.0, (m_DFiber / 2.) - m_GFiber, (m_HFiber / 2.), 0.0, CLHEP::twopi);
0274
0275   G4LogicalVolume* single_fiber_log = new G4LogicalVolume(single_fiber_solid, PMMA, G4String("single_fiber_log"), nullptr, nullptr, nullptr);
0276   m_FiberLogicalVolSet.insert(single_fiber_log);
0277   GetDisplayAction()->AddVolume(single_fiber_log, "Fiber");
0278
0279   /* Rotation Matrix for fibers */
0280   G4RotationMatrix* FiberRotation = new G4RotationMatrix();
0281   FiberRotation->rotateX(90. * deg);
0282
0283   /* Place fibers in the fiber plate */
0284   G4double fiber_XPos = -m_WFiber / 2.;
0285   G4double fiber_step = m_DFiber / 2.;
0286   int Nfiber = m_WFiber / m_DFiber;
0287
0288   for (int i = 0; i < Nfiber; i++)
0289   {
0290     fiber_XPos += fiber_step;
0291     int copyno = i;
0292
0293     new G4PVPlacement(FiberRotation, G4ThreeVector(fiber_XPos, 0.0, 0.0),
0294                       single_fiber_log,
0295                       G4String("single_fiber_scint"),
0296                       fiber_plate_log,
0297                       false, copyno, OverlapCheck());
0298     fiber_XPos += fiber_step;
0299   }
0300
0301   /* Rotation for plates in ZDC */
0302   G4RotationMatrix* PlateRotation = new G4RotationMatrix();
0303
0304   PlateRotation->rotateX(-m_Angle);
0305
0306   /* construct ZDC */
0307   G4double Plate_Step = m_TPlate * RTT;
0308   G4double Absorber_Step = m_TAbsorber * RTT;
0309   G4double Gap_Step = TGap * RTT;
0310   G4double SMD_Step = m_TSMD * RTT;
0311   G4double ZPos = First_Pos - Plate_Step / 2.;
0312   G4double Gap_YPos = (m_HFiber - m_HAbsorber) / 2. * sin(m_Angle);
0313   G4double Gap_ZPos = (m_HFiber - m_HAbsorber) / 2. * cos(m_Angle);
0314
0315   G4double Plate_YPos = (m_HPlate - m_HAbsorber) / 2. * sin(m_Angle);
0316   G4double Plate_ZPos = (m_HPlate - m_HAbsorber) / 2. * cos(m_Angle);
0317
0318   G4double SMD_YPos = (m_HSMD - m_HAbsorber) / 2. * sin(m_Angle);
0319   G4double SMD_ZPos = (m_HSMD - m_HAbsorber) / 2. * cos(m_Angle);
0320
0321   /* start the loop: for every module ---  front plate-absorber-fiber plate-absorber-.....-fiber plate-back plate */
0322   int copyno_plate = 0;
0323   for (int i = 0; i < m_NMod; i++)
0324   {
0325     // place the SMD in between the 1st and 2nd module
0326     if (i == 1)
0327     {
0328       ZPos += (SMD_Step / 2.);
0329       new G4PVPlacement(PlateRotation, G4ThreeVector(0.0, SMD_YPos, SMD_ZPos + ZPos),
0330                         SMD_log,
0331                         G4String("SMD"),
0332                         ZDC_envelope_log,
0333                         false, 0, OverlapCheck());  // using copy number for now, need to find a better way
0334       ZPos += (SMD_Step / 2.);
0335     }
0336     /* place the front plate */
0337     ZPos += (Plate_Step / 2.);
0338     new G4PVPlacement(PlateRotation, G4ThreeVector(0.0, Plate_YPos, Plate_ZPos + ZPos),
0339                       fb_plate_log,
0340                       G4String("front_plate"),
0341                       ZDC_envelope_log,
0342                       false, copyno_plate, OverlapCheck());
0343     ZPos += (Plate_Step / 2.);
0344     copyno_plate++;
0345     for (int j = 0; j < m_NLay; j++)
0346     {
0347       /* place the Absorber */
0348       ZPos += (Absorber_Step / 2.);
0349       new G4PVPlacement(PlateRotation, G4ThreeVector(0.0, 0.0, ZPos),
0350                         absorber_log,
0351                         G4String("single_absorber"),
0352                         ZDC_envelope_log,
0353                         false, i * 100 + j, OverlapCheck());
0354       ZPos += (Absorber_Step / 2.);
0355
0356       /* place the fiber plate */
0357       ZPos += (Gap_Step / 2.);
0358       int copyno = 27 * i + j;
0359       std::string name_fiber_plate = "Fiber_Plate_" + std::to_string(copyno);
0360       new G4PVPlacement(PlateRotation, G4ThreeVector(0.0, Gap_YPos, Gap_ZPos + ZPos),
0361                         fiber_plate_log,
0362                         name_fiber_plate,
0363                         ZDC_envelope_log,
0364                         false, copyno, OverlapCheck());
0365       ZPos += (Gap_Step / 2.);
0366     }
0367     /* place the back plate */
0368     ZPos += (Plate_Step / 2.);
0369     new G4PVPlacement(PlateRotation, G4ThreeVector(0.0, Plate_YPos, Plate_ZPos + ZPos),
0370                       fb_plate_log,
0371                       G4String("back_plate"),
0372                       ZDC_envelope_log,
0373                       false, copyno_plate, OverlapCheck());
0374     copyno_plate++;
0375     ZPos += (Plate_Step / 2.);
0376   }
0377
0378   /* Place envelope cone in simulation */
0379
0380   if (m_Layer == PHG4ZDCDefs::NORTH)
0381   {
0382     new G4PVPlacement(G4Transform3D(ZDC_rotm, G4ThreeVector(m_Params->get_double_param("place_x") * cm, m_Params->get_double_param("place_y") * cm, m_Params->get_double_param("place_z") * cm)),
0383                       ZDC_envelope_log, "ZDC_Envelope_North", logicWorld, false, PHG4ZDCDefs::NORTH, OverlapCheck());
0384   }
0385   else if (m_Layer == PHG4ZDCDefs::SOUTH)
0386   {
0387     ZDC_rotm.rotateY(180 * deg);
0388     new G4PVPlacement(G4Transform3D(ZDC_rotm, G4ThreeVector(m_Params->get_double_param("place_x") * cm, m_Params->get_double_param("place_y") * cm, -m_Params->get_double_param("place_z") * cm)),
0389                       ZDC_envelope_log, "ZDC_Envelope_South", logicWorld, false, PHG4ZDCDefs::SOUTH, OverlapCheck());
0390   }
0391   return;
0392 }