PHG4CylinderGeom_Spacalv3.cc

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4CylinderGeom_Spacalv3.cc

// $$Id: PHG4CylinderGeom_Spacalv3.cc,v 1.3 2014/08/28 22:18:35 jinhuang Exp $$

#include "PHG4CylinderGeom_Spacalv3.h"
#include "PHG4CylinderGeom_Spacalv1.h"  // for PHG4CylinderGeom_Spacalv1::...

#include <phparameter/PHParameters.h>

#include <Geant4/G4PhysicalConstants.hh>

#include <CLHEP/Units/SystemOfUnits.h>  // for twopi

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdlib>  // for exit
#include <iostream>
#include <limits>  // std::numeric_limits
#include <map>
#include <sstream>


PHG4CylinderGeom_Spacalv3::PHG4CylinderGeom_Spacalv3()
{
  SetDefault();
}

PHG4CylinderGeom_Spacalv3::~PHG4CylinderGeom_Spacalv3()
{
  sector_tower_map.erase(sector_tower_map.begin(), sector_tower_map.end());
}

void PHG4CylinderGeom_Spacalv3::identify(std::ostream& os) const
{
  os << "PHG4CylinderGeom_Spacalv3: layer: " << layer
     //
     << ", radius: " << radius
     //
     << ", thickness: " << thickness
     //
     << ", zmin: " << zmin
     //
     << ", zmax: " << zmax << ", num scint: " << nscint << ", num sector: "
     << azimuthal_n_sec << ", unique tower: " << sector_tower_map.size()
     << std::endl;
}

void PHG4CylinderGeom_Spacalv3::Print(Option_t* opt) const
{
  PHG4CylinderGeom_Spacalv2::Print(opt);

  std::cout << "\t"
       << "get_sidewall_outer_torr() = " << get_sidewall_outer_torr()
       << std::endl;
  std::cout << "\t"
       << "get_sidewall_thickness() = " << get_sidewall_thickness()
       << std::endl;
  std::cout << "\t"
       << "get_sidewall_mat() = " << get_sidewall_mat() << std::endl;
  std::cout << "\t"
       << "get_max_phi_bin_in_sec() = " << get_max_phi_bin_in_sec()
       << std::endl;

  subtower_consistency_check();
  std::cout << "\t"
       << "get_n_subtower_eta() = " << get_n_subtower_eta() << std::endl;
  std::cout << "\t"
       << "get_n_subtower_phi() = " << get_n_subtower_phi() << std::endl;

  std::cout << "\t"
       << "get_max_lightguide_height() = "
       << get_max_lightguide_height() << std::endl;
  std::cout << "\t"
       << "Containing " << sector_tower_map.size()
       << " unique towers per sector." << std::endl;

  if (get_construction_verbose() >= 2)
    for (const auto& it : sector_tower_map)
    {
      std::cout << "\t";
      std::cout << "\t";
      it.second.identify(std::cout);
    }
}

void PHG4CylinderGeom_Spacalv3::SetDefault()
{
  PHG4CylinderGeom_Spacalv2::SetDefault();

  //  radius = 90.000000;
  //  thickness = 26.130000;
  //  zmin = 149.470000;
  //  zmax = -zmin;
  //  azimuthal_n_sec = 32;
  //  polar_taper_ratio = 1.;
  //  assembly_spacing = 0.002500;
  sidewall_thickness = 0.075000;
  sidewall_outer_torr = 0.030000;
  sidewall_mat = "SS310";
  max_phi_bin_in_sec = 8;
  divider_mat = "G4_AIR";
  divider_width = 14.5;
}

