Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:20:45 

0001 #include "PHG4SpacalPrototypeDetector.h"
0002 
0003 #include <g4detectors/PHG4CylinderGeomContainer.h>
0004 #include <g4detectors/PHG4CylinderGeom_Spacalv1.h>  // for PHG4CylinderGeom_Spacalv1:...
0005 
0006 #include <phparameter/PHParameters.h>
0007 
0008 #include <g4detectors/PHG4CylinderGeom_Spacalv3.h>  // for PHG4CylinderGeom_...
0009 #include <g4main/PHG4Detector.h>                    // for PHG4Detector
0010 #include <g4main/PHG4Utils.h>
0011 
0012 #include <phool/PHCompositeNode.h>
0013 #include <phool/PHIODataNode.h>
0014 #include <phool/PHNode.h>          // for PHNode
0015 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0016 #include <phool/PHObject.h>        // for PHObject
0017 #include <phool/getClass.h>
0018 
0019 #include <Geant4/G4Box.hh>
0020 #include <Geant4/G4DisplacedSolid.hh>
0021 #include <Geant4/G4LogicalVolume.hh>
0022 #include <Geant4/G4Material.hh>
0023 #include <Geant4/G4PVPlacement.hh>
0024 #include <Geant4/G4PhysicalConstants.hh>
0025 #include <Geant4/G4String.hh>  // for G4String
0026 #include <Geant4/G4SubtractionSolid.hh>
0027 #include <Geant4/G4SystemOfUnits.hh>
0028 #include <Geant4/G4ThreeVector.hh>  // for G4ThreeVector
0029 #include <Geant4/G4Transform3D.hh>  // for G4Transform3D, G4Translate3D
0030 #include <Geant4/G4Trap.hh>
0031 #include <Geant4/G4Tubs.hh>
0032 #include <Geant4/G4Types.hh>  // for G4double
0033 #include <Geant4/G4UserLimits.hh>
0034 #include <Geant4/G4Vector3D.hh>
0035 #include <Geant4/G4VisAttributes.hh>
0036 
0037 #include <boost/foreach.hpp>
0038 
0039 #include <algorithm>
0040 #include <cassert>
0041 #include <cmath>
0042 #include <iostream>  // for operator<<, basic_ostream
0043 #include <numeric>   // std::accumulate
0044 #include <sstream>
0045 #include <string>  // std::string, std::to_string
0046 #include <vector>  // for vector
0047 
0048 class PHG4CylinderGeom;
0049 
0050 using namespace std;
0051 
0052 //_______________________________________________________________
0053 //note this inactive thickness is ~1.5% of a radiation length
0054 PHG4SpacalPrototypeDetector::PHG4SpacalPrototypeDetector(PHG4Subsystem* subsys, PHCompositeNode* Node, PHParameters* parameters, const std::string& dnam)
0055   : PHG4Detector(subsys, Node, dnam)
0056   , construction_params(parameters)
0057   , cylinder_solid(nullptr)
0058   , cylinder_logic(nullptr)
0059   , cylinder_physi(nullptr)
0060   ,  //
0061   active(0)
0062   , absorberactive(0)
0063   ,  //
0064   step_limits(nullptr)
0065   , clading_step_limits(nullptr)
0066   , fiber_core_step_limits(nullptr)
0067   ,  //
0068   _geom(nullptr)
0069 {
0070 }
0071 
0072 PHG4SpacalPrototypeDetector::~PHG4SpacalPrototypeDetector(void)
0073 {
0074   // deleting nullptr pointers is allowed (results in NOOP)
0075   // so checking for not null before deleting is not needed
0076   delete step_limits;
0077   delete clading_step_limits;
0078   delete fiber_core_step_limits;
0079   delete _geom;
0080 }
0081 
0082 //_______________________________________________________________
0083 int PHG4SpacalPrototypeDetector::IsInCylinderActive(
0084     const G4VPhysicalVolume* volume)
0085 {
0086   //  cout << "checking detector" << endl;
0087   if (active && fiber_core_vol.find(volume) != fiber_core_vol.end())
0088   {
0089     //      return fiber_core_vol.find(volume)->second;
0090     return FIBER_CORE;
0091   }
0092   if (absorberactive)
0093   {
0094     if (fiber_vol.find(volume) != fiber_vol.end())
0095       return FIBER_CLADING;
0096 
0097     if (block_vol.find(volume) != block_vol.end())
0098       return ABSORBER;
0099 
0100     if (calo_vol.find(volume) != calo_vol.end())
0101       return SUPPORT;
0102   }
0103   return INACTIVE;
0104 }
0105 
0106 //_______________________________________________________________
0107 void PHG4SpacalPrototypeDetector::ConstructMe(G4LogicalVolume* logicWorld)
0108 {
0109   PHNodeIterator iter(topNode());
0110   PHCompositeNode* parNode = dynamic_cast<PHCompositeNode*>(iter.findFirst(
0111       "PHCompositeNode", "RUN"));
0112   assert(parNode);
0113 
0114   if (!_geom)
0115   {
0116     _geom = new SpacalGeom_t();
0117   }
0118   assert(_geom);
0119 
0120   _geom->set_config(
0121       SpacalGeom_t::kFullProjective_2DTaper_SameLengthFiberPerTower);
0122   //  _geom->load_demo_sector_tower_map4();
0123   //    _geom->set_construction_verbose(2);
0124   _geom->ImportParameters(*construction_params);
0125 
0126   _geom->subtower_consistency_check();
0127 
0128   //  step_limits = new G4UserLimits(_geom->get_calo_step_size() * cm);
0129   //
0130   //  clading_step_limits = new G4UserLimits(
0131   //      _geom->get_fiber_clading_step_size() * cm);
0132   step_limits = nullptr;
0133   clading_step_limits = nullptr;
0134 
0135   fiber_core_step_limits = new G4UserLimits(
0136       _geom->get_fiber_core_step_size() * cm);
0137 
0138   const G4double z_rotation = construction_params->get_double_param(
0139                                   "z_rotation_degree") *
0140                               degree;
0141   const G4double enclosure_x = construction_params->get_double_param(
0142                                    "enclosure_x") *
0143                                cm;
0144   const G4double enclosure_x_shift = construction_params->get_double_param(
0145                                          "enclosure_x_shift") *
0146                                      cm;
0147   const G4double enclosure_y = construction_params->get_double_param(
0148                                    "enclosure_y") *
0149                                cm;
0150   const G4double enclosure_z = construction_params->get_double_param(
0151                                    "enclosure_z") *
0152                                cm;
0153 
0154   const G4double box_x_shift = (_geom->get_radius() + 0.5 * _geom->get_thickness()) * cm + enclosure_x_shift;
0155   const G4double box_z_shift = 0.5 * (_geom->get_zmin() + _geom->get_zmax()) * cm;
0156 
0157   //  G4Tubs* _cylinder_solid = new G4Tubs(G4String(GetName()),
0158   //      _geom->get_radius() * cm, _geom->get_max_radius() * cm,
0159   //      _geom->get_length() * cm / 2.0, 0, twopi);
0160   //  G4VSolid* cylinder_solid = new G4Box(G4String(GetName()),
0161   //      _geom->get_thickness() * 1.1 * 0.5 * cm, //
0162   //      twopi / _geom->get_azimuthal_n_sec() * _geom->get_sector_map().size()
0163   //          * 0.5 * _geom->get_max_radius() * cm, //
0164   //      _geom->get_length() * cm / 2.0);
0165   G4VSolid* cylinder_solid = new G4Box(G4String(GetName()),
0166                                        enclosure_x * 0.5,  //
0167                                        enclosure_y * 0.5,  //
0168                                        enclosure_z * 0.5);
0169 
0170   cylinder_solid = new G4DisplacedSolid(G4String(GetName() + "_displaced"),
0171                                         cylinder_solid, 0,  //
0172                                         G4ThreeVector(box_x_shift, 0, box_z_shift));
0173 
0174   G4Material* cylinder_mat = G4Material::GetMaterial("G4_AIR");
0175   assert(cylinder_mat);
0176 
0177   cylinder_logic = new G4LogicalVolume(cylinder_solid, cylinder_mat,
0178                                        G4String(GetName()), 0, 0, 0);
0179   G4VisAttributes* VisAtt = new G4VisAttributes();
0180   PHG4Utils::SetColour(VisAtt, "W_Epoxy");
0181   VisAtt->SetVisibility(true);
0182   VisAtt->SetForceSolid(
0183       (not _geom->is_virualize_fiber()) and (not _geom->is_azimuthal_seg_visible()));
0184   cylinder_logic->SetVisAttributes(VisAtt);
0185 
0186   G4Transform3D cylinder_place(
0187       G4Translate3D(_geom->get_xpos() * cm, _geom->get_ypos() * cm,
0188                     _geom->get_zpos() * cm)  //
0189       * G4RotateZ3D(z_rotation)              //
0190       * G4Translate3D(
0191             -(_geom->get_radius() + 0.5 * _geom->get_thickness()) * cm, 0,
0192             0));
0193 
0194   cylinder_physi = new G4PVPlacement(cylinder_place, cylinder_logic,
0195                                      G4String(GetName()), logicWorld, false, 0, OverlapCheck());
0196 
0197   // install sectors
0198   std::pair<G4LogicalVolume*, G4Transform3D> psec = Construct_AzimuthalSeg();
0199   G4LogicalVolume* sec_logic = psec.first;
0200   const G4Transform3D& sec_trans = psec.second;
0201   BOOST_FOREACH (const SpacalGeom_t::sector_map_t::value_type& val, _geom->get_sector_map())
0202   {
0203     const int sec = val.first;
0204     const double rot = val.second;
0205 
0206     G4Transform3D sec_place = G4RotateZ3D(rot) * sec_trans;
0207 
0208     ostringstream name;
0209     name << GetName() << "_sec" << sec;
0210 
0211     G4PVPlacement* calo_phys = new G4PVPlacement(sec_place, sec_logic,
0212                                                  G4String(name.str()), cylinder_logic, false, sec,
0213                                                  OverlapCheck());
0214     calo_vol[calo_phys] = sec;
0215   }
0216   _geom->set_nscint(_geom->get_nscint() * _geom->get_sector_map().size());
0217 
0218   // install electronics
0219   const G4double electronics_thickness = construction_params->get_double_param(
0220                                              "electronics_thickness") *
0221                                          cm;
0222   const string electronics_material = construction_params->get_string_param(
0223       "electronics_material");
0224 
0225   G4VSolid* electronics_solid = new G4Box(G4String(GetName()),
0226                                           electronics_thickness / 2.0,  //
0227                                           sin(
0228                                               twopi / _geom->get_azimuthal_n_sec() * _geom->get_sector_map().size() / 2) *
0229                                               (_geom->get_radius() - _geom->get_max_lightguide_height()) * cm,  //
0230                                           _geom->get_length() * cm / 2.0);
0231 
0232   electronics_solid = new G4DisplacedSolid(G4String(GetName() + "_displaced"),
0233                                            electronics_solid,
0234                                            0,  //
0235                                            G4ThreeVector(
0236                                                cos(
0237                                                    twopi / _geom->get_azimuthal_n_sec() * _geom->get_sector_map().size() / 2) *
0238                                                        (_geom->get_radius() - _geom->get_max_lightguide_height()) * cm -
0239                                                    electronics_thickness,
0240                                                0, box_z_shift));
0241 
0242   G4Material* electronics_mat = G4Material::GetMaterial(electronics_material);
0243   assert(electronics_mat);
0244 
0245   G4LogicalVolume* electronics_logic = new G4LogicalVolume(electronics_solid,
0246                                                            electronics_mat, G4String(GetName()) + "_Electronics", 0, 0, 0);
0247   VisAtt = new G4VisAttributes();
0248   PHG4Utils::SetColour(VisAtt, electronics_material);
0249   VisAtt->SetVisibility(true);
0250   VisAtt->SetForceSolid(not _geom->is_virualize_fiber());
0251   electronics_logic->SetVisAttributes(VisAtt);
0252 
0253   new G4PVPlacement(G4Transform3D::Identity, electronics_logic,
0254                     G4String(GetName()) + "_Electronics", cylinder_logic, false, 0,
0255                     OverlapCheck());
0256 
0257   // install the enclosure box
0258   const G4double enclosure_thickness = construction_params->get_double_param(
0259                                            "enclosure_thickness") *
0260                                        cm;
0261   const string enclosure_material = construction_params->get_string_param(
0262       "enclosure_material");
0263 
0264   G4VSolid* outer_enclosur_solid = new G4Box(
0265       G4String(GetName()) + "_outer_enclosur_solid", enclosure_x * 0.5,  //
0266       enclosure_y * 0.5,                                                 //
0267       enclosure_z * 0.5);
0268   G4VSolid* inner_enclosur_solid = new G4Box(
0269       G4String(GetName()) + "_inner_enclosur_solid",
0270       enclosure_x * 0.5 - enclosure_thickness,  //
0271       enclosure_y * 0.5 - enclosure_thickness,  //
0272       enclosure_z * 0.5 - enclosure_thickness);
0273   G4VSolid* enclosure_solid = new G4SubtractionSolid(
0274       G4String(GetName()) + "_enclosure_solid",  //
0275       outer_enclosur_solid, inner_enclosur_solid);
0276   enclosure_solid = new G4DisplacedSolid(
0277       G4String(GetName() + "_enclosure_solid_displaced"), enclosure_solid, 0,  //
0278       G4ThreeVector(box_x_shift, 0, box_z_shift));
0279 
0280   G4Material* enclosure_mat = G4Material::GetMaterial(enclosure_material);
0281   assert(enclosure_mat);
0282 
0283   G4LogicalVolume* enclosure_logic = new G4LogicalVolume(enclosure_solid,
0284                                                          enclosure_mat, G4String(GetName()) + "_enclosure", 0, 0, 0);
0285   VisAtt = new G4VisAttributes();
0286   PHG4Utils::SetColour(VisAtt, enclosure_material);
0287   VisAtt->SetVisibility(true);
0288   VisAtt->SetForceSolid(not _geom->is_virualize_fiber());
0289   enclosure_logic->SetVisAttributes(VisAtt);
0290 
0291   new G4PVPlacement(G4Transform3D::Identity, enclosure_logic,
0292                     G4String(GetName()) + "_enclosure", cylinder_logic, false, 0,
0293                     OverlapCheck());
0294 
0295   if (active)
0296   {
0297     ostringstream geonode;
0298     if (superdetector != "NONE")
0299     {
0300       geonode << "CYLINDERGEOM_" << superdetector;
0301     }
0302     else
0303     {
0304       geonode << "CYLINDERGEOM_" << detector_type;
0305     }
0306     PHG4CylinderGeomContainer* geo = findNode::getClass<
0307         PHG4CylinderGeomContainer>(topNode(), geonode.str());
0308     if (!geo)
0309     {
0310       geo = new PHG4CylinderGeomContainer();
0311       PHNodeIterator iter(topNode());
0312       PHCompositeNode* runNode =
0313           dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode",
0314                                                         "RUN"));
0315       PHIODataNode<PHObject>* newNode = new PHIODataNode<PHObject>(geo,
0316                                                                    geonode.str(), "PHObject");
0317       runNode->addNode(newNode);
0318     }
0319     // here in the detector class we have internal units, convert to cm
0320     // before putting into the geom object
0321     PHG4CylinderGeom* mygeom = new SpacalGeom_t(*_geom);
0322     geo->AddLayerGeom(0, mygeom);
0323     if (_geom->get_construction_verbose() >= 1)
0324     {
0325       cout << "PHG4SpacalPrototypeDetector::Construct::" << GetName()
0326            << " - Print Layer Geometry:" << endl;
0327       geo->identify();
0328     }
0329   }
0330 
0331   if (absorberactive)
0332   {
0333     ostringstream geonode;
0334     if (superdetector != "NONE")
0335     {
0336       geonode << "CYLINDERGEOM_ABSORBER_" << superdetector;
0337     }
0338     else
0339     {
0340       geonode << "CYLINDERGEOM_ABSORBER_" << detector_type << "_" << 0;
0341     }
0342     PHG4CylinderGeomContainer* geo = findNode::getClass<PHG4CylinderGeomContainer>(topNode(), geonode.str());
0343     if (!geo)
0344     {
0345       geo = new PHG4CylinderGeomContainer();
0346       PHNodeIterator iter(topNode());
0347       PHCompositeNode* runNode =
0348           dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode",
0349                                                         "RUN"));
0350       PHIODataNode<PHObject>* newNode = new PHIODataNode<PHObject>(geo, geonode.str(), "PHObject");
0351       runNode->addNode(newNode);
0352     }
0353     // here in the detector class we have internal units, convert to cm
0354     // before putting into the geom object
0355     PHG4CylinderGeom* mygeom = new SpacalGeom_t(*_geom);
0356     geo->AddLayerGeom(0, mygeom);
0357     if (_geom->get_construction_verbose() >= 1)
0358     {
0359       cout << "PHG4SpacalPrototypeDetector::Construct::" << GetName()
0360            << " - Print Absorber Layer Geometry:" << endl;
0361       geo->identify();
0362     }
0363   }
0364 
0365   if (_geom->get_construction_verbose() >= 1)
0366   {
0367     cout << "PHG4SpacalPrototypeDetector::Construct::" << GetName()
0368          << " - Completed. Print G4 Geometry:" << endl;
0369     Print();
0370     cout << "box_x_shift = " << box_x_shift << endl;
0371     cout << "box_z_shift = " << box_z_shift << endl;
0372   }
0373 }
0374 
0375 std::pair<G4LogicalVolume*, G4Transform3D>
0376 PHG4SpacalPrototypeDetector::Construct_AzimuthalSeg()
0377 {
0378   assert(_geom);
0379   assert(_geom->get_azimuthal_n_sec() > 4);
0380 
0381   G4VSolid* sec_solid = new G4Tubs(G4String(GetName() + string("_sec")),
0382                                    (_geom->get_radius() - _geom->get_max_lightguide_height()) * cm,
0383                                    _geom->get_max_radius() * cm, _geom->get_length() * cm / 2.0,
0384                                    halfpi - pi / _geom->get_azimuthal_n_sec(),
0385                                    twopi / _geom->get_azimuthal_n_sec());
0386 
0387   const G4double box_z_shift = 0.5 * (_geom->get_zmin() + _geom->get_zmax()) * cm;
0388   sec_solid = new G4DisplacedSolid(G4String(GetName() + "_displaced" + string("_sec")),
0389                                    sec_solid, 0,  //
0390                                    G4ThreeVector(0, 0, box_z_shift));
0391 
0392   G4Material* cylinder_mat = G4Material::GetMaterial("G4_AIR");
0393   assert(cylinder_mat);
0394 
0395   G4LogicalVolume* sec_logic = new G4LogicalVolume(sec_solid, cylinder_mat,
0396                                                    G4String(G4String(GetName() + string("_sec"))), 0, 0, step_limits);
0397 
0398   G4VisAttributes* VisAtt = new G4VisAttributes();
0399   VisAtt->SetColor(.5, .9, .5, .1);
0400   VisAtt->SetVisibility(
0401       _geom->is_azimuthal_seg_visible() or _geom->is_virualize_fiber());
0402   VisAtt->SetForceSolid(false);
0403   VisAtt->SetForceWireframe(true);
0404   sec_logic->SetVisAttributes(VisAtt);
0405 
0406   //  // construct walls
0407   //  G4Material * wall_mat = G4Material::GetMaterial(_geom->get_sidewall_mat());
0408   //  assert(wall_mat);
0409   //
0410   //  G4VisAttributes* wall_VisAtt = new G4VisAttributes();
0411   //  VisAtt->SetColor(.5, .9, .5, .5);
0412   //  wall_VisAtt->SetVisibility(
0413   //      _geom->is_azimuthal_seg_visible() and (not _geom->is_virualize_fiber()));
0414   //  wall_VisAtt->SetForceSolid(true);
0415   //
0416   assert(_geom->get_sidewall_thickness() == 0);
0417   //section skin not used in prototype construction
0418   //  if (_geom->get_sidewall_thickness() > 0)
0419   //    {
0420   //      // end walls
0421   //      if (_geom->get_construction_verbose() >= 1)
0422   //        {
0423   //          cout << "PHG4SpacalPrototypeDetector::Construct_AzimuthalSeg::"
0424   //              << GetName() << " - construct end walls." << endl;
0425   //        }
0426   //      G4Tubs* wall_solid = new G4Tubs(G4String(GetName() + string("_EndWall")),
0427   //          _geom->get_radius() * cm + _geom->get_sidewall_outer_torr() * cm,
0428   //          _geom->get_max_radius() * cm - _geom->get_sidewall_outer_torr() * cm,
0429   //          _geom->get_sidewall_thickness() * cm / 2.0,
0430   //          halfpi - pi / _geom->get_azimuthal_n_sec(),
0431   //          twopi / _geom->get_azimuthal_n_sec());
0432   //
0433   //      G4LogicalVolume * wall_logic = new G4LogicalVolume(wall_solid, wall_mat,
0434   //          G4String(G4String(GetName() + string("_EndWall"))), 0, 0,
0435   //          step_limits);
0436   //      wall_logic->SetVisAttributes(wall_VisAtt);
0437   //
0438   //      typedef map<int, double> z_locations_t;
0439   //      z_locations_t z_locations;
0440   //      z_locations[1000] = _geom->get_sidewall_thickness() * cm / 2.0
0441   //          + _geom->get_assembly_spacing() * cm;
0442   //      z_locations[1001] = _geom->get_length() * cm / 2.0
0443   //          - (_geom->get_sidewall_thickness() * cm / 2.0
0444   //              + _geom->get_assembly_spacing() * cm);
0445   //      z_locations[1100] = -(_geom->get_sidewall_thickness() * cm / 2.0
0446   //          + _geom->get_assembly_spacing() * cm);
0447   //      z_locations[1101] = -(_geom->get_length() * cm / 2.0
0448   //          - (_geom->get_sidewall_thickness() * cm / 2.0
0449   //              + _geom->get_assembly_spacing() * cm));
0450   //
0451   //      BOOST_FOREACH(z_locations_t::value_type& val, z_locations)
0452   //        {
0453   //          if (_geom->get_construction_verbose() >= 2)
0454   //            cout << "PHG4SpacalPrototypeDetector::Construct_AzimuthalSeg::"
0455   //                << GetName() << " - constructed End Wall ID " << val.first
0456   //                << " @ Z = " << val.second << endl;
0457   //
0458   //          G4Transform3D wall_trans = G4TranslateZ3D(val.second);
0459   //
0460   //          G4PVPlacement * wall_phys = new G4PVPlacement(wall_trans, wall_logic,
0461   //              G4String(GetName()) + G4String("_EndWall"), sec_logic,
0462   //              false, val.first, OverlapCheck());
0463   //
0464   //          calo_vol[wall_phys] = val.first;
0465   //        }
0466   //    }
0467   //
0468   //  if (_geom->get_sidewall_thickness() > 0)
0469   //    {
0470   //      // side walls
0471   //      if (_geom->get_construction_verbose() >= 1)
0472   //        {
0473   //          cout << "PHG4SpacalPrototypeDetector::Construct_AzimuthalSeg::"
0474   //              << GetName() << " - construct side walls." << endl;
0475   //        }
0476   //      G4Box* wall_solid = new G4Box(G4String(GetName() + string("_SideWall")),
0477   //          _geom->get_sidewall_thickness() * cm / 2.0,
0478   //          _geom->get_thickness() * cm / 2.
0479   //              - 2 * _geom->get_sidewall_outer_torr() * cm,
0480   //          (_geom->get_length() / 2.
0481   //              - 2
0482   //                  * (_geom->get_sidewall_thickness()
0483   //                      + 2. * _geom->get_assembly_spacing())) * cm * .5);
0484   //
0485   //      G4LogicalVolume * wall_logic = new G4LogicalVolume(wall_solid, wall_mat,
0486   //          G4String(G4String(GetName() + string("_SideWall"))), 0, 0,
0487   //          step_limits);
0488   //      wall_logic->SetVisAttributes(wall_VisAtt);
0489   //
0490   //      typedef map<int, pair<int, int> > sign_t;
0491   //      sign_t signs;
0492   //      signs[2000] = make_pair(+1, +1);
0493   //      signs[2001] = make_pair(+1, -1);
0494   //      signs[2100] = make_pair(-1, +1);
0495   //      signs[2101] = make_pair(-1, -1);
0496   //
0497   //      BOOST_FOREACH(sign_t::value_type& val, signs)
0498   //        {
0499   //          const int sign_z = val.second.first;
0500   //          const int sign_azimuth = val.second.second;
0501   //
0502   //          if (_geom->get_construction_verbose() >= 2)
0503   //            cout << "PHG4SpacalPrototypeDetector::Construct_AzimuthalSeg::"
0504   //                << GetName() << " - constructed Side Wall ID " << val.first
0505   //                << " with" << " Shift X = "
0506   //                << sign_azimuth
0507   //                    * (_geom->get_sidewall_thickness() * cm / 2.0
0508   //                        + _geom->get_sidewall_outer_torr() * cm)
0509   //                << " Rotation Z = "
0510   //                << sign_azimuth * pi / _geom->get_azimuthal_n_sec()
0511   //                << " Shift Z = " << sign_z * (_geom->get_length() * cm / 4)
0512   //                << endl;
0513   //
0514   //          G4Transform3D wall_trans = G4RotateZ3D(
0515   //              sign_azimuth * pi / _geom->get_azimuthal_n_sec())
0516   //              * G4TranslateZ3D(sign_z * (_geom->get_length() * cm / 4))
0517   //              * G4TranslateY3D(
0518   //                  _geom->get_radius() * cm + _geom->get_thickness() * cm / 2.)
0519   //              * G4TranslateX3D(
0520   //                  sign_azimuth
0521   //                      * (_geom->get_sidewall_thickness() * cm / 2.0
0522   //                          + _geom->get_sidewall_outer_torr() * cm));
0523   //
0524   //          G4PVPlacement * wall_phys = new G4PVPlacement(wall_trans, wall_logic,
0525   //              G4String(GetName()) + G4String("_EndWall"), sec_logic,
0526   //              false, val.first, OverlapCheck());
0527   //
0528   //          calo_vol[wall_phys] = val.first;
0529   //        }
0530   //    }
0531 
0532   // construct towers
0533 
0534   BOOST_FOREACH (const SpacalGeom_t::tower_map_t::value_type& val, _geom->get_sector_tower_map())
0535   {
0536     const SpacalGeom_t::geom_tower& g_tower = val.second;
0537     G4LogicalVolume* LV_tower = Construct_Tower(g_tower);
0538 
0539     G4Transform3D block_trans = G4TranslateX3D(g_tower.centralX * cm) * G4TranslateY3D(g_tower.centralY * cm) * G4TranslateZ3D(g_tower.centralZ * cm) * G4RotateX3D(g_tower.pRotationAngleX * rad);
0540 
0541     const bool overlapcheck_block = OverlapCheck() and (_geom->get_construction_verbose() >= 2);
0542 
0543     G4PVPlacement* block_phys = new G4PVPlacement(block_trans, LV_tower,
0544                                                   LV_tower->GetName(), sec_logic, false, g_tower.id,
0545                                                   overlapcheck_block);
0546     block_vol[block_phys] = g_tower.id;
0547 
0548     if (g_tower.LightguideHeight > 0)
0549     {
0550       // also build a light guide
0551 
0552       for (int ix = 0; ix < g_tower.NSubtowerX; ix++)
0553       //  int ix = 0;
0554       {
0555         for (int iy = 0; iy < g_tower.NSubtowerY; iy++)
0556         //        int iy = 0;
0557         {
0558           G4LogicalVolume* LV_lg = Construct_LightGuide(g_tower, ix,
0559                                                         iy);
0560 
0561           new G4PVPlacement(block_trans, LV_lg, LV_lg->GetName(),
0562                             sec_logic, false, g_tower.id, overlapcheck_block);
0563         }
0564       }
0565     }
0566   }
0567 
0568   cout << "PHG4SpacalPrototypeDetector::Construct_AzimuthalSeg::" << GetName()
0569        << " - constructed " << _geom->get_sector_tower_map().size()
0570        << " unique towers" << endl;
0571 
0572   return make_pair(sec_logic, G4Transform3D::Identity);
0573 }
0574 
0575 G4LogicalVolume*
0576 PHG4SpacalPrototypeDetector::Construct_Fiber(const G4double length,
0577                                              const string& id)
0578 {
0579   G4Tubs* fiber_solid = new G4Tubs(G4String(GetName() + string("_fiber") + id),
0580                                    0, _geom->get_fiber_outer_r() * cm, length / 2.0, 0, twopi);
0581 
0582   G4Material* clading_mat = G4Material::GetMaterial(
0583       _geom->get_fiber_clading_mat());
0584   assert(clading_mat);
0585 
0586   G4LogicalVolume* fiber_logic = new G4LogicalVolume(fiber_solid, clading_mat,
0587                                                      G4String(G4String(GetName() + string("_fiber") + id)), 0, 0,
0588                                                      clading_step_limits);
0589 
0590   {
0591     G4VisAttributes* VisAtt = new G4VisAttributes();
0592     PHG4Utils::SetColour(VisAtt, "G4_POLYSTYRENE");
0593     VisAtt->SetVisibility(_geom->is_virualize_fiber());
0594     VisAtt->SetForceSolid(_geom->is_virualize_fiber());
0595     fiber_logic->SetVisAttributes(VisAtt);
0596   }
0597 
0598   G4Tubs* core_solid = new G4Tubs(
0599       G4String(GetName() + string("_fiber_core") + id), 0,
0600       _geom->get_fiber_core_diameter() * cm / 2, length / 2.0, 0, twopi);
0601 
0602   G4Material* core_mat = G4Material::GetMaterial(_geom->get_fiber_core_mat());
0603   assert(core_mat);
0604 
0605   G4LogicalVolume* core_logic = new G4LogicalVolume(core_solid, core_mat,
0606                                                     G4String(G4String(GetName() + string("_fiber_core") + id)), 0, 0,
0607                                                     fiber_core_step_limits);
0608 
0609   {
0610     G4VisAttributes* VisAtt = new G4VisAttributes();
0611     PHG4Utils::SetColour(VisAtt, "G4_POLYSTYRENE");
0612     VisAtt->SetVisibility(false);
0613     VisAtt->SetForceSolid(false);
0614     core_logic->SetVisAttributes(VisAtt);
0615   }
0616 
0617   const bool overlapcheck_fiber = OverlapCheck() and (_geom->get_construction_verbose() >= 3);
0618   G4PVPlacement* core_physi = new G4PVPlacement(0, G4ThreeVector(), core_logic,
0619                                                 G4String(G4String(GetName() + string("_fiber_core") + id)), fiber_logic,
0620                                                 false, 0, overlapcheck_fiber);
0621   fiber_core_vol[core_physi] = 0;
0622 
0623   return fiber_logic;
0624 }
0625 
0626 void PHG4SpacalPrototypeDetector::Print(const std::string& /*what*/) const
0627 {
0628   assert(_geom);
0629 
0630   cout << "PHG4SpacalPrototypeDetector::Print::" << GetName()
0631        << " - Print Geometry:" << endl;
0632   _geom->Print();
0633 
0634   return;
0635 }
0636 
0637 //! Fully projective spacal with 2D tapered modules. To speed up construction, same-length fiber is used cross one tower
0638 int PHG4SpacalPrototypeDetector::Construct_Fibers_SameLengthFiberPerTower(
0639     const PHG4SpacalPrototypeDetector::SpacalGeom_t::geom_tower& g_tower,
0640     G4LogicalVolume* LV_tower)
0641 {
0642   assert(_geom);
0643 
0644   // construct fibers
0645 
0646   // first check out the fibers geometry
0647 
0648   typedef map<int, pair<G4Vector3D, G4Vector3D> > fiber_par_map;
0649   fiber_par_map fiber_par;
0650   G4double min_fiber_length = g_tower.pDz * cm * 4;
0651 
0652   G4Vector3D v_zshift = G4Vector3D(tan(g_tower.pTheta) * cos(g_tower.pPhi),
0653                                    tan(g_tower.pTheta) * sin(g_tower.pPhi), 1) *
0654                         g_tower.pDz;
0655   //  int fiber_ID = 0;
0656   for (int ix = 0; ix < g_tower.NFiberX; ix++)
0657   //  int ix = 0;
0658   {
0659     const double weighted_ix = static_cast<double>(ix) / (g_tower.NFiberX - 1.);
0660 
0661     const double weighted_pDx1 = (g_tower.pDx1 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0662     const double weighted_pDx2 = (g_tower.pDx2 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0663 
0664     const double weighted_pDx3 = (g_tower.pDx3 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0665     const double weighted_pDx4 = (g_tower.pDx4 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_ix * 2 - 1);
0666 
0667     for (int iy = 0; iy < g_tower.NFiberY; iy++)
0668     //        int iy = 0;
0669     {
0670       if ((ix + iy) % 2 == 1)
0671         continue;  // make a triangle pattern
0672 
0673       const double weighted_iy = static_cast<double>(iy) / (g_tower.NFiberY - 1.);
0674 
0675       const double weighted_pDy1 = (g_tower.pDy1 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
0676       const double weighted_pDy2 = (g_tower.pDy2 - g_tower.ModuleSkinThickness - _geom->get_fiber_outer_r()) * (weighted_iy * 2 - 1);
0677 
0678       const double weighted_pDx12 = weighted_pDx1 * (1 - weighted_iy) + weighted_pDx2 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp1);
0679       const double weighted_pDx34 = weighted_pDx3 * (1 - weighted_iy) + weighted_pDx4 * (weighted_iy) + weighted_pDy1 * tan(g_tower.pAlp2);
0680 
0681       G4Vector3D v1 = G4Vector3D(weighted_pDx12, weighted_pDy1, 0) - v_zshift;
0682       G4Vector3D v2 = G4Vector3D(weighted_pDx34, weighted_pDy2, 0) + v_zshift;
0683 
0684       G4Vector3D vector_fiber = (v2 - v1);
0685       vector_fiber *= (vector_fiber.mag() - _geom->get_fiber_outer_r()) / vector_fiber.mag();  // shrink by fiber boundary protection
0686       G4Vector3D center_fiber = (v2 + v1) / 2;
0687 
0688       // convert to Geant4 units
0689       vector_fiber *= cm;
0690       center_fiber *= cm;
0691 
0692       const int fiber_ID = g_tower.compose_fiber_id(ix, iy);
0693       fiber_par[fiber_ID] = make_pair(vector_fiber, center_fiber);
0694 
0695       const G4double fiber_length = vector_fiber.mag();
0696 
0697       min_fiber_length = min(fiber_length, min_fiber_length);
0698 
0699       //          ++fiber_ID;
0700     }
0701   }
0702 
0703   int fiber_count = 0;
0704 
0705   const G4double fiber_length = min_fiber_length;
0706   vector<G4double> fiber_cut;
0707 
0708   ostringstream ss;
0709   ss << string("_Tower") << g_tower.id;
0710   G4LogicalVolume* fiber_logic = Construct_Fiber(fiber_length, ss.str());
0711 
0712   BOOST_FOREACH (const fiber_par_map::value_type& val, fiber_par)
0713   {
0714     const int fiber_ID = val.first;
0715     G4Vector3D vector_fiber = val.second.first;
0716     G4Vector3D center_fiber = val.second.second;
0717     const G4double optimal_fiber_length = vector_fiber.mag();
0718 
0719     const G4Vector3D v1 = center_fiber - 0.5 * vector_fiber;
0720 
0721     // keep a statistics
0722     assert(optimal_fiber_length - fiber_length >= 0);
0723     fiber_cut.push_back(optimal_fiber_length - fiber_length);
0724 
0725     center_fiber += (fiber_length / optimal_fiber_length - 1) * 0.5 * vector_fiber;
0726     vector_fiber *= fiber_length / optimal_fiber_length;
0727 
0728     //      const G4Vector3D v1_new = center_fiber - 0.5 *vector_fiber;
0729 
0730     if (_geom->get_construction_verbose() >= 3)
0731       cout
0732           << "PHG4SpacalPrototypeDetector::Construct_Fibers_SameLengthFiberPerTower::"
0733           << GetName() << " - constructed fiber " << fiber_ID << ss.str()
0734           //
0735           << ", Length = " << optimal_fiber_length << "-"
0736           << (optimal_fiber_length - fiber_length) << "mm, "  //
0737           << "x = " << center_fiber.x() << "mm, "             //
0738           << "y = " << center_fiber.y() << "mm, "             //
0739           << "z = " << center_fiber.z() << "mm, "             //
0740           << "vx = " << vector_fiber.x() << "mm, "            //
0741           << "vy = " << vector_fiber.y() << "mm, "            //
0742           << "vz = " << vector_fiber.z() << "mm, "            //
0743           << endl;
0744 
0745     const G4double rotation_angle = G4Vector3D(0, 0, 1).angle(vector_fiber);
0746     const G4Vector3D rotation_axis =
0747         rotation_angle == 0 ? G4Vector3D(1, 0, 0) : G4Vector3D(0, 0, 1).cross(vector_fiber);
0748 
0749     G4Transform3D fiber_place(
0750         G4Translate3D(center_fiber.x(), center_fiber.y(), center_fiber.z()) * G4Rotate3D(rotation_angle, rotation_axis));
0751 
0752     ostringstream name;
0753     name << GetName() + string("_Tower") << g_tower.id << "_fiber"
0754          << ss.str();
0755 
0756     const bool overlapcheck_fiber = OverlapCheck() and (_geom->get_construction_verbose() >= 3);
0757     G4PVPlacement* fiber_physi = new G4PVPlacement(fiber_place, fiber_logic,
0758                                                    G4String(name.str()), LV_tower, false, fiber_ID,
0759                                                    overlapcheck_fiber);
0760     fiber_vol[fiber_physi] = fiber_ID;
0761 
0762     fiber_count++;
0763   }
0764 
0765   if (_geom->get_construction_verbose() >= 2)
0766     cout
0767         << "PHG4SpacalPrototypeDetector::Construct_Fibers_SameLengthFiberPerTower::"
0768         << GetName() << " - constructed tower ID " << g_tower.id << " with "
0769         << fiber_count << " fibers. Average fiber length cut = "
0770         << accumulate(fiber_cut.begin(), fiber_cut.end(), 0.0) / fiber_cut.size() << " mm" << endl;
0771 
0772   return fiber_count;
0773 }
0774 
0775 //! a block along z axis built with G4Trd that is slightly tapered in x dimension
0776 G4LogicalVolume*
0777 PHG4SpacalPrototypeDetector::Construct_Tower(
0778     const PHG4SpacalPrototypeDetector::SpacalGeom_t::geom_tower& g_tower)
0779 {
0780   assert(_geom);
0781 
0782   if (_geom->get_construction_verbose() >= 2)
0783   {
0784     cout << "PHG4SpacalPrototypeDetector::Construct_Tower::" << GetName()
0785          << " - constructed tower ID " << g_tower.id
0786          << " with geometry parameter: ";
0787     g_tower.identify(cout);
0788   }
0789 
0790   std::ostringstream sout;
0791   sout << "_" << g_tower.id;
0792   const G4String sTowerID(sout.str());
0793 
0794   //Processed PostionSeeds 1 from 1 1
0795 
0796   G4Trap* block_solid = new G4Trap(
0797       /*const G4String& pName*/ G4String(GetName()) + sTowerID,
0798       g_tower.pDz * cm,                                         // G4double pDz,
0799       g_tower.pTheta * rad, g_tower.pPhi * rad,                 // G4double pTheta, G4double pPhi,
0800       g_tower.pDy1 * cm, g_tower.pDx1 * cm, g_tower.pDx2 * cm,  // G4double pDy1, G4double pDx1, G4double pDx2,
0801       g_tower.pAlp1 * rad,                                      // G4double pAlp1,
0802       g_tower.pDy2 * cm, g_tower.pDx3 * cm, g_tower.pDx4 * cm,  // G4double pDy2, G4double pDx3, G4double pDx4,
0803       g_tower.pAlp2 * rad                                       // G4double pAlp2 //
0804   );
0805 
0806   G4Material* cylinder_mat = G4Material::GetMaterial(
0807       _geom->get_absorber_mat());
0808   assert(cylinder_mat);
0809 
0810   G4LogicalVolume* block_logic = new G4LogicalVolume(block_solid, cylinder_mat,
0811                                                      G4String(G4String(GetName()) + string("_Tower") + sTowerID), 0, 0,
0812                                                      step_limits);
0813 
0814   G4VisAttributes* VisAtt = new G4VisAttributes();
0815   //  PHG4Utils::SetColour(VisAtt, "W_Epoxy");
0816   VisAtt->SetColor(.3, .3, .3, .3);
0817   VisAtt->SetVisibility(
0818       _geom->is_azimuthal_seg_visible() or _geom->is_virualize_fiber());
0819   VisAtt->SetForceSolid(not _geom->is_virualize_fiber());
0820   block_logic->SetVisAttributes(VisAtt);
0821 
0822   // construct fibers
0823 
0824   int fiber_count = 0;
0825 
0826   fiber_count = Construct_Fibers_SameLengthFiberPerTower(g_tower, block_logic);
0827 
0828   if (_geom->get_construction_verbose() >= 2)
0829     cout << "PHG4SpacalPrototypeDetector::Construct_Tower::" << GetName()
0830          << " - constructed tower ID " << g_tower.id << " with " << fiber_count
0831          << " fibers using Construct_Fibers_SameLengthFiberPerTower" << endl;
0832 
0833   return block_logic;
0834 }
0835 
0836 G4LogicalVolume*
0837 PHG4SpacalPrototypeDetector::Construct_LightGuide(
0838     const PHG4SpacalPrototypeDetector::SpacalGeom_t::geom_tower& g_tower,
0839     const int index_x, const int index_y)
0840 {
0841   assert(_geom);
0842 
0843   std::ostringstream sout;
0844   sout << "_Lightguide_" << g_tower.id << "_" << index_x << "_" << index_y;
0845   const G4String sTowerID(sout.str());
0846 
0847   assert(g_tower.LightguideHeight > 0);
0848 
0849   // light guide parameters in PHENIX units
0850   const double weight_x1 = 1 - (double) index_y / g_tower.NSubtowerY;
0851   const double weight_x2 = 1 - (double) (index_y + 1) / g_tower.NSubtowerY;
0852   const double weight_xcenter = 1 - (double) (index_y + 0.5) / g_tower.NSubtowerY;
0853 
0854   assert(weight_x1 >= 0 and weight_x1 <= 1);
0855   assert(weight_x2 >= 0 and weight_x2 <= 1);
0856   assert(weight_xcenter >= 0 and weight_xcenter <= 1);
0857 
0858   const double lg_pDx1 = (g_tower.pDx1 * weight_x1  //
0859                           + g_tower.pDx2 * (1 - weight_x1)) /
0860                          g_tower.NSubtowerX;
0861   const double lg_pDx2 = (g_tower.pDx1 * weight_x2  //
0862                           + g_tower.pDx2 * (1 - weight_x2)) /
0863                          g_tower.NSubtowerX;
0864   const double lg_pDy1 = g_tower.pDy1 / g_tower.NSubtowerY;
0865   const double lg_Alp1 = atan(
0866       (g_tower.pDx2 - g_tower.pDx1) * (-g_tower.NSubtowerX + 1. + 2 * index_x) / (double) (g_tower.NSubtowerX) / (2. * g_tower.pDy1) + tan(g_tower.pAlp1));
0867 
0868   const double shift_xcenter = (g_tower.pDx1 * weight_xcenter           //
0869                                 + g_tower.pDx2 * (1 - weight_xcenter))  //
0870                                *                                        //
0871                                (-g_tower.NSubtowerX + 1. + 2 * index_x) / (double) (g_tower.NSubtowerX);
0872   const double shift_ycenter = g_tower.pDy1  //
0873                                *             //
0874                                (-g_tower.NSubtowerY + 1. + 2 * index_y) / (double) (g_tower.NSubtowerY);
0875 
0876   G4VSolid* block_solid = new G4Trap(
0877       /*const G4String& pName*/ G4String(GetName()) + sTowerID,
0878       0.5 * g_tower.LightguideHeight * cm,  // G4double pDz,
0879       0 * rad, 0 * rad,                     // G4double pTheta, G4double pPhi,
0880       g_tower.LightguideTaperRatio * lg_pDy1 * cm,
0881       g_tower.LightguideTaperRatio * lg_pDx1 * cm,
0882       g_tower.LightguideTaperRatio * lg_pDx2 * cm,  // G4double pDy1, G4double pDx1, G4double pDx2,
0883       lg_Alp1 * rad,                                // G4double pAlp1,
0884       lg_pDy1 * cm, lg_pDx1 * cm, lg_pDx2 * cm,     // G4double pDy2, G4double pDx3, G4double pDx4,
0885       lg_Alp1 * rad                                 // G4double pAlp2 //
0886   );
0887 
0888   block_solid = new G4DisplacedSolid(G4String(GetName() + "_displaced"),
0889                                      block_solid, 0,                                           //
0890                                      G4ThreeVector(                                            //
0891                                          tan(g_tower.pTheta * rad) * cos(g_tower.pPhi * rad),  //
0892                                          tan(g_tower.pTheta * rad) * sin(g_tower.pPhi * rad),  //
0893                                          1) *                                                  // G4ThreeVector
0894                                              -(g_tower.pDz) *
0895                                              cm                                                       //
0896                                          + G4ThreeVector(shift_xcenter * cm, shift_ycenter * cm, 0)   // shit in subtower direction
0897                                          + G4ThreeVector(0, 0, -0.5 * g_tower.LightguideHeight * cm)  //shift in the light guide height
0898   );
0899 
0900   G4Material* cylinder_mat = G4Material::GetMaterial(
0901       g_tower.LightguideMaterial);
0902   assert(cylinder_mat);
0903 
0904   G4LogicalVolume* block_logic = new G4LogicalVolume(block_solid, cylinder_mat,
0905                                                      G4String(G4String(GetName()) + string("_Tower") + sTowerID), 0, 0,
0906                                                      step_limits);
0907 
0908   G4VisAttributes* VisAtt = new G4VisAttributes();
0909   PHG4Utils::SetColour(VisAtt, g_tower.LightguideMaterial);
0910   //  VisAtt->SetColor(.3, .3, .3, .3);
0911   VisAtt->SetVisibility(
0912       _geom->is_azimuthal_seg_visible() or _geom->is_virualize_fiber());
0913   VisAtt->SetForceSolid(not _geom->is_virualize_fiber());
0914   block_logic->SetVisAttributes(VisAtt);
0915 
0916   return block_logic;
0917 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!