g4simulation/g4detectors/PHG4CylinderGeom_Spacalv3.cc

0001 // $$Id: PHG4CylinderGeom_Spacalv3.cc,v 1.3 2014/08/28 22:18:35 jinhuang Exp $$
0002
0003 /*!
0004  * \file ${file_name}
0005  * \brief
0006  * \author Jin Huang <jhuang@bnl.gov>
0007  * \version $$Revision: 1.3 $$
0008  * \date $$Date: 2014/08/28 22:18:35 $$
0009  */
0010
0011 #include "PHG4CylinderGeom_Spacalv3.h"
0012 #include "PHG4CylinderGeom_Spacalv1.h"  // for PHG4CylinderGeom_Spacalv1::...
0013
0014 #include <phparameter/PHParameters.h>
0015
0016 #include <Geant4/G4PhysicalConstants.hh>
0017
0018 #include <CLHEP/Units/SystemOfUnits.h>  // for twopi
0019
0020 #include <algorithm>
0021 #include <cassert>
0022 #include <cmath>
0023 #include <cstdlib>  // for exit
0024 #include <iostream>
0025 #include <limits>  // std::numeric_limits
0026 #include <map>
0027 #include <sstream>
0028
0029 PHG4CylinderGeom_Spacalv3::PHG4CylinderGeom_Spacalv3()
0030 {
0031   SetDefault();
0032 }
0033
0034 PHG4CylinderGeom_Spacalv3::~PHG4CylinderGeom_Spacalv3()
0035 {
0036   sector_tower_map.erase(sector_tower_map.begin(), sector_tower_map.end());
0037 }
0038
0039 void PHG4CylinderGeom_Spacalv3::identify(std::ostream& os) const
0040 {
0041   os << "PHG4CylinderGeom_Spacalv3: layer: " << layer
0042      //
0043      << ", radius: " << radius
0044      //
0045      << ", thickness: " << thickness
0046      //
0047      << ", zmin: " << zmin
0048      //
0049      << ", zmax: " << zmax << ", num scint: " << nscint << ", num sector: "
0050      << azimuthal_n_sec << ", unique tower: " << sector_tower_map.size()
0051      << std::endl;
0052 }
0053
0054 void PHG4CylinderGeom_Spacalv3::Print(Option_t* opt) const
0055 {
0056   PHG4CylinderGeom_Spacalv2::Print(opt);
0057
0058   std::cout << "\t"
0059             << "get_sidewall_outer_torr() = " << get_sidewall_outer_torr()
0060             << std::endl;
0061   std::cout << "\t"
0062             << "get_sidewall_thickness() = " << get_sidewall_thickness()
0063             << std::endl;
0064   std::cout << "\t"
0065             << "get_sidewall_mat() = " << get_sidewall_mat() << std::endl;
0066   std::cout << "\t"
0067             << "get_max_phi_bin_in_sec() = " << get_max_phi_bin_in_sec()
0068             << std::endl;
0069
0070   subtower_consistency_check();
0071   std::cout << "\t"
0072             << "get_n_subtower_eta() = " << get_n_subtower_eta() << std::endl;
0073   std::cout << "\t"
0074             << "get_n_subtower_phi() = " << get_n_subtower_phi() << std::endl;
0075
0076   std::cout << "\t"
0077             << "get_max_lightguide_height() = "
0078             << get_max_lightguide_height() << std::endl;
0079   std::cout << "\t"
0080             << "Containing " << sector_tower_map.size()
0081             << " unique towers per sector." << std::endl;
0082
0083   if (get_construction_verbose() >= 2)
0084   {
0085     for (const auto& it : sector_tower_map)
0086     {
0087       std::cout << "\t";
0088       std::cout << "\t";
0089       it.second.identify(std::cout);
0090     }
0091   }
0092 }
0093
0094 void PHG4CylinderGeom_Spacalv3::SetDefault()
0095 {
0096   PHG4CylinderGeom_Spacalv2::SetDefault();
0097
0098   //  radius = 90.000000;
0099   //  thickness = 26.130000;
0100   //  zmin = 149.470000;
0101   //  zmax = -zmin;
0102   //  azimuthal_n_sec = 32;
0103   //  polar_taper_ratio = 1.;
0104   //  assembly_spacing = 0.002500;
0105   sidewall_thickness = 0.075000;
0106   sidewall_outer_torr = 0.030000;
0107   sidewall_mat = "SS310";
0108   max_phi_bin_in_sec = 8;
0109   divider_mat = "G4_AIR";
0110   divider_width = 14.5;
0111 }
0112
0113 void PHG4CylinderGeom_Spacalv3::ImportParameters(const PHParameters& param)
0114 {
0115   PHG4CylinderGeom_Spacalv2::ImportParameters(param);
0116
0117   if (param.exist_double_param("sidewall_thickness"))
0118   {
0119     sidewall_thickness = param.get_double_param("sidewall_thickness");
0120   }
0121   if (param.exist_double_param("sidewall_outer_torr"))
0122   {
0123     sidewall_outer_torr = param.get_double_param("sidewall_outer_torr");
0124   }
0125   if (param.exist_string_param("sidewall_mat"))
0126   {
0127     sidewall_mat = param.get_string_param("sidewall_mat");
0128   }
0129   if (param.exist_int_param("max_phi_bin_in_sec"))
0130   {
0131     max_phi_bin_in_sec = param.get_int_param("max_phi_bin_in_sec");
0132   }
0133   if (param.exist_string_param("divider_mat"))
0134   {
0135     divider_mat = param.get_string_param("divider_mat");
0136   }
0137   if (param.exist_double_param("divider_width"))
0138   {
0139     divider_width = param.get_double_param("divider_width");
0140   }
0141
0142   // load sector_tower_map
0143   if (param.exist_int_param("sector_tower_map_size"))
0144   {
0145     sector_tower_map.clear();
0146
0147     const int n = param.get_int_param("sector_tower_map_size");
0148
0149     for (int i = 0; i < n; i++)
0150     {
0151       std::stringstream prefix;
0152       prefix << "sector_tower_map";
0153       prefix << "[" << i << "]"
0154              << ".";
0155
0156       geom_tower t;
0157       t.ImportParameters(param, prefix.str());
0158
0159       sector_tower_map[t.id] = t;
0160     }
0161   }
0162
0163   return;
0164 }
0165
0166 PHG4CylinderGeom_Spacalv3::geom_tower::geom_tower()
0167   : id(std::numeric_limits<int>::min())
0168   , pDz(std::numeric_limits<double>::quiet_NaN())
0169   , pDy1(std::numeric_limits<double>::quiet_NaN())
0170   , pDx1(std::numeric_limits<double>::quiet_NaN())
0171   , pDx2(std::numeric_limits<double>::quiet_NaN())
0172   , pDy2(std::numeric_limits<double>::quiet_NaN())
0173   , pDx3(std::numeric_limits<double>::quiet_NaN())
0174   , pDx4(std::numeric_limits<double>::quiet_NaN())
0175   , pTheta(std::numeric_limits<double>::quiet_NaN())
0176   , pPhi(std::numeric_limits<double>::quiet_NaN())
0177   , pAlp1(std::numeric_limits<double>::quiet_NaN())
0178   , pAlp2(std::numeric_limits<double>::quiet_NaN())
0179   , pRotationAngleX(std::numeric_limits<double>::quiet_NaN())
0180   , centralX(std::numeric_limits<double>::quiet_NaN())
0181   , centralY(std::numeric_limits<double>::quiet_NaN())
0182   , centralZ(std::numeric_limits<double>::quiet_NaN())
0183   , ModuleSkinThickness(std::numeric_limits<double>::quiet_NaN())
0184   , NFiberX(std::numeric_limits<int>::min())
0185   , NFiberY(std::numeric_limits<int>::min())
0186   , NSubtowerX(1)
0187   , NSubtowerY(1)
0188   , LightguideHeight(0)
0189   , LightguideTaperRatio(std::numeric_limits<double>::quiet_NaN())
0190   , LightguideMaterial("PMMA")
0191 {
0192 }
0193
0194 //! fiber layout -> fiber_id
0195 int PHG4CylinderGeom_Spacalv3::geom_tower::compose_fiber_id(int index_x,
0196                                                             int index_y) const
0197 {
0198   return NFiberY * index_x + index_y;
0199 }
0200
0201 //! fiber_id -> sub tower ID x: 0 ... NSubtowerX -1
0202 int PHG4CylinderGeom_Spacalv3::geom_tower::get_sub_tower_ID_x(int fiber_id) const
0203 {
0204   const int index_x = fiber_id / NFiberY;
0205   assert(index_x < NFiberX);
0206   assert(index_x >= 0);
0207
0208   const double sub_tower_width_x = (double) NFiberX / NSubtowerX;
0209   const int tower_ID_x = (NSubtowerX - 1) - floor(index_x / sub_tower_width_x);  //! x is negative azimuthal direction
0210   assert(tower_ID_x < NSubtowerX);
0211   assert(tower_ID_x >= 0);
0212
0213   return tower_ID_x;
0214 }
0215
0216 //! fiber_id -> sub tower ID y: 0 ... NSubtowerY -1
0217 int PHG4CylinderGeom_Spacalv3::geom_tower::get_sub_tower_ID_y(int fiber_id) const
0218 {
0219   assert(fiber_id >= 0);
0220   const int index_y = fiber_id % NFiberY;
0221
0222   const double sub_tower_width_y = (double) NFiberY / NSubtowerY;
0223
0224   assert(pRotationAngleX < 0);
0225   const int tower_ID_y = (NSubtowerY - 1) - floor(index_y / sub_tower_width_y);  //! y is negative polar direction
0226   assert(tower_ID_y < NSubtowerY);
0227   assert(tower_ID_y >= 0);
0228
0229   return tower_ID_y;
0230 }
0231
0232 double
0233 PHG4CylinderGeom_Spacalv3::geom_tower::get_position_fraction_x_in_sub_tower(int fiber_id) const
0234 {
0235   const int index_x = fiber_id / NFiberY;
0236   assert(index_x < NFiberX);
0237   assert(index_x >= 0);
0238
0239   const double sub_tower_width_x = (double) NFiberX / NSubtowerX;
0240   const double x_in_sub_tower = (fmod(index_x, sub_tower_width_x) + 0.5) / sub_tower_width_x;  //! x is negative azimuthal direction
0241   assert(x_in_sub_tower <= 1);
0242   assert(x_in_sub_tower >= 0);
0243
0244   return x_in_sub_tower;
0245 }
0246
0247 double
0248 PHG4CylinderGeom_Spacalv3::geom_tower::get_position_fraction_y_in_sub_tower(int fiber_id) const
0249 {
0250   assert(fiber_id >= 0);
0251   const int index_y = fiber_id % NFiberY;
0252
0253   const double sub_tower_width_y = (double) NFiberY / NSubtowerY;
0254
0255   assert(pRotationAngleX < 0);
0256   const double y_in_sub_tower = (fmod(index_y, sub_tower_width_y) + 0.5) / sub_tower_width_y;  //! y is negative polar direction
0257   assert(y_in_sub_tower <= 1);
0258   assert(y_in_sub_tower >= 0);
0259
0260   return y_in_sub_tower;
0261 }
0262
0263 void PHG4CylinderGeom_Spacalv3::geom_tower::identify(std::ostream& os) const
0264 {
0265   os << "PHG4CylinderGeom_Spacalv3::geom_super_tower"
0266      << "[" << id << "]"
0267      << " @ <Azimuthal, R, z> = " << centralX << ", " << centralY << ", "
0268      << centralZ << " cm"
0269      //
0270      << " with "
0271      //
0272      << "Half length = " << pDz << ", " << pDy1 << ", " << pDx1 << ", " << pDx2
0273      << ", " << pDy2 << ", " << pDx3 << ", " << pDx4 << ", "
0274      //
0275      << "Angles = " << pTheta << ", " << pPhi << ", " << pAlp1 << ", " << pAlp2
0276      << ", "                              //
0277      << "Rotation = " << pRotationAngleX  //
0278      << std::endl;
0279 }
0280
0281 void PHG4CylinderGeom_Spacalv3::geom_tower::ImportParameters(
0282     const PHParameters& param, const std::string& param_prefix)
0283 {
0284   id = param.get_int_param(param_prefix + "id");
0285   pDz = param.get_double_param(param_prefix + "pDz");
0286
0287   pDy1 = param.get_double_param(param_prefix + "pDy1");
0288   pDx1 = param.get_double_param(param_prefix + "pDx1");
0289   pDx2 = param.get_double_param(param_prefix + "pDx2");
0290   pDy2 = param.get_double_param(param_prefix + "pDy2");
0291   pDx3 = param.get_double_param(param_prefix + "pDx3");
0292   pDx4 = param.get_double_param(param_prefix + "pDx4");
0293
0294   pTheta = param.get_double_param(param_prefix + "pTheta");
0295   pPhi = param.get_double_param(param_prefix + "pPhi");
0296   pAlp1 = param.get_double_param(param_prefix + "pAlp1");
0297   pAlp2 = param.get_double_param(param_prefix + "pAlp2");
0298
0299   pRotationAngleX = param.get_double_param(param_prefix + "pRotationAngleX");
0300   centralX = param.get_double_param(param_prefix + "centralX");
0301   centralY = param.get_double_param(param_prefix + "centralY");
0302   centralZ = param.get_double_param(param_prefix + "centralZ");
0303
0304   ModuleSkinThickness = param.get_double_param(
0305       param_prefix + "ModuleSkinThickness");
0306   NFiberX = param.get_int_param(param_prefix + "NFiberX");
0307   NFiberY = param.get_int_param(param_prefix + "NFiberY");
0308   NSubtowerX = param.get_int_param(param_prefix + "NSubtowerX");
0309   NSubtowerY = param.get_int_param(param_prefix + "NSubtowerY");
0310
0311   LightguideHeight = param.get_double_param(param_prefix + "LightguideHeight");
0312   LightguideTaperRatio = param.get_double_param(
0313       param_prefix + "LightguideTaperRatio");
0314   LightguideMaterial = param.get_string_param(
0315       param_prefix + "LightguideMaterial");
0316 }
0317
0318 PHG4CylinderGeom_Spacalv3::scint_id_coder::scint_id_coder(int scint_id)
0319   : scint_ID(scint_id)
0320 {
0321   // NOLINTNEXTLINE(hicpp-signed-bitwise)
0322   sector_ID = (scint_ID >> (kfiber_bit + ktower_bit)) & ((1 << ksector_bit) - 1);
0323   // NOLINTNEXTLINE(hicpp-signed-bitwise)
0324   tower_ID = (scint_ID >> kfiber_bit) & ((1 << ktower_bit) - 1);
0325   // NOLINTNEXTLINE(hicpp-signed-bitwise)
0326   fiber_ID = (scint_ID) & ((1 << kfiber_bit) - 1);
0327 }
0328
0329 PHG4CylinderGeom_Spacalv3::scint_id_coder::scint_id_coder(int sector_id,
0330                                                           int tower_id, int fiber_id)
0331   : sector_ID(sector_id)
0332   , tower_ID(tower_id)
0333   , fiber_ID(fiber_id)
0334 {
0335   // NOLINTNEXTLINE(hicpp-signed-bitwise)
0336   assert(fiber_ID < (1 << kfiber_bit) && fiber_ID >= 0);
0337   // NOLINTNEXTLINE(hicpp-signed-bitwise)
0338   assert(tower_ID < (1 << ktower_bit) && tower_ID >= 0);
0339   // NOLINTNEXTLINE(hicpp-signed-bitwise)
0340   assert(sector_ID < (1 << ksector_bit) && sector_ID >= 0);
0341
0342   // NOLINTNEXTLINE(hicpp-signed-bitwise)
0343   scint_ID = (((sector_ID << ktower_bit) | tower_ID) << kfiber_bit) | fiber_ID;
0344 }
0345
0346 std::pair<int, int>
0347 PHG4CylinderGeom_Spacalv3::get_tower_z_phi_ID(const int tower_ID,
0348                                               const int sector_ID) const
0349 {
0350   // tower_ID to eta/z within a sector
0351   // NOLINTNEXTLINE(bugprone-integer-division)
0352   int z_bin = floor(tower_ID / 10);
0353
0354   int phi_bin_in_sec = -1;
0355
0356   if (get_config() == kFullProjective_2DTaper_SameLengthFiberPerTower or get_config() == kFullProjective_2DTaper)
0357   {
0358     // colume ID is from -x to +x at the top of the detector, which is reverse of the phi bin direction.
0359     phi_bin_in_sec = max_phi_bin_in_sec - (tower_ID % 10);
0360   }
0361   else if (get_config() == kFullProjective_2DTaper_Tilted_SameLengthFiberPerTower or get_config() == kFullProjective_2DTaper_Tilted)
0362   {
0363     phi_bin_in_sec = (tower_ID % 10);
0364   }
0365
0366   if (!(phi_bin_in_sec < max_phi_bin_in_sec and phi_bin_in_sec >= 0))
0367   {
0368     std::cout
0369         << "PHG4CylinderGeom_Spacalv3::get_tower_z_phi_ID - Fatal Error - invalid in put with "
0370         << "tower_ID = " << tower_ID << ", sector_ID = " << sector_ID << ", phi_bin_in_sec = " << phi_bin_in_sec
0371         << ". Dump object:" << std::endl;
0372     Print();
0373   }
0374
0375   assert(phi_bin_in_sec < max_phi_bin_in_sec and phi_bin_in_sec >= 0);
0376
0377   int phi_bin = sector_ID * max_phi_bin_in_sec + phi_bin_in_sec;
0378
0379   return std::make_pair(z_bin, phi_bin);
0380 }
0381
0382 //! get approximate radial position of tower
0383 double PHG4CylinderGeom_Spacalv3::
0384     get_tower_radial_position(const PHG4CylinderGeom_Spacalv3::geom_tower& tower) const
0385 {
0386   if (get_config() == kFullProjective_2DTaper_SameLengthFiberPerTower or get_config() == kFullProjective_2DTaper)
0387   {
0388     return tower.centralY;
0389   }
0390   else if (get_config() == kFullProjective_2DTaper_Tilted_SameLengthFiberPerTower or get_config() == kFullProjective_2DTaper_Tilted)
0391   {
0392     const double outter_wall_shift = get_sidewall_thickness() + get_sidewall_outer_torr() + get_assembly_spacing();
0393     assert(outter_wall_shift >= 0);
0394     const double tilted_radial_shift = outter_wall_shift / sin(M_PI / get_azimuthal_n_sec());
0395     assert(tilted_radial_shift >= 0);
0396     const double tower_radial =  //
0397         tilted_radial_shift * cos(get_azimuthal_tilt()) +
0398         get_half_radius() * sin(get_azimuthal_tilt()) * sin(get_azimuthal_tilt()) +
0399         tower.centralY * cos(get_azimuthal_tilt());
0400
0401     if (get_construction_verbose() >= 2)
0402     {
0403       std::cout
0404           << "PHG4CylinderGeom_Spacalv3::get_tower_radial_position - tower radial adjustment: "
0405              "from "
0406           << tower.centralY << " to " << tower_radial
0407           << std::endl;
0408     }
0409
0410     return tower_radial;
0411   }
0412   else
0413   {
0414     std::cout
0415         << "PHG4CylinderGeom_Spacalv3::get_tower_radial_position - ERROR - "
0416            "unsupported configuration!"
0417         << std::endl;
0418     Print();
0419     exit(10);
0420   }
0421   return NAN;
0422 }
0423
0424 //! check that all towers has consistent sub-tower divider
0425 void PHG4CylinderGeom_Spacalv3::subtower_consistency_check() const
0426 {
0427   if (sector_tower_map.begin() == sector_tower_map.end())
0428   {
0429     return;
0430   }
0431
0432   for (tower_map_t::const_iterator it = sector_tower_map.begin();
0433        it != sector_tower_map.end(); ++it)
0434   {
0435     assert(get_n_subtower_eta() == it->second.NSubtowerY);
0436     assert(get_n_subtower_phi() == it->second.NSubtowerX);
0437   }
0438
0439   if (get_construction_verbose())
0440   {
0441     std::cout
0442         << "PHG4CylinderGeom_Spacalv3::subtower_consistency_check - Passed with get_n_subtower_phi() = "
0443         << get_n_subtower_phi() << " and get_n_subtower_eta()"
0444         << get_n_subtower_eta() << std::endl;
0445   }
0446 }
0447 //! sub-tower divider along the polar direction
0448 int PHG4CylinderGeom_Spacalv3::get_n_subtower_eta() const
0449 {
0450   if (sector_tower_map.begin() == sector_tower_map.end())
0451   {
0452     return 0;
0453   }
0454   assert(sector_tower_map.begin() != sector_tower_map.end());
0455   return sector_tower_map.begin()->second.NSubtowerY;
0456 }
0457 //! sub-tower divider along the azimuthal direction
0458 int PHG4CylinderGeom_Spacalv3::get_n_subtower_phi() const
0459 {
0460   if (sector_tower_map.begin() == sector_tower_map.end())
0461   {
0462     return 0;
0463   }
0464   assert(sector_tower_map.begin() != sector_tower_map.end());
0465   return sector_tower_map.begin()->second.NSubtowerX;
0466 }
0467
0468 double
0469 PHG4CylinderGeom_Spacalv3::get_max_lightguide_height() const
0470 {
0471   double max_height = 0;
0472
0473   for (const auto& it : sector_tower_map)
0474   {
0475     const double h = it.second.LightguideHeight;
0476     max_height = std::max(max_height, h);
0477   }
0478
0479   return max_height;
0480 }
0481
0482 void PHG4CylinderGeom_Spacalv3::load_demo_sector_tower_map1()
0483 {
0484   std::cout << "PHG4CylinderGeom_Spacalv3::load_demo_sector_tower_map1 - "
0485             << "load four example central towers" << std::endl;
0486
0487   // Chris Cullen 2D spacal design July 2015
0488   radius = 90.000000;
0489   thickness = 26.130000;
0490   zmin = 149.470000;
0491   zmax = -zmin;
0492   azimuthal_n_sec = 32;
0493   max_phi_bin_in_sec = 8;
0494   sector_map.clear();
0495   sector_map[0] = 0;  // only install one sector
0496
0497   azimuthal_tilt = 0;
0498   azimuthal_seg_visible = false;
0499   polar_taper_ratio = 1.;
0500   assembly_spacing = 0.002500;
0501   sidewall_thickness = 0.075000;
0502   sidewall_outer_torr = 0.030000;
0503   sector_tower_map.clear();
0504
0505   {
0506     // tower 1023 based Row/Col = 102/3
0507     geom_tower geom;
0508     geom.id = 1023;
0509     geom.pDz = 6.751948;
0510     geom.pTheta = 0.038660;
0511     geom.pPhi = -2.829992;
0512     geom.pAlp1 = -0.000872;
0513     geom.pAlp2 = -0.000872;
0514     geom.pDy1 = 1.121195;
0515     geom.pDx1 = 1.191868;
0516     geom.pDx2 = 1.192521;
0517     geom.pDy2 = 1.281451;
0518     geom.pDx3 = 1.357679;
0519     geom.pDx4 = 1.358425;
0520     geom.centralX = -3.829796;
0521     geom.centralY = 105.060369;
0522     geom.centralZ = 3.686651;
0523     geom.pRotationAngleX = -1.547057;
0524     geom.ModuleSkinThickness = 0.010000;
0525     geom.NFiberX = 30;
0526     geom.NFiberY = 48;
0527     sector_tower_map[geom.id] = geom;
0528   }
0529
0530   {
0531     // tower 1024 based Row/Col = 102/4
0532     geom_tower geom;
0533     geom.id = 1024;
0534     geom.pDz = 6.751948;
0535     geom.pTheta = 0.017060;
0536     geom.pPhi = -2.373142;
0537     geom.pAlp1 = -0.000290;
0538     geom.pAlp2 = -0.000290;
0539     geom.pDy1 = 1.121195;
0540     geom.pDx1 = 1.190432;
0541     geom.pDx2 = 1.191082;
0542     geom.pDy2 = 1.281451;
0543     geom.pDx3 = 1.356043;
0544     geom.pDx4 = 1.356787;
0545     geom.centralX = -1.276086;
0546     geom.centralY = 105.060369;
0547     geom.centralZ = 3.686651;
0548     geom.pRotationAngleX = -1.547057;
0549     geom.ModuleSkinThickness = 0.010000;
0550     geom.NFiberX = 30;
0551     geom.NFiberY = 48;
0552     sector_tower_map[geom.id] = geom;
0553   }
0554 }
0555
0556 void PHG4CylinderGeom_Spacalv3::load_demo_sector_tower_map2()
0557 {
0558   std::cout << "PHG4CylinderGeom_Spacalv3::load_demo_sector_tower_map2 - "
0559             << "load one row of example forward towers" << std::endl;
0560
0561   // Chris Cullen 2D spacal design July 2015
0562   radius = 90.000000;
0563   thickness = 26.130000;
0564   zmin = 149.470000;
0565   zmax = -zmin;
0566   azimuthal_n_sec = 32;
0567   max_phi_bin_in_sec = 8;
0568   sector_map.clear();
0569   sector_map[0] = 0;  // only install one sector
0570   azimuthal_tilt = 0;
0571   azimuthal_seg_visible = false;
0572   polar_taper_ratio = 1.;
0573   assembly_spacing = 0.002500;
0574   sidewall_thickness = 0.075000;
0575   sidewall_outer_torr = 0.030000;
0576   sector_tower_map.clear();
0577
0578   {
0579     // tower 541 based Row/Col = 54/1
0580     geom_tower geom;
0581     geom.id = 541;
0582     geom.pDz = 8.326697;
0583     geom.pTheta = 0.053578;
0584     geom.pPhi = -3.015910;
0585     geom.pAlp1 = 0.068143;
0586     geom.pAlp2 = 0.068127;
0587     geom.pDy1 = 1.116997;
0588     geom.pDx1 = 1.235953;
0589     geom.pDx2 = 1.214052;
0590     geom.pDy2 = 1.231781;
0591     geom.pDx3 = 1.363889;
0592     geom.pDx4 = 1.339743;
0593     geom.centralX = -9.019696;
0594     geom.centralY = 105.782078;
0595     geom.centralZ = -134.504208;
0596     geom.pRotationAngleX = -2.482741;
0597     geom.ModuleSkinThickness = 0.010000;
0598     geom.NFiberX = 30;
0599     geom.NFiberY = 48;
0600     sector_tower_map[geom.id] = geom;
0601   }
0602   {
0603     // tower 542 based Row/Col = 54/2
0604     geom_tower geom;
0605     geom.id = 542;
0606     geom.pDz = 8.326697;
0607     geom.pTheta = 0.038558;
0608     geom.pPhi = -2.966436;
0609     geom.pAlp1 = 0.048643;
0610     geom.pAlp2 = 0.048632;
0611     geom.pDy1 = 1.116997;
0612     geom.pDx1 = 1.234208;
0613     geom.pDx2 = 1.212398;
0614     geom.pDy2 = 1.231781;
0615     geom.pDx3 = 1.361965;
0616     geom.pDx4 = 1.337918;
0617     geom.centralX = -6.439665;
0618     geom.centralY = 105.782078;
0619     geom.centralZ = -134.504208;
0620     geom.pRotationAngleX = -2.482741;
0621     geom.ModuleSkinThickness = 0.010000;
0622     geom.NFiberX = 30;
0623     geom.NFiberY = 48;
0624     sector_tower_map[geom.id] = geom;
0625   }
0626   {
0627     // tower 543 based Row/Col = 54/3
0628     geom_tower geom;
0629     geom.id = 543;
0630     geom.pDz = 8.326697;
0631     geom.pTheta = 0.023752;
0632     geom.pPhi = -2.854692;
0633     geom.pAlp1 = 0.029174;
0634     geom.pAlp2 = 0.029167;
0635     geom.pDy1 = 1.116997;
0636     geom.pDx1 = 1.233047;
0637     geom.pDx2 = 1.211297;
0638     geom.pDy2 = 1.231781;
0639     geom.pDx3 = 1.360684;
0640     geom.pDx4 = 1.336703;
0641     geom.centralX = -3.862610;
0642     geom.centralY = 105.782078;
0643     geom.centralZ = -134.504208;
0644     geom.pRotationAngleX = -2.482741;
0645     geom.ModuleSkinThickness = 0.010000;
0646     geom.NFiberX = 30;
0647     geom.NFiberY = 48;
0648     sector_tower_map[geom.id] = geom;
0649   }
0650   {
0651     // tower 544 based Row/Col = 54/4
0652     geom_tower geom;
0653     geom.id = 544;
0654     geom.pDz = 8.326697;
0655     geom.pTheta = 0.010142;
0656     geom.pPhi = -2.416982;
0657     geom.pAlp1 = 0.009723;
0658     geom.pAlp2 = 0.009720;
0659     geom.pDy1 = 1.116997;
0660     geom.pDx1 = 1.232467;
0661     geom.pDx2 = 1.210746;
0662     geom.pDy2 = 1.231781;
0663     geom.pDx3 = 1.360044;
0664     geom.pDx4 = 1.336096;
0665     geom.centralX = -1.287339;
0666     geom.centralY = 105.782078;
0667     geom.centralZ = -134.504208;
0668     geom.pRotationAngleX = -2.482741;
0669     geom.ModuleSkinThickness = 0.010000;
0670     geom.NFiberX = 30;
0671     geom.NFiberY = 48;
0672     sector_tower_map[geom.id] = geom;
0673   }
0674   {
0675     // tower 545 based Row/Col = 54/5
0676     geom_tower geom;
0677     geom.id = 545;
0678     geom.pDz = 8.326697;
0679     geom.pTheta = 0.010142;
0680     geom.pPhi = -0.724610;
0681     geom.pAlp1 = -0.009723;
0682     geom.pAlp2 = -0.009720;
0683     geom.pDy1 = 1.116997;
0684     geom.pDx1 = 1.232467;
0685     geom.pDx2 = 1.210746;
0686     geom.pDy2 = 1.231781;
0687     geom.pDx3 = 1.360044;
0688     geom.pDx4 = 1.336096;
0689     geom.centralX = 1.287339;
0690     geom.centralY = 105.782078;
0691     geom.centralZ = -134.504208;
0692     geom.pRotationAngleX = -2.482741;
0693     geom.ModuleSkinThickness = 0.010000;
0694     geom.NFiberX = 30;
0695     geom.NFiberY = 48;
0696     sector_tower_map[geom.id] = geom;
0697   }
0698   {
0699     // tower 546 based Row/Col = 54/6
0700     geom_tower geom;
0701     geom.id = 546;
0702     geom.pDz = 8.326697;
0703     geom.pTheta = 0.023752;
0704     geom.pPhi = -0.286901;
0705     geom.pAlp1 = -0.029174;
0706     geom.pAlp2 = -0.029167;
0707     geom.pDy1 = 1.116997;
0708     geom.pDx1 = 1.233047;
0709     geom.pDx2 = 1.211297;
0710     geom.pDy2 = 1.231781;
0711     geom.pDx3 = 1.360684;
0712     geom.pDx4 = 1.336703;
0713     geom.centralX = 3.862610;
0714     geom.centralY = 105.782078;
0715     geom.centralZ = -134.504208;
0716     geom.pRotationAngleX = -2.482741;
0717     geom.ModuleSkinThickness = 0.010000;
0718     geom.NFiberX = 30;
0719     geom.NFiberY = 48;
0720     sector_tower_map[geom.id] = geom;
0721   }
0722   {
0723     // tower 547 based Row/Col = 54/7
0724     geom_tower geom;
0725     geom.id = 547;
0726     geom.pDz = 8.326697;
0727     geom.pTheta = 0.038558;
0728     geom.pPhi = -0.175156;
0729     geom.pAlp1 = -0.048643;
0730     geom.pAlp2 = -0.048632;
0731     geom.pDy1 = 1.116997;
0732     geom.pDx1 = 1.234208;
0733     geom.pDx2 = 1.212398;
0734     geom.pDy2 = 1.231781;
0735     geom.pDx3 = 1.361965;
0736     geom.pDx4 = 1.337918;
0737     geom.centralX = 6.439665;
0738     geom.centralY = 105.782078;
0739     geom.centralZ = -134.504208;
0740     geom.pRotationAngleX = -2.482741;
0741     geom.ModuleSkinThickness = 0.010000;
0742     geom.NFiberX = 30;
0743     geom.NFiberY = 48;
0744     sector_tower_map[geom.id] = geom;
0745   }
0746   {
0747     // tower 548 based Row/Col = 54/8
0748     geom_tower geom;
0749     geom.id = 548;
0750     geom.pDz = 8.326697;
0751     geom.pTheta = 0.053578;
0752     geom.pPhi = -0.125683;
0753     geom.pAlp1 = -0.068143;
0754     geom.pAlp2 = -0.068127;
0755     geom.pDy1 = 1.116997;
0756     geom.pDx1 = 1.235953;
0757     geom.pDx2 = 1.214052;
0758     geom.pDy2 = 1.231781;
0759     geom.pDx3 = 1.363889;
0760     geom.pDx4 = 1.339743;
0761     geom.centralX = 9.019696;
0762     geom.centralY = 105.782078;
0763     geom.centralZ = -134.504208;
0764     geom.pRotationAngleX = -2.482741;
0765     geom.ModuleSkinThickness = 0.010000;
0766     geom.NFiberX = 30;
0767     geom.NFiberY = 48;
0768     sector_tower_map[geom.id] = geom;
0769   }
0770 }
0771
0772 void PHG4CylinderGeom_Spacalv3::load_demo_sector_tower_map4()
0773 {
0774   std::cout << "PHG4CylinderGeom_Spacalv3::load_demo_sector_tower_map4 - "
0775             << "FermiLab test beam 2014. Need to move to calibration database"
0776             << std::endl;
0777
0778   // From Oleg's documents
0779
0780   const int ny = 3;
0781   const int nx = 3;
0782
0783   const double inch_to_cm = 2.54;
0784
0785   const double screen_size_y = 2.108 * inch_to_cm;
0786   const double screen_size_x = 1.049 * inch_to_cm;
0787
0788   const double z_screen_6_7 = 0.5 * (6.6 + 6.75) * inch_to_cm;
0789   const double angle_screen_6_7 = 64.78 / 180. * M_PI;
0790   const double nawrrow_width_x = screen_size_x * sin(angle_screen_6_7);
0791
0792   const double z_screen_1_2 = 0.5 * (1.35 + 1.6) * inch_to_cm;
0793   const double angle_screen_1_2 = 90.56 / 180. * M_PI;
0794   const double wide_width_x = screen_size_x * sin(angle_screen_1_2);
0795
0796   const double module_length = z_screen_6_7 - z_screen_1_2;
0797
0798   // NOLINTNEXTLINE(hicpp-static-assert,misc-static-assert)
0799   assert(module_length > 0);
0800
0801   // tapering, dxwidth/dlength
0802   const double tapering_ratio = (wide_width_x - nawrrow_width_x) / module_length;
0803   assert(tapering_ratio < 1);
0804   assert(tapering_ratio > 0);
0805
0806   fiber_clading_thickness = 0.003;
0807   fiber_core_diameter = 0.050 - fiber_clading_thickness * 2;
0808
0809   assembly_spacing = 0.002500;
0810
0811   radius = (nawrrow_width_x) / tapering_ratio;
0812   thickness = module_length * 1.5;  // keep a large torlerence space
0813   zmin = -(0.5 * ny * screen_size_y + 2 * assembly_spacing * (ny + 1));
0814   zmax = -zmin;
0815   azimuthal_n_sec = floor(2 * M_PI / atan(tapering_ratio));
0816   max_phi_bin_in_sec = 1;
0817
0818   const double nawrrow_width_x_construction = radius * 2 * tan(M_PI / azimuthal_n_sec) - 2 * assembly_spacing;
0819   const double wide_width_x_construction = (radius + module_length) * 2 * tan(M_PI / azimuthal_n_sec) - 2 * assembly_spacing;
0820
0821   std::cout << "PHG4CylinderGeom_Spacalv3::load_demo_sector_tower_map4 - "
0822
0823             << "Adjust wide end width by ratio of "
0824             << wide_width_x_construction / wide_width_x
0825             << " and narrow end by ratio of "
0826             << nawrrow_width_x_construction / nawrrow_width_x << std::endl;
0827
0828   sector_map.clear();
0829   for (int sec = 0; sec < nx; ++sec)
0830   {
0831     const double rot = twopi / (double) (get_azimuthal_n_sec()) * ((double) (sec) -1);
0832
0833     sector_map[sec] = rot;
0834   }
0835
0836   azimuthal_tilt = 0;
0837   azimuthal_seg_visible = true;
0838   polar_taper_ratio = 1.;
0839   sidewall_thickness = 0;
0840   sidewall_outer_torr = 0;
0841   sector_tower_map.clear();
0842
0843   for (int y = 0; y < ny; y++)
0844   {
0845     geom_tower geom;
0846     geom.id = y;
0847     geom.pDz = module_length / 2;
0848     geom.pTheta = 0.;
0849     geom.pPhi = 0.;
0850     geom.pAlp1 = 0.;
0851     geom.pAlp2 = 0.;
0852     geom.pDy1 = 0.5 * screen_size_y;
0853     geom.pDx1 = 0.5 * nawrrow_width_x_construction;
0854     geom.pDx2 = 0.5 * nawrrow_width_x_construction;
0855     geom.pDy2 = 0.5 * screen_size_y;
0856     geom.pDx3 = 0.5 * wide_width_x_construction;
0857     geom.pDx4 = 0.5 * wide_width_x_construction;
0858
0859     geom.centralX = 0.;
0860     geom.centralY = module_length * 0.5 + radius + assembly_spacing;
0861     geom.centralZ = screen_size_y * (y - 1);
0862
0863     geom.pRotationAngleX = -M_PI / 2.;
0864     geom.ModuleSkinThickness = 0.010000;
0865     geom.NFiberX = 30;
0866     geom.NFiberY = 26 * 2 * 2;
0867     sector_tower_map[geom.id] = geom;
0868   }
0869 }