void PHG4CylinderGeom_Spacalv3::ImportParameters(const PHParameters& param)
{
  PHG4CylinderGeom_Spacalv2::ImportParameters(param);

  if (param.exist_double_param("sidewall_thickness"))
    sidewall_thickness = param.get_double_param("sidewall_thickness");
  if (param.exist_double_param("sidewall_outer_torr"))
    sidewall_outer_torr = param.get_double_param("sidewall_outer_torr");
  if (param.exist_string_param("sidewall_mat"))
    sidewall_mat = param.get_string_param("sidewall_mat");
  if (param.exist_int_param("max_phi_bin_in_sec"))
    max_phi_bin_in_sec = param.get_int_param("max_phi_bin_in_sec");
  if (param.exist_string_param("divider_mat"))
    divider_mat = param.get_string_param("divider_mat");
  if (param.exist_double_param("divider_width"))
    divider_width = param.get_double_param("divider_width");

  // load sector_tower_map
  if (param.exist_int_param("sector_tower_map_size"))
  {
    sector_tower_map.clear();

    const int n = param.get_int_param("sector_tower_map_size");

    for (int i = 0; i < n; i++)
    {
      std::stringstream prefix;
      prefix << "sector_tower_map";
      prefix << "[" << i << "]"
             << ".";

      geom_tower t;
      t.ImportParameters(param, prefix.str());

      sector_tower_map[t.id] = t;
    }
  }

  return;
}

PHG4CylinderGeom_Spacalv3::geom_tower::geom_tower()
  : id(std::numeric_limits<int>::min())
  , pDz(std::numeric_limits<double>::quiet_NaN())
  , pDy1(std::numeric_limits<double>::quiet_NaN())
  , pDx1(std::numeric_limits<double>::quiet_NaN())
  , pDx2(std::numeric_limits<double>::quiet_NaN())
  , pDy2(std::numeric_limits<double>::quiet_NaN())
  , pDx3(std::numeric_limits<double>::quiet_NaN())
  , pDx4(std::numeric_limits<double>::quiet_NaN())
  , pTheta(std::numeric_limits<double>::quiet_NaN())
  , pPhi(std::numeric_limits<double>::quiet_NaN())
  , pAlp1(std::numeric_limits<double>::quiet_NaN())
  , pAlp2(std::numeric_limits<double>::quiet_NaN())
  , pRotationAngleX(std::numeric_limits<double>::quiet_NaN())
  , centralX(std::numeric_limits<double>::quiet_NaN())
  , centralY(std::numeric_limits<double>::quiet_NaN())
  , centralZ(std::numeric_limits<double>::quiet_NaN())
  , ModuleSkinThickness(std::numeric_limits<double>::quiet_NaN())
  , NFiberX(std::numeric_limits<int>::min())
  , NFiberY(std::numeric_limits<int>::min())
  , NSubtowerX(1)
  , NSubtowerY(1)
  , LightguideHeight(0)
  , LightguideTaperRatio(std::numeric_limits<double>::quiet_NaN())
  , LightguideMaterial("PMMA")
{
}

int PHG4CylinderGeom_Spacalv3::geom_tower::compose_fiber_id(int index_x,
                                                            int index_y) const
{
  return NFiberY * index_x + index_y;
}

int PHG4CylinderGeom_Spacalv3::geom_tower::get_sub_tower_ID_x(int fiber_id) const
{
  const int index_x = fiber_id / NFiberY;
  assert(index_x < NFiberX);
  assert(index_x >= 0);

  const double sub_tower_width_x = (double) NFiberX / NSubtowerX;
  const int tower_ID_x = (NSubtowerX - 1) - floor(index_x / sub_tower_width_x);  
  assert(tower_ID_x < NSubtowerX);
  assert(tower_ID_x >= 0);

  return tower_ID_x;
}

int PHG4CylinderGeom_Spacalv3::geom_tower::get_sub_tower_ID_y(int fiber_id) const
{
  assert(fiber_id >= 0);
  const int index_y = fiber_id % NFiberY;

  const double sub_tower_width_y = (double) NFiberY / NSubtowerY;

  assert(pRotationAngleX < 0);
  const int tower_ID_y = (NSubtowerY - 1) - floor(index_y / sub_tower_width_y);  
  assert(tower_ID_y < NSubtowerY);
  assert(tower_ID_y >= 0);

  return tower_ID_y;
}

double
PHG4CylinderGeom_Spacalv3::geom_tower::get_position_fraction_x_in_sub_tower(int fiber_id) const
{
  const int index_x = fiber_id / NFiberY;
  assert(index_x < NFiberX);
  assert(index_x >= 0);

  const double sub_tower_width_x = (double) NFiberX / NSubtowerX;
  const double x_in_sub_tower = (fmod(index_x, sub_tower_width_x) + 0.5) / sub_tower_width_x;  
  assert(x_in_sub_tower <= 1);
  assert(x_in_sub_tower >= 0);

  return x_in_sub_tower;
}

double
PHG4CylinderGeom_Spacalv3::geom_tower::get_position_fraction_y_in_sub_tower(int fiber_id) const
{
  assert(fiber_id >= 0);
  const int index_y = fiber_id % NFiberY;

  const double sub_tower_width_y = (double) NFiberY / NSubtowerY;

  assert(pRotationAngleX < 0);
  const double y_in_sub_tower = (fmod(index_y, sub_tower_width_y) + 0.5) / sub_tower_width_y;  
  assert(y_in_sub_tower <= 1);
  assert(y_in_sub_tower >= 0);

  return y_in_sub_tower;
}

void PHG4CylinderGeom_Spacalv3::geom_tower::identify(std::ostream& os) const
{
  os << "PHG4CylinderGeom_Spacalv3::geom_super_tower"
     << "[" << id << "]"
     << " @ <Azimuthal, R, z> = " << centralX << ", " << centralY << ", "
     << centralZ << " cm"
     //
     << " with "
     //
     << "Half length = " << pDz << ", " << pDy1 << ", " << pDx1 << ", " << pDx2
     << ", " << pDy2 << ", " << pDx3 << ", " << pDx4 << ", "
     //
     << "Angles = " << pTheta << ", " << pPhi << ", " << pAlp1 << ", " << pAlp2
     << ", "                              //
     << "Rotation = " << pRotationAngleX  //
     << std::endl;
}

void PHG4CylinderGeom_Spacalv3::geom_tower::ImportParameters(
    const PHParameters& param, const std::string& param_prefix)
{
  id = param.get_int_param(param_prefix + "id");
  pDz = param.get_double_param(param_prefix + "pDz");

  pDy1 = param.get_double_param(param_prefix + "pDy1");
  pDx1 = param.get_double_param(param_prefix + "pDx1");
  pDx2 = param.get_double_param(param_prefix + "pDx2");
  pDy2 = param.get_double_param(param_prefix + "pDy2");
  pDx3 = param.get_double_param(param_prefix + "pDx3");
  pDx4 = param.get_double_param(param_prefix + "pDx4");

  pTheta = param.get_double_param(param_prefix + "pTheta");
  pPhi = param.get_double_param(param_prefix + "pPhi");
  pAlp1 = param.get_double_param(param_prefix + "pAlp1");
  pAlp2 = param.get_double_param(param_prefix + "pAlp2");

  pRotationAngleX = param.get_double_param(param_prefix + "pRotationAngleX");
  centralX = param.get_double_param(param_prefix + "centralX");
  centralY = param.get_double_param(param_prefix + "centralY");
  centralZ = param.get_double_param(param_prefix + "centralZ");

  ModuleSkinThickness = param.get_double_param(
      param_prefix + "ModuleSkinThickness");
  NFiberX = param.get_int_param(param_prefix + "NFiberX");
  NFiberY = param.get_int_param(param_prefix + "NFiberY");
  NSubtowerX = param.get_int_param(param_prefix + "NSubtowerX");
  NSubtowerY = param.get_int_param(param_prefix + "NSubtowerY");

  LightguideHeight = param.get_double_param(param_prefix + "LightguideHeight");
  LightguideTaperRatio = param.get_double_param(
      param_prefix + "LightguideTaperRatio");
  LightguideMaterial = param.get_string_param(
      param_prefix + "LightguideMaterial");
}

PHG4CylinderGeom_Spacalv3::scint_id_coder::scint_id_coder(int scint_id)
  : scint_ID(scint_id)
{
  sector_ID = (scint_ID >> (kfiber_bit + ktower_bit)) & ((1 << ksector_bit) - 1);
  tower_ID = (scint_ID >> kfiber_bit) & ((1 << ktower_bit) - 1);
  fiber_ID = (scint_ID) & ((1 << kfiber_bit) - 1);
}

PHG4CylinderGeom_Spacalv3::scint_id_coder::scint_id_coder(int sector_id,
                                                          int tower_id, int fiber_id)
  : sector_ID(sector_id)
  , tower_ID(tower_id)
  , fiber_ID(fiber_id)
{
  assert(fiber_ID < (1 << kfiber_bit) and fiber_ID >= 0);
  assert(tower_ID < (1 << ktower_bit) and tower_ID >= 0);
  assert(sector_ID < (1 << ksector_bit) and sector_ID >= 0);

  scint_ID = (((sector_ID << ktower_bit) | tower_ID) << kfiber_bit) | fiber_ID;
}

std::pair<int, int>
PHG4CylinderGeom_Spacalv3::get_tower_z_phi_ID(const int tower_ID,
                                              const int sector_ID) const
{
  // tower_ID to eta/z within a sector
  int z_bin = floor(tower_ID / 10);

  int phi_bin_in_sec = -1;

  if (get_config() == kFullProjective_2DTaper_SameLengthFiberPerTower or get_config() == kFullProjective_2DTaper)
    // colume ID is from -x to +x at the top of the detector, which is reverse of the phi bin direction.
    phi_bin_in_sec = max_phi_bin_in_sec - (tower_ID % 10);
  else if (get_config() == kFullProjective_2DTaper_Tilted_SameLengthFiberPerTower or get_config() == kFullProjective_2DTaper_Tilted)
    phi_bin_in_sec = (tower_ID % 10);

  if (!(phi_bin_in_sec < max_phi_bin_in_sec and phi_bin_in_sec >= 0))
  {
    std::cout
        << "PHG4CylinderGeom_Spacalv3::get_tower_z_phi_ID - Fatal Error - invalid in put with "
        << "tower_ID = " << tower_ID << ", sector_ID = " << sector_ID << ", phi_bin_in_sec = " << phi_bin_in_sec
        << ". Dump object:" << std::endl;
    Print();
  }

  assert(phi_bin_in_sec < max_phi_bin_in_sec and phi_bin_in_sec >= 0);

  int phi_bin = sector_ID * max_phi_bin_in_sec + phi_bin_in_sec;

  return std::make_pair(z_bin, phi_bin);
}

double PHG4CylinderGeom_Spacalv3::
    get_tower_radial_position(const PHG4CylinderGeom_Spacalv3::geom_tower& tower) const
{
  if (get_config() == kFullProjective_2DTaper_SameLengthFiberPerTower or get_config() == kFullProjective_2DTaper)
    return tower.centralY;
  else if (get_config() == kFullProjective_2DTaper_Tilted_SameLengthFiberPerTower or get_config() == kFullProjective_2DTaper_Tilted)
  {
    const double outter_wall_shift = get_sidewall_thickness() + get_sidewall_outer_torr() + get_assembly_spacing();
    assert(outter_wall_shift >= 0);
    const double tilted_radial_shift = outter_wall_shift / sin(M_PI / get_azimuthal_n_sec());
    assert(tilted_radial_shift >= 0);
    const double tower_radial =  //
        tilted_radial_shift * cos(get_azimuthal_tilt()) +
        get_half_radius() * sin(get_azimuthal_tilt()) * sin(get_azimuthal_tilt()) +
        tower.centralY * cos(get_azimuthal_tilt());

    if (get_construction_verbose() >= 2)
    {
      std::cout
          << "PHG4CylinderGeom_Spacalv3::get_tower_radial_position - tower radial adjustment: "
             "from "
          << tower.centralY << " to " << tower_radial
          << std::endl;
    }

    return tower_radial;
  }
  else
  {
    std::cout
        << "PHG4CylinderGeom_Spacalv3::get_tower_radial_position - ERROR - "
           "unsupported configuration!"
        << std::endl;
    Print();
    exit(10);
  }
  return NAN;
}

void PHG4CylinderGeom_Spacalv3::subtower_consistency_check() const
{
  if (sector_tower_map.begin() == sector_tower_map.end()) return;

  for (tower_map_t::const_iterator it = sector_tower_map.begin();
       it != sector_tower_map.end(); ++it)
  {
    assert(get_n_subtower_eta() == it->second.NSubtowerY);
    assert(get_n_subtower_phi() == it->second.NSubtowerX);
  }

  if (get_construction_verbose())
  {
    std::cout
        << "PHG4CylinderGeom_Spacalv3::subtower_consistency_check - Passed with get_n_subtower_phi() = "
        << get_n_subtower_phi() << " and get_n_subtower_eta()"
        << get_n_subtower_eta() << std::endl;
  }
}
int PHG4CylinderGeom_Spacalv3::get_n_subtower_eta() const
{
  if (sector_tower_map.begin() == sector_tower_map.end()) return 0;
  assert(sector_tower_map.begin() != sector_tower_map.end());
  return sector_tower_map.begin()->second.NSubtowerY;
}
int PHG4CylinderGeom_Spacalv3::get_n_subtower_phi() const
{
  if (sector_tower_map.begin() == sector_tower_map.end()) return 0;
  assert(sector_tower_map.begin() != sector_tower_map.end());
  return sector_tower_map.begin()->second.NSubtowerX;
}

double
PHG4CylinderGeom_Spacalv3::get_max_lightguide_height() const
{
  double max_height = 0;

  for (const auto& it : sector_tower_map)
  {
    const double h = it.second.LightguideHeight;
    max_height = std::max(max_height, h);
  }

  return max_height;
}

void PHG4CylinderGeom_Spacalv3::load_demo_sector_tower_map1()
{
  std::cout << "PHG4CylinderGeom_Spacalv3::load_demo_sector_tower_map1 - "
       << "load four example central towers" << std::endl;

  // Chris Cullen 2D spacal design July 2015
  radius = 90.000000;
  thickness = 26.130000;
  zmin = 149.470000;
  zmax = -zmin;
  azimuthal_n_sec = 32;
  max_phi_bin_in_sec = 8;
  sector_map.clear();
  sector_map[0] = 0;  // only install one sector

  azimuthal_tilt = 0;
  azimuthal_seg_visible = false;
  polar_taper_ratio = 1.;
  assembly_spacing = 0.002500;
  sidewall_thickness = 0.075000;
  sidewall_outer_torr = 0.030000;
  sector_tower_map.clear();

  {
    // tower 1023 based Row/Col = 102/3
    geom_tower geom;
    geom.id = 1023;
    geom.pDz = 6.751948;
    geom.pTheta = 0.038660;
    geom.pPhi = -2.829992;
    geom.pAlp1 = -0.000872;
    geom.pAlp2 = -0.000872;
    geom.pDy1 = 1.121195;
    geom.pDx1 = 1.191868;
    geom.pDx2 = 1.192521;
    geom.pDy2 = 1.281451;
    geom.pDx3 = 1.357679;
    geom.pDx4 = 1.358425;
    geom.centralX = -3.829796;
    geom.centralY = 105.060369;
    geom.centralZ = 3.686651;
    geom.pRotationAngleX = -1.547057;
    geom.ModuleSkinThickness = 0.010000;
    geom.NFiberX = 30;
    geom.NFiberY = 48;
    sector_tower_map[geom.id] = geom;
  }

  {
    // tower 1024 based Row/Col = 102/4
    geom_tower geom;
    geom.id = 1024;
    geom.pDz = 6.751948;
    geom.pTheta = 0.017060;
    geom.pPhi = -2.373142;
    geom.pAlp1 = -0.000290;
    geom.pAlp2 = -0.000290;
    geom.pDy1 = 1.121195;
    geom.pDx1 = 1.190432;
    geom.pDx2 = 1.191082;
    geom.pDy2 = 1.281451;
    geom.pDx3 = 1.356043;
    geom.pDx4 = 1.356787;
    geom.centralX = -1.276086;
    geom.centralY = 105.060369;
    geom.centralZ = 3.686651;
    geom.pRotationAngleX = -1.547057;
    geom.ModuleSkinThickness = 0.010000;
    geom.NFiberX = 30;
    geom.NFiberY = 48;
    sector_tower_map[geom.id] = geom;
  }
}

void PHG4CylinderGeom_Spacalv3::load_demo_sector_tower_map2()
{
  std::cout << "PHG4CylinderGeom_Spacalv3::load_demo_sector_tower_map2 - "
       << "load one row of example forward towers" << std::endl;

  // Chris Cullen 2D spacal design July 2015
  radius = 90.000000;
  thickness = 26.130000;
  zmin = 149.470000;
  zmax = -zmin;
  azimuthal_n_sec = 32;
  max_phi_bin_in_sec = 8;
  sector_map.clear();
  sector_map[0] = 0;  // only install one sector
  azimuthal_tilt = 0;
  azimuthal_seg_visible = false;
  polar_taper_ratio = 1.;
  assembly_spacing = 0.002500;
  sidewall_thickness = 0.075000;
  sidewall_outer_torr = 0.030000;
  sector_tower_map.clear();

  {
    // tower 541 based Row/Col = 54/1
    geom_tower geom;
    geom.id = 541;
    geom.pDz = 8.326697;
    geom.pTheta = 0.053578;
    geom.pPhi = -3.015910;
    geom.pAlp1 = 0.068143;
    geom.pAlp2 = 0.068127;
    geom.pDy1 = 1.116997;
    geom.pDx1 = 1.235953;
    geom.pDx2 = 1.214052;
    geom.pDy2 = 1.231781;
    geom.pDx3 = 1.363889;
    geom.pDx4 = 1.339743;
    geom.centralX = -9.019696;
    geom.centralY = 105.782078;
    geom.centralZ = -134.504208;
    geom.pRotationAngleX = -2.482741;
    geom.ModuleSkinThickness = 0.010000;
    geom.NFiberX = 30;
    geom.NFiberY = 48;
    sector_tower_map[geom.id] = geom;
  }
  {
    // tower 542 based Row/Col = 54/2
    geom_tower geom;
    geom.id = 542;
    geom.pDz = 8.326697;
    geom.pTheta = 0.038558;
    geom.pPhi = -2.966436;
    geom.pAlp1 = 0.048643;
    geom.pAlp2 = 0.048632;
    geom.pDy1 = 1.116997;
    geom.pDx1 = 1.234208;
    geom.pDx2 = 1.212398;
    geom.pDy2 = 1.231781;
    geom.pDx3 = 1.361965;
    geom.pDx4 = 1.337918;
    geom.centralX = -6.439665;
    geom.centralY = 105.782078;
    geom.centralZ = -134.504208;
    geom.pRotationAngleX = -2.482741;
    geom.ModuleSkinThickness = 0.010000;
    geom.NFiberX = 30;
    geom.NFiberY = 48;
    sector_tower_map[geom.id] = geom;
  }
  {
    // tower 543 based Row/Col = 54/3
    geom_tower geom;
    geom.id = 543;
    geom.pDz = 8.326697;
    geom.pTheta = 0.023752;
    geom.pPhi = -2.854692;
    geom.pAlp1 = 0.029174;
    geom.pAlp2 = 0.029167;
    geom.pDy1 = 1.116997;
    geom.pDx1 = 1.233047;
    geom.pDx2 = 1.211297;
    geom.pDy2 = 1.231781;
    geom.pDx3 = 1.360684;
    geom.pDx4 = 1.336703;
    geom.centralX = -3.862610;
    geom.centralY = 105.782078;
    geom.centralZ = -134.504208;
    geom.pRotationAngleX = -2.482741;
    geom.ModuleSkinThickness = 0.010000;
    geom.NFiberX = 30;
    geom.NFiberY = 48;
    sector_tower_map[geom.id] = geom;
  }
  {
    // tower 544 based Row/Col = 54/4
    geom_tower geom;
    geom.id = 544;
    geom.pDz = 8.326697;
    geom.pTheta = 0.010142;
    geom.pPhi = -2.416982;
    geom.pAlp1 = 0.009723;
    geom.pAlp2 = 0.009720;
    geom.pDy1 = 1.116997;
    geom.pDx1 = 1.232467;
    geom.pDx2 = 1.210746;
    geom.pDy2 = 1.231781;
    geom.pDx3 = 1.360044;
    geom.pDx4 = 1.336096;
    geom.centralX = -1.287339;
    geom.centralY = 105.782078;
    geom.centralZ = -134.504208;
    geom.pRotationAngleX = -2.482741;
    geom.ModuleSkinThickness = 0.010000;
    geom.NFiberX = 30;
    geom.NFiberY = 48;
    sector_tower_map[geom.id] = geom;
  }
  {
    // tower 545 based Row/Col = 54/5
    geom_tower geom;
    geom.id = 545;
    geom.pDz = 8.326697;
    geom.pTheta = 0.010142;
    geom.pPhi = -0.724610;
    geom.pAlp1 = -0.009723;
    geom.pAlp2 = -0.009720;
    geom.pDy1 = 1.116997;
    geom.pDx1 = 1.232467;
    geom.pDx2 = 1.210746;
    geom.pDy2 = 1.231781;
    geom.pDx3 = 1.360044;
    geom.pDx4 = 1.336096;
    geom.centralX = 1.287339;
    geom.centralY = 105.782078;
    geom.centralZ = -134.504208;
    geom.pRotationAngleX = -2.482741;
    geom.ModuleSkinThickness = 0.010000;
    geom.NFiberX = 30;
    geom.NFiberY = 48;
    sector_tower_map[geom.id] = geom;
  }
  {
    // tower 546 based Row/Col = 54/6
    geom_tower geom;
    geom.id = 546;
    geom.pDz = 8.326697;
    geom.pTheta = 0.023752;
    geom.pPhi = -0.286901;
    geom.pAlp1 = -0.029174;
    geom.pAlp2 = -0.029167;
    geom.pDy1 = 1.116997;
    geom.pDx1 = 1.233047;
    geom.pDx2 = 1.211297;
    geom.pDy2 = 1.231781;
    geom.pDx3 = 1.360684;
    geom.pDx4 = 1.336703;
    geom.centralX = 3.862610;
    geom.centralY = 105.782078;
    geom.centralZ = -134.504208;
    geom.pRotationAngleX = -2.482741;
    geom.ModuleSkinThickness = 0.010000;
    geom.NFiberX = 30;
    geom.NFiberY = 48;
    sector_tower_map[geom.id] = geom;
  }
  {
    // tower 547 based Row/Col = 54/7
    geom_tower geom;
    geom.id = 547;
    geom.pDz = 8.326697;
    geom.pTheta = 0.038558;
    geom.pPhi = -0.175156;
    geom.pAlp1 = -0.048643;
    geom.pAlp2 = -0.048632;
    geom.pDy1 = 1.116997;
    geom.pDx1 = 1.234208;
    geom.pDx2 = 1.212398;
    geom.pDy2 = 1.231781;
    geom.pDx3 = 1.361965;
    geom.pDx4 = 1.337918;
    geom.centralX = 6.439665;
    geom.centralY = 105.782078;
    geom.centralZ = -134.504208;
    geom.pRotationAngleX = -2.482741;
    geom.ModuleSkinThickness = 0.010000;
    geom.NFiberX = 30;
    geom.NFiberY = 48;
    sector_tower_map[geom.id] = geom;
  }
  {
    // tower 548 based Row/Col = 54/8
    geom_tower geom;
    geom.id = 548;
    geom.pDz = 8.326697;
    geom.pTheta = 0.053578;
    geom.pPhi = -0.125683;
    geom.pAlp1 = -0.068143;
    geom.pAlp2 = -0.068127;
    geom.pDy1 = 1.116997;
    geom.pDx1 = 1.235953;
    geom.pDx2 = 1.214052;
    geom.pDy2 = 1.231781;
    geom.pDx3 = 1.363889;
    geom.pDx4 = 1.339743;
    geom.centralX = 9.019696;
    geom.centralY = 105.782078;
    geom.centralZ = -134.504208;
    geom.pRotationAngleX = -2.482741;
    geom.ModuleSkinThickness = 0.010000;
    geom.NFiberX = 30;
    geom.NFiberY = 48;
    sector_tower_map[geom.id] = geom;
  }
}

void PHG4CylinderGeom_Spacalv3::load_demo_sector_tower_map4()
{
  std::cout << "PHG4CylinderGeom_Spacalv3::load_demo_sector_tower_map4 - "
       << "FermiLab test beam 2014. Need to move to calibration database"
       << std::endl;

  // From Oleg's documents

  const int ny = 3;
  const int nx = 3;

  const double inch_to_cm = 2.54;

  const double screen_size_y = 2.108 * inch_to_cm;
  const double screen_size_x = 1.049 * inch_to_cm;

  const double z_screen_6_7 = 0.5 * (6.6 + 6.75) * inch_to_cm;
  const double angle_screen_6_7 = 64.78 / 180. * M_PI;
  const double nawrrow_width_x = screen_size_x * sin(angle_screen_6_7);

  const double z_screen_1_2 = 0.5 * (1.35 + 1.6) * inch_to_cm;
  const double angle_screen_1_2 = 90.56 / 180. * M_PI;
  const double wide_width_x = screen_size_x * sin(angle_screen_1_2);

  const double module_length = z_screen_6_7 - z_screen_1_2;
  assert(module_length > 0);

  // tapering, dxwidth/dlength
  const double tapering_ratio = (wide_width_x - nawrrow_width_x) / module_length;
  assert(tapering_ratio < 1);
  assert(tapering_ratio > 0);

  fiber_clading_thickness = 0.003;
  fiber_core_diameter = 0.050 - fiber_clading_thickness * 2;

  assembly_spacing = 0.002500;

  radius = (nawrrow_width_x) / tapering_ratio;
  thickness = module_length * 1.5;  // keep a large torlerence space
  zmin = -(0.5 * ny * screen_size_y + 2 * assembly_spacing * (ny + 1));
  zmax = -zmin;
  azimuthal_n_sec = floor(2 * M_PI / atan(tapering_ratio));
  max_phi_bin_in_sec = 1;

  const double nawrrow_width_x_construction = radius * 2 * tan(M_PI / azimuthal_n_sec) - 2 * assembly_spacing;
  const double wide_width_x_construction = (radius + module_length) * 2 * tan(M_PI / azimuthal_n_sec) - 2 * assembly_spacing;

  std::cout << "PHG4CylinderGeom_Spacalv3::load_demo_sector_tower_map4 - "

       << "Adjust wide end width by ratio of "
       << wide_width_x_construction / wide_width_x
       << " and narrow end by ratio of "
       << nawrrow_width_x_construction / nawrrow_width_x << std::endl;

  sector_map.clear();
  for (int sec = 0; sec < nx; ++sec)
  {
    const double rot = twopi / (double) (get_azimuthal_n_sec()) * ((double) (sec) -1);

    sector_map[sec] = rot;
  }

  azimuthal_tilt = 0;
  azimuthal_seg_visible = true;
  polar_taper_ratio = 1.;
  sidewall_thickness = 0;
  sidewall_outer_torr = 0;
  sector_tower_map.clear();

  for (int y = 0; y < ny; y++)
  {
    geom_tower geom;
    geom.id = y;
    geom.pDz = module_length / 2;
    geom.pTheta = 0.;
    geom.pPhi = 0.;
    geom.pAlp1 = 0.;
    geom.pAlp2 = 0.;
    geom.pDy1 = 0.5 * screen_size_y;
    geom.pDx1 = 0.5 * nawrrow_width_x_construction;
    geom.pDx2 = 0.5 * nawrrow_width_x_construction;
    geom.pDy2 = 0.5 * screen_size_y;
    geom.pDx3 = 0.5 * wide_width_x_construction;
    geom.pDx4 = 0.5 * wide_width_x_construction;

    geom.centralX = 0.;
    geom.centralY = module_length * 0.5 + radius + assembly_spacing;
    geom.centralZ = screen_size_y * (y - 1);

    geom.pRotationAngleX = -M_PI / 2.;
    geom.ModuleSkinThickness = 0.010000;
    geom.NFiberX = 30;
    geom.NFiberY = 26 * 2 * 2;
    sector_tower_map[geom.id] = geom;
  }
}