g4simulation/g4detectors/PHG4SpacalDetector.cc

0001 /*!
0002  * \file ${file_name}
0003  * \brief
0004  * \author Jin Huang <jhuang@bnl.gov>
0005  * \version $$Revision: 1.7 $$
0006  * \date $$Date: 2015/02/10 15:39:26 $$
0007  */
0008 #include "PHG4SpacalDetector.h"
0009 #include "PHG4CylinderGeom_Spacalv3.h"
0010
0011 #include "PHG4CylinderCellGeom_Spacalv1.h"
0012 #include "PHG4CylinderGeom_Spacalv1.h"  // for PHG4CylinderGeom_Spaca...
0013
0014 #include "PHG4CellDefs.h"
0015 #include "PHG4CylinderCellGeom.h"
0016 #include "PHG4CylinderCellGeomContainer.h"
0017 #include "PHG4CylinderGeomContainer.h"
0018 #include "PHG4SpacalDisplayAction.h"
0019
0020 #include <g4main/PHG4Detector.h>       // for PHG4Detector
0021 #include <g4main/PHG4DisplayAction.h>  // for PHG4DisplayAction
0022 #include <g4main/PHG4Subsystem.h>
0023 #include <g4main/PHG4Utils.h>
0024
0025 #include <phool/PHCompositeNode.h>
0026 #include <phool/PHIODataNode.h>
0027 #include <phool/PHNode.h>          // for PHNode
0028 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0029 #include <phool/PHObject.h>        // for PHObject
0030 #include <phool/getClass.h>
0031 #include <phool/recoConsts.h>
0032
0033 #include <g4gdml/PHG4GDMLConfig.hh>
0034 #include <g4gdml/PHG4GDMLUtility.hh>
0035
0036 #include <calobase/RawTowerDefs.h>           // for convert_name_...
0037 #include <calobase/RawTowerGeom.h>           // for RawTowerGeom
0038 #include <calobase/RawTowerGeomContainer.h>  // for RawTowerGeomC...
0039 #include <calobase/RawTowerGeomContainer_Cylinderv1.h>
0040 #include <calobase/RawTowerGeomv1.h>
0041
0042 #include <TSystem.h>
0043
0044 #include <Geant4/G4LogicalVolume.hh>
0045 #include <Geant4/G4Material.hh>
0046 #include <Geant4/G4PVPlacement.hh>
0047 #include <Geant4/G4PhysicalConstants.hh>
0048 #include <Geant4/G4String.hh>  // for G4String
0049 #include <Geant4/G4SystemOfUnits.hh>
0050 #include <Geant4/G4ThreeVector.hh>  // for G4ThreeVector
0051 #include <Geant4/G4Transform3D.hh>  // for G4Transform3D, G4RotateZ3D
0052 #include <Geant4/G4Tubs.hh>
0053 #include <Geant4/G4Types.hh>  // for G4double
0054 #include <Geant4/G4UserLimits.hh>
0055
0056 #include <cassert>
0057 #include <cstdlib>   // for exit
0058 #include <iostream>  // for operator<<, basic_ostream
0059 #include <sstream>
0060
0061 class PHG4CylinderGeom;
0062
0063 //_______________________________________________________________
0064 // note this inactive thickness is ~1.5% of a radiation length
0065 PHG4SpacalDetector::PHG4SpacalDetector(PHG4Subsystem *subsys,
0066                                        PHCompositeNode *Node,
0067                                        const std::string &dnam,
0068                                        PHParameters *parameters,
0069                                        const int lyr,
0070                                        bool init_geom)
0071   : PHG4Detector(subsys, Node, dnam)
0072   , m_DisplayAction(dynamic_cast<PHG4SpacalDisplayAction *>(subsys->GetDisplayAction()))
0073   , layer(lyr)
0074 {
0075   if (init_geom)
0076   {
0077     _geom = new SpacalGeom_t();
0078     if (_geom == nullptr)
0079     {
0080       std::cout << "PHG4SpacalDetector::Constructor - Fatal Error - invalid geometry object!" << std::endl;
0081       gSystem->Exit(1);
0082       exit(1);
0083     }
0084     assert(parameters);
0085     _geom->ImportParameters(*parameters);
0086
0087     //  _geom->Print();
0088   }
0089
0090   gdml_config = PHG4GDMLUtility::GetOrMakeConfigNode(Node);
0091   assert(gdml_config);
0092 }
0093
0094 PHG4SpacalDetector::~PHG4SpacalDetector()
0095 {
0096   // deleting nullptr pointers is allowed (results in NOOP)
0097   // so checking for not null before deleting is not needed
0098   delete fiber_core_step_limits;
0099   delete _geom;
0100 }
0101
0102 //_______________________________________________________________
0103 int PHG4SpacalDetector::IsInCylinderActive(const G4VPhysicalVolume *volume)
0104 {
0105   //  std::cout << "checking detector" << std::endl;
0106   if (active && fiber_core_vol.find(volume) != fiber_core_vol.end())
0107   {
0108     //      return fiber_core_vol.find(volume)->second;
0109     return FIBER_CORE;
0110   }
0111   if (absorberactive)
0112   {
0113     if (fiber_vol.find(volume) != fiber_vol.end())
0114     {
0115       return FIBER_CLADING;
0116     }
0117
0118     if (block_vol.find(volume) != block_vol.end())
0119     {
0120       return ABSORBER;
0121     }
0122
0123     if (calo_vol.find(volume) != calo_vol.end())
0124     {
0125       return SUPPORT;
0126     }
0127   }
0128   return INACTIVE;
0129 }
0130
0131 //_______________________________________________________________
0132 void PHG4SpacalDetector::ConstructMe(G4LogicalVolume *logicWorld)
0133 {
0134   assert(_geom);
0135   fiber_core_step_limits = new G4UserLimits(_geom->get_fiber_core_step_size() * cm);
0136
0137   Verbosity(_geom->get_construction_verbose());
0138
0139   if ((Verbosity() > 0))
0140   {
0141     std::cout << "PHG4SpacalDetector::Construct::" << GetName()
0142               << " - Start. Print Geometry:" << std::endl;
0143     Print();
0144   }
0145
0146   if ((_geom->get_zmin() * cm + _geom->get_zmax() * cm) / 2 != _geom->get_zpos() * cm)
0147   {
0148     std::cout << "PHG4SpacalDetector::Construct - ERROR - not yet support unsymmetric system. Let me know if you need it. - Jin" << std::endl;
0149     _geom->Print();
0150     gSystem->Exit(-1);
0151   }
0152   if (_geom->get_zmin() * cm >= _geom->get_zmax() * cm)
0153   {
0154     std::cout << "PHG4SpacalDetector::Construct - ERROR - zmin >= zmax!" << std::endl;
0155     _geom->Print();
0156     gSystem->Exit(-1);
0157   }
0158
0159   G4Tubs *cylinder_solid = new G4Tubs(G4String(GetName()),
0160                                       _geom->get_radius() * cm, _geom->get_max_radius() * cm,
0161                                       _geom->get_length() * cm / 2.0, 0, twopi);
0162
0163   recoConsts *rc = recoConsts::instance();
0164   G4Material *cylinder_mat = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
0165   assert(cylinder_mat);
0166
0167   G4LogicalVolume *cylinder_logic = new G4LogicalVolume(cylinder_solid, cylinder_mat, GetName(), nullptr, nullptr, nullptr);
0168   GetDisplayAction()->AddVolume(cylinder_logic, "SpacalCylinder");
0169   if (!m_CosmicSetupFlag)
0170   {
0171     new G4PVPlacement(nullptr, G4ThreeVector(_geom->get_xpos() * cm, _geom->get_ypos() * cm, _geom->get_zpos() * cm),
0172                       cylinder_logic, GetName(),
0173                       logicWorld, false, 0, OverlapCheck());
0174   }
0175   // install sectors
0176   if (_geom->get_sector_map().size() == 0)
0177   {
0178     _geom->init_default_sector_map();
0179   }
0180
0181   if ((Verbosity() > 0))
0182   {
0183     std::cout << "PHG4SpacalDetector::Construct::" << GetName()
0184               << " - start constructing " << _geom->get_sector_map().size() << " sectors in total. " << std::endl;
0185     Print();
0186   }
0187
0188   std::pair<G4LogicalVolume *, G4Transform3D> psec = Construct_AzimuthalSeg();
0189   G4LogicalVolume *sec_logic = psec.first;
0190   const G4Transform3D &sec_trans = psec.second;
0191
0192   for (const SpacalGeom_t::sector_map_t::value_type &val : _geom->get_sector_map())
0193   {
0194     const int sec = val.first;
0195     const double rot = val.second;
0196
0197     G4Transform3D sec_place = G4RotateZ3D(rot) * sec_trans;
0198
0199     std::ostringstream name;
0200     name << GetName() << "_sec" << sec;
0201     G4PVPlacement *calo_phys = nullptr;
0202     if (m_CosmicSetupFlag)
0203     {
0204       calo_phys = new G4PVPlacement(nullptr, G4ThreeVector(0, -(_geom->get_radius()) * cm, 0), sec_logic,
0205                                     G4String(name.str()), logicWorld, false, sec,
0206                                     OverlapCheck());
0207     }
0208     else
0209     {
0210       calo_phys = new G4PVPlacement(sec_place, sec_logic,
0211                                     G4String(name.str()), cylinder_logic, false, sec,
0212                                     OverlapCheck());
0213     }
0214     calo_vol[calo_phys] = sec;
0215
0216     assert(gdml_config);
0217     gdml_config->exclude_physical_vol(calo_phys);
0218   }
0219   _geom->set_nscint(_geom->get_nscint() * _geom->get_sector_map().size());
0220
0221   if (active)
0222   {
0223     std::ostringstream geonode;
0224     if (superdetector != "NONE")
0225     {
0226       geonode << "CYLINDERGEOM_" << superdetector;
0227     }
0228     else
0229     {
0230       geonode << "CYLINDERGEOM_" << detector_type << "_" << layer;
0231     }
0232     PHG4CylinderGeomContainer *geo = findNode::getClass<PHG4CylinderGeomContainer>(topNode(), geonode.str());
0233     if (!geo)
0234     {
0235       geo = new PHG4CylinderGeomContainer();
0236       PHNodeIterator iter(topNode());
0237       PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0238       PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(geo,
0239                                                                    geonode.str(), "PHObject");
0240       runNode->addNode(newNode);
0241     }
0242     // here in the detector class we have internal units, convert to cm
0243     // before putting into the geom object
0244     PHG4CylinderGeom *mygeom = clone_geom();
0245     geo->AddLayerGeom(layer, mygeom);
0246     //    geo->identify();
0247   }
0248
0249   if (absorberactive)
0250   {
0251     std::ostringstream geonode;
0252     if (superdetector != "NONE")
0253     {
0254       geonode << "CYLINDERGEOM_ABSORBER_" << superdetector;
0255     }
0256     else
0257     {
0258       geonode << "CYLINDERGEOM_ABSORBER_" << detector_type << "_" << layer;
0259     }
0260     PHG4CylinderGeomContainer *geo = findNode::getClass<PHG4CylinderGeomContainer>(topNode(), geonode.str());
0261     if (!geo)
0262     {
0263       geo = new PHG4CylinderGeomContainer();
0264       PHNodeIterator iter(topNode());
0265       PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0266       PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(geo,
0267                                                                    geonode.str(), "PHObject");
0268       runNode->addNode(newNode);
0269     }
0270     // here in the detector class we have internal units, convert to cm
0271     // before putting into the geom object
0272     PHG4CylinderGeom *mygeom = clone_geom();
0273     geo->AddLayerGeom(layer, mygeom);
0274     //    geo->identify();
0275   }
0276
0277   if ((Verbosity() > 0))
0278   {
0279     std::cout << "PHG4SpacalDetector::Construct::" << GetName()
0280               << " - Completed. Print Geometry:" << std::endl;
0281     Print();
0282   }
0283 }
0284
0285 std::pair<G4LogicalVolume *, G4Transform3D>
0286 PHG4SpacalDetector::Construct_AzimuthalSeg()
0287 {
0288   G4Tubs *sec_solid = new G4Tubs(G4String(GetName() + std::string("_sec")),
0289                                  _geom->get_radius() * cm, _geom->get_max_radius() * cm,
0290                                  _geom->get_length() * cm / 2.0, 0, twopi / _geom->get_azimuthal_n_sec());
0291
0292   G4Material *cylinder_mat = GetDetectorMaterial(_geom->get_absorber_mat());
0293   assert(cylinder_mat);
0294
0295   G4LogicalVolume *sec_logic = new G4LogicalVolume(sec_solid, cylinder_mat,
0296                                                    G4String(G4String(GetName() + std::string("_sec"))), nullptr, nullptr, nullptr);
0297
0298   GetDisplayAction()->AddVolume(sec_logic, "AzimuthSegment");
0299
0300   const double fiber_length = _geom->get_thickness() * cm - 2 * _geom->get_fiber_outer_r() * cm;
0301   G4LogicalVolume *fiber_logic = Construct_Fiber(fiber_length, std::string(""));
0302
0303   int fiber_count = 0;
0304   //  double z_step = _geom->get_fiber_distance() * cm * sqrt(3) / 2.;
0305   double z_step = _geom->get_z_distance() * cm;
0306   G4double z = _geom->get_zmin() * cm - _geom->get_zpos() * cm + z_step;
0307   while (z < _geom->get_zmax() * cm - _geom->get_zpos() * cm - z_step)
0308   {
0309     const double rot = twopi / _geom->get_azimuthal_n_sec() * ((fiber_count % 2 == 0) ? 1. / 4. : 3. / 4.);
0310
0311     G4Transform3D fiber_place(G4RotateZ3D(rot) * G4TranslateZ3D(z) * G4TranslateX3D(_geom->get_half_radius() * cm) * G4RotateY3D(halfpi));
0312
0313     std::ostringstream name;
0314     name << GetName() << "_fiber_" << fiber_count;
0315
0316     G4PVPlacement *fiber_physi = new G4PVPlacement(fiber_place, fiber_logic,
0317                                                    G4String(name.str()), sec_logic, false, fiber_count,
0318                                                    OverlapCheck());
0319     fiber_vol[fiber_physi] = fiber_count;
0320     assert(gdml_config);
0321     gdml_config->exclude_physical_vol(fiber_physi);
0322
0323     z += z_step;
0324     fiber_count++;
0325   }
0326   _geom->set_nscint(fiber_count);
0327
0328   if (Verbosity() > 0)
0329   {
0330     std::cout << "PHG4SpacalDetector::Construct_AzimuthalSeg::" << GetName()
0331               << " - constructed " << fiber_count << " fibers" << std::endl;
0332     std::cout << "\t"
0333               << "_geom->get_fiber_distance() = " << _geom->get_fiber_distance()
0334               << std::endl;
0335     std::cout << "\t"
0336               << "fiber_length = " << fiber_length
0337               << "_geom->get_azimuthal_distance() = "
0338               << _geom->get_azimuthal_distance() << std::endl;
0339     std::cout << "\t"
0340               << "_geom->is_virualize_fiber() = " << _geom->is_virualize_fiber()
0341               << std::endl;
0342   }
0343
0344   return std::make_pair(sec_logic, G4Transform3D::Identity);
0345 }
0346
0347 G4LogicalVolume *
0348 PHG4SpacalDetector::Construct_Fiber(const G4double length, const std::string &id)
0349 {
0350   G4Tubs *fiber_solid = new G4Tubs(G4String(GetName() + std::string("_fiber") + id),
0351                                    0, _geom->get_fiber_outer_r() * cm, length / 2.0, 0, twopi);
0352
0353   G4Material *clading_mat = GetDetectorMaterial(_geom->get_fiber_clading_mat());
0354   assert(clading_mat);
0355
0356   G4LogicalVolume *fiber_logic = new G4LogicalVolume(fiber_solid, clading_mat,
0357                                                      G4String(G4String(GetName() + std::string("_fiber") + id)), nullptr, nullptr,
0358                                                      nullptr);
0359
0360   {
0361     GetDisplayAction()->AddVolume(fiber_logic, "Fiber");
0362   }
0363
0364   G4Tubs *core_solid = new G4Tubs(
0365       G4String(GetName() + std::string("_fiber_core") + id), 0,
0366       _geom->get_fiber_core_diameter() * cm / 2, length / 2.0, 0, twopi);
0367
0368   G4Material *core_mat = GetDetectorMaterial(_geom->get_fiber_core_mat());
0369   assert(core_mat);
0370
0371   G4LogicalVolume *core_logic = new G4LogicalVolume(core_solid, core_mat,
0372                                                     G4String(G4String(GetName() + std::string("_fiber_core") + id)), nullptr, nullptr,
0373                                                     fiber_core_step_limits);
0374
0375   {
0376     GetDisplayAction()->AddVolume(core_logic, "FiberCore");
0377   }
0378
0379   const bool overlapcheck_fiber = OverlapCheck() and (Verbosity() >= 3);
0380   G4PVPlacement *core_physi = new G4PVPlacement(nullptr, G4ThreeVector(), core_logic,
0381                                                 G4String(G4String(GetName() + std::string("_fiber_core") + id)), fiber_logic,
0382                                                 false, 0, overlapcheck_fiber);
0383   fiber_core_vol[core_physi] = 0;
0384
0385   return fiber_logic;
0386 }
0387
0388 void PHG4SpacalDetector::Print(const std::string & /*what*/) const
0389 {
0390   std::cout << "PHG4SpacalDetector::Print::" << GetName() << " - Print Geometry:" << std::endl;
0391   _geom->Print();
0392
0393   return;
0394 }
0395 // This is dulplicated code, we can get rid of it when we have the code to make towergeom for real data reco.
0396 void PHG4SpacalDetector::AddTowerGeometryNode()
0397 {
0398   PHNodeIterator iter(topNode());
0399   PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0400   if (!runNode)
0401   {
0402     std::cout << PHWHERE << "Run Node missing, doing nothing." << std::endl;
0403     throw std::runtime_error("Failed to find Run node in PHG4SpacalDetector::AddGeometryNode");
0404   }
0405
0406   PHNodeIterator runIter(runNode);
0407   PHCompositeNode *RunDetNode = dynamic_cast<PHCompositeNode *>(runIter.findFirst("PHCompositeNode", superdetector));
0408   if (!RunDetNode)
0409   {
0410     RunDetNode = new PHCompositeNode(superdetector);
0411     runNode->addNode(RunDetNode);
0412   }
0413
0414   const RawTowerDefs::CalorimeterId caloid = RawTowerDefs::convert_name_to_caloid(superdetector);
0415
0416   // get the cell geometry and build up the tower geometry object
0417   std::string geonodename = "CYLINDERCELLGEOM_" + superdetector;
0418   PHG4CylinderCellGeomContainer *cellgeos = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode(), geonodename);
0419   if (!cellgeos)
0420   {
0421     std::cout << PHWHERE << " " << geonodename
0422               << " Node missing, doing nothing." << std::endl;
0423     throw std::runtime_error(
0424         "Failed to find " + geonodename + " node in PHG4SpacalDetector::AddGeometryNode");
0425   }
0426   m_TowerGeomNodeName = "TOWERGEOM_" + superdetector;
0427   m_RawTowerGeomContainer = findNode::getClass<RawTowerGeomContainer>(topNode(), m_TowerGeomNodeName);
0428   if (!m_RawTowerGeomContainer)
0429   {
0430     m_RawTowerGeomContainer = new RawTowerGeomContainer_Cylinderv1(caloid);
0431     PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(m_RawTowerGeomContainer, m_TowerGeomNodeName, "PHObject");
0432     RunDetNode->addNode(newNode);
0433   }
0434   // fill the number of layers in the calorimeter
0435   m_NumLayers = cellgeos->get_NLayers();
0436
0437   // Create the tower nodes on the tree
0438   PHG4CylinderCellGeomContainer::ConstIterator miter;
0439   PHG4CylinderCellGeomContainer::ConstRange begin_end =
0440       cellgeos->get_begin_end();
0441   int ifirst = 1;
0442   int first_layer = -1;
0443   PHG4CylinderCellGeom *first_cellgeo = nullptr;
0444   double inner_radius = 0;
0445   double thickness = 0;
0446   for (miter = begin_end.first; miter != begin_end.second; ++miter)
0447   {
0448     PHG4CylinderCellGeom *cellgeo = miter->second;
0449
0450     if (Verbosity())
0451     {
0452       cellgeo->identify();
0453     }
0454     thickness += cellgeo->get_thickness();
0455     if (ifirst)
0456     {
0457       first_cellgeo = miter->second;
0458       m_CellBinning = cellgeo->get_binning();
0459       m_NumPhiBins = cellgeo->get_phibins();
0460       m_PhiMin = cellgeo->get_phimin();
0461       m_PhiStep = cellgeo->get_phistep();
0462       if (m_CellBinning == PHG4CellDefs::etaphibinning || m_CellBinning == PHG4CellDefs::etaslatbinning)
0463       {
0464         m_NumEtaBins = cellgeo->get_etabins();
0465         m_EtaMin = cellgeo->get_etamin();
0466         m_EtaStep = cellgeo->get_etastep();
0467       }
0468       else if (m_CellBinning == PHG4CellDefs::sizebinning)
0469       {
0470         m_NumEtaBins = cellgeo->get_zbins();  // bin eta in the same number of z bins
0471       }
0472       else if (m_CellBinning == PHG4CellDefs::spacalbinning)
0473       {
0474         // use eta definiton for each row of towers
0475         m_NumEtaBins = cellgeo->get_etabins();
0476       }
0477       else
0478       {
0479         std::cout << "PHG4SpacalDetector::AddGeometryNode::" << GetName()
0480                   << " - Fatal Error - unsupported cell binning method "
0481                   << m_CellBinning << std::endl;
0482       }
0483       inner_radius = cellgeo->get_radius();
0484       first_layer = cellgeo->get_layer();
0485       ifirst = 0;
0486     }
0487     else
0488     {
0489       if (m_CellBinning != cellgeo->get_binning())
0490       {
0491         std::cout << "inconsistent binning method from " << m_CellBinning
0492                   << " layer " << cellgeo->get_layer() << ": "
0493                   << cellgeo->get_binning() << std::endl;
0494         exit(1);
0495       }
0496       if (inner_radius > cellgeo->get_radius())
0497       {
0498         std::cout << "radius of layer " << cellgeo->get_layer() << " is "
0499                   << cellgeo->get_radius() << " which smaller than radius "
0500                   << inner_radius << " of first layer in list " << first_layer
0501                   << std::endl;
0502       }
0503       if (m_NumPhiBins != cellgeo->get_phibins())
0504       {
0505         std::cout << "mixing number of phibins, fisrt layer: " << m_NumPhiBins
0506                   << " layer " << cellgeo->get_layer() << ": "
0507                   << cellgeo->get_phibins() << std::endl;
0508         exit(1);
0509       }
0510       if (m_PhiMin != cellgeo->get_phimin())
0511       {
0512         std::cout << "mixing number of phimin, fisrt layer: " << m_PhiMin
0513                   << " layer " << cellgeo->get_layer() << ": "
0514                   << cellgeo->get_phimin() << std::endl;
0515         exit(1);
0516       }
0517       if (m_PhiStep != cellgeo->get_phistep())
0518       {
0519         std::cout << "mixing phisteps first layer: " << m_PhiStep << " layer "
0520                   << cellgeo->get_layer() << ": " << cellgeo->get_phistep()
0521                   << " diff: " << m_PhiStep - cellgeo->get_phistep() << std::endl;
0522         exit(1);
0523       }
0524       if (m_CellBinning == PHG4CellDefs::etaphibinning || m_CellBinning == PHG4CellDefs::etaslatbinning)
0525       {
0526         if (m_NumEtaBins != cellgeo->get_etabins())
0527         {
0528           std::cout << "mixing number of EtaBins , first layer: "
0529                     << m_NumEtaBins << " layer " << cellgeo->get_layer() << ": "
0530                     << cellgeo->get_etabins() << std::endl;
0531           exit(1);
0532         }
0533         if (fabs(m_EtaMin - cellgeo->get_etamin()) > 1e-9)
0534         {
0535           std::cout << "mixing etamin, fisrt layer: " << m_EtaMin << " layer "
0536                     << cellgeo->get_layer() << ": " << cellgeo->get_etamin()
0537                     << " diff: " << m_EtaMin - cellgeo->get_etamin() << std::endl;
0538           exit(1);
0539         }
0540         if (fabs(m_EtaStep - cellgeo->get_etastep()) > 1e-9)
0541         {
0542           std::cout << "mixing eta steps first layer: " << m_EtaStep
0543                     << " layer " << cellgeo->get_layer() << ": "
0544                     << cellgeo->get_etastep() << " diff: "
0545                     << m_EtaStep - cellgeo->get_etastep() << std::endl;
0546           exit(1);
0547         }
0548       }
0549
0550       else if (m_CellBinning == PHG4CellDefs::sizebinning)
0551       {
0552         if (m_NumEtaBins != cellgeo->get_zbins())
0553         {
0554           std::cout << "mixing number of z bins , first layer: " << m_NumEtaBins
0555                     << " layer " << cellgeo->get_layer() << ": "
0556                     << cellgeo->get_zbins() << std::endl;
0557           exit(1);
0558         }
0559       }
0560     }
0561   }
0562   m_RawTowerGeomContainer->set_radius(inner_radius);
0563   m_RawTowerGeomContainer->set_thickness(thickness);
0564   m_RawTowerGeomContainer->set_phibins(m_NumPhiBins);
0565   //  m_RawTowerGeomContainer->set_phistep(m_PhiStep);
0566   //  m_RawTowerGeomContainer->set_phimin(m_PhiMin);
0567   m_RawTowerGeomContainer->set_etabins(m_NumEtaBins);
0568
0569   if (!first_cellgeo)
0570   {
0571     std::cout << "PHG4SpacalDetector::AddGeometryNode - ERROR - can not find first layer of cells "
0572               << std::endl;
0573
0574     exit(1);
0575   }
0576
0577   for (int ibin = 0; ibin < first_cellgeo->get_phibins(); ibin++)
0578   {
0579     const std::pair<double, double> range = first_cellgeo->get_phibounds(ibin);
0580
0581     m_RawTowerGeomContainer->set_phibounds(ibin, range);
0582   }
0583
0584   if (m_CellBinning == PHG4CellDefs::etaphibinning || m_CellBinning == PHG4CellDefs::etaslatbinning || m_CellBinning == PHG4CellDefs::spacalbinning)
0585   {
0586     const double r = inner_radius;
0587
0588     for (int ibin = 0; ibin < first_cellgeo->get_etabins(); ibin++)
0589     {
0590       const std::pair<double, double> range = first_cellgeo->get_etabounds(ibin);
0591
0592       m_RawTowerGeomContainer->set_etabounds(ibin, range);
0593     }
0594
0595     // setup location of all towers
0596     for (int iphi = 0; iphi < m_RawTowerGeomContainer->get_phibins(); iphi++)
0597     {
0598       for (int ieta = 0; ieta < m_RawTowerGeomContainer->get_etabins(); ieta++)
0599       {
0600         const RawTowerDefs::keytype key =
0601             RawTowerDefs::encode_towerid(caloid, ieta, iphi);
0602
0603         const double x(r * cos(m_RawTowerGeomContainer->get_phicenter(iphi)));
0604         const double y(r * sin(m_RawTowerGeomContainer->get_phicenter(iphi)));
0605         const double z(r / tan(PHG4Utils::get_theta(m_RawTowerGeomContainer->get_etacenter(ieta))));
0606
0607         RawTowerGeom *tg = m_RawTowerGeomContainer->get_tower_geometry(key);
0608         if (tg)
0609         {
0610           if (Verbosity() > 0)
0611           {
0612             std::cout << "PHG4SpacalDetector::AddGeometryNode - Tower geometry " << key << " already exists" << std::endl;
0613           }
0614
0615           if (fabs(tg->get_center_x() - x) > 1e-4)
0616           {
0617             std::cout << "PHG4SpacalDetector::AddGeometryNode - Fatal Error - duplicated Tower geometry " << key << " with existing x = " << tg->get_center_x() << " and expected x = " << x
0618                       << std::endl;
0619
0620             exit(1);
0621           }
0622           if (fabs(tg->get_center_y() - y) > 1e-4)
0623           {
0624             std::cout << "PHG4SpacalDetector::AddGeometryNode - Fatal Error - duplicated Tower geometry " << key << " with existing y = " << tg->get_center_y() << " and expected y = " << y
0625                       << std::endl;
0626             exit(1);
0627           }
0628           if (fabs(tg->get_center_z() - z) > 1e-4)
0629           {
0630             std::cout << "PHG4SpacalDetector::AddGeometryNode - Fatal Error - duplicated Tower geometry " << key << " with existing z= " << tg->get_center_z() << " and expected z = " << z
0631                       << std::endl;
0632             exit(1);
0633           }
0634         }
0635         else
0636         {
0637           if (Verbosity() > 0)
0638           {
0639             std::cout << "PHG4SpacalDetector::AddGeometryNode - building tower geometry " << key << "" << std::endl;
0640           }
0641
0642           tg = new RawTowerGeomv1(key);
0643
0644           tg->set_center_x(x);
0645           tg->set_center_y(y);
0646           tg->set_center_z(z);
0647           m_RawTowerGeomContainer->add_tower_geometry(tg);
0648         }
0649       }
0650     }
0651   }
0652   else if (m_CellBinning == PHG4CellDefs::sizebinning)
0653   {
0654     const double r = inner_radius;
0655
0656     for (int ibin = 0; ibin < first_cellgeo->get_zbins(); ibin++)
0657     {
0658       const std::pair<double, double> z_range = first_cellgeo->get_zbounds(ibin);
0659       //          const double r = first_cellgeo->get_radius();
0660       const double eta1 = -log(tan(atan2(r, z_range.first) / 2.));
0661       const double eta2 = -log(tan(atan2(r, z_range.second) / 2.));
0662       m_RawTowerGeomContainer->set_etabounds(ibin, std::make_pair(eta1, eta2));
0663     }
0664
0665     // setup location of all towers
0666     for (int iphi = 0; iphi < m_RawTowerGeomContainer->get_phibins(); iphi++)
0667     {
0668       for (int ibin = 0; ibin < first_cellgeo->get_zbins(); ibin++)
0669       {
0670         const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(caloid, ibin, iphi);
0671
0672         const double x(r * cos(m_RawTowerGeomContainer->get_phicenter(iphi)));
0673         const double y(r * sin(m_RawTowerGeomContainer->get_phicenter(iphi)));
0674         const double z(first_cellgeo->get_zcenter(ibin));
0675
0676         RawTowerGeom *tg = m_RawTowerGeomContainer->get_tower_geometry(key);
0677         if (tg)
0678         {
0679           if (Verbosity() > 0)
0680           {
0681             std::cout << "PHG4SpacalDetector::AddGeometryNode - Tower geometry " << key << " already exists" << std::endl;
0682           }
0683
0684           if (fabs(tg->get_center_x() - x) > 1e-4)
0685           {
0686             std::cout << "PHG4SpacalDetector::AddGeometryNode - Fatal Error - duplicated Tower geometry " << key << " with existing x = " << tg->get_center_x() << " and expected x = " << x
0687                       << std::endl;
0688
0689             exit(1);
0690           }
0691           if (fabs(tg->get_center_y() - y) > 1e-4)
0692           {
0693             std::cout << "PHG4SpacalDetector::AddGeometryNode - Fatal Error - duplicated Tower geometry " << key << " with existing y = " << tg->get_center_y() << " and expected y = " << y
0694                       << std::endl;
0695             exit(1);
0696           }
0697           if (fabs(tg->get_center_z() - z) > 1e-4)
0698           {
0699             std::cout << "PHG4SpacalDetector::AddGeometryNode - Fatal Error - duplicated Tower geometry " << key << " with existing z= " << tg->get_center_z() << " and expected z = " << z
0700                       << std::endl;
0701             exit(1);
0702           }
0703         }
0704         else
0705         {
0706           if (Verbosity() > 0)
0707           {
0708             std::cout << "PHG4SpacalDetector::AddGeometryNode - building tower geometry " << key << "" << std::endl;
0709           }
0710
0711           tg = new RawTowerGeomv1(key);
0712
0713           tg->set_center_x(x);
0714           tg->set_center_y(y);
0715           tg->set_center_z(z);
0716           m_RawTowerGeomContainer->add_tower_geometry(tg);
0717         }
0718       }
0719     }
0720   }
0721   else
0722   {
0723     std::cout << "PHG4SpacalDetector::AddGeometryNode - ERROR - unsupported cell geometry "
0724               << m_CellBinning << std::endl;
0725     exit(1);
0726   }
0727
0728   if (Verbosity() >= 1)
0729   {
0730     m_RawTowerGeomContainer->identify();
0731   }
0732
0733   return;
0734 }
0735
0736 void PHG4SpacalDetector::AddCellGeometryNode()
0737 {
0738   PHNodeIterator iter(topNode());
0739   std::string detector = superdetector;
0740   // Looking for the DST node
0741   PHCompositeNode *dstNode;
0742   dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0743   if (!dstNode)
0744   {
0745     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0746     exit(1);
0747   }
0748
0749   std::string geonodename = "CYLINDERGEOM_" + detector;
0750   PHG4CylinderGeomContainer *geo = findNode::getClass<PHG4CylinderGeomContainer>(topNode(), geonodename);
0751   if (!geo)
0752   {
0753     std::cout << "PHG4SpacalDetector::AddCellGeometryNode - Could not locate geometry node "
0754               << geonodename << std::endl;
0755     topNode()->print();
0756     exit(1);
0757   }
0758   if (Verbosity() > 0)
0759   {
0760     std::cout << "PHG4SpacalDetector::AddCellGeometryNode- incoming geometry:"
0761               << std::endl;
0762     geo->identify();
0763   }
0764   std::string seggeonodename = "CYLINDERCELLGEOM_" + detector;
0765   PHG4CylinderCellGeomContainer *seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode(), seggeonodename);
0766   if (!seggeo)
0767   {
0768     seggeo = new PHG4CylinderCellGeomContainer();
0769     PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0770     PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(seggeo, seggeonodename, "PHObject");
0771     runNode->addNode(newNode);
0772   }
0773
0774   const PHG4CylinderGeom *layergeom_raw = geo->GetFirstLayerGeom();
0775   assert(layergeom_raw);
0776
0777   // a special implimentation of PHG4CylinderGeom is required here.
0778   const PHG4CylinderGeom_Spacalv3 *layergeom =
0779       dynamic_cast<const PHG4CylinderGeom_Spacalv3 *>(layergeom_raw);
0780
0781   if (!layergeom)
0782   {
0783     std::cout << "PHG4SpacalDetector::AddCellGeometryNode- Fatal Error -"
0784               << " require to work with a version of SPACAL geometry of PHG4CylinderGeom_Spacalv3 or higher. "
0785               << "However the incoming geometry has version "
0786               << layergeom_raw->ClassName() << std::endl;
0787     exit(1);
0788   }
0789   if (Verbosity() > 1)
0790   {
0791     layergeom->identify();
0792   }
0793
0794   layergeom->subtower_consistency_check();
0795
0796   //      int layer = layergeom->get_layer();
0797
0798   // create geo object and fill with variables common to all binning methods
0799   PHG4CylinderCellGeom_Spacalv1 *layerseggeo = new PHG4CylinderCellGeom_Spacalv1();
0800   layerseggeo->set_layer(layergeom->get_layer());
0801   layerseggeo->set_radius(layergeom->get_radius());
0802   layerseggeo->set_thickness(layergeom->get_thickness());
0803   layerseggeo->set_binning(PHG4CellDefs::spacalbinning);
0804
0805   // construct a map to convert tower_ID into the older eta bins.
0806
0807   const PHG4CylinderGeom_Spacalv3::tower_map_t &tower_map = layergeom->get_sector_tower_map();
0808   const PHG4CylinderGeom_Spacalv3::sector_map_t &sector_map = layergeom->get_sector_map();
0809   const int nphibin = layergeom->get_azimuthal_n_sec()       // sector
0810                       * layergeom->get_max_phi_bin_in_sec()  // blocks per sector
0811                       * layergeom->get_n_subtower_phi();     // subtower per block
0812   const double deltaphi = 2. * M_PI / nphibin;
0813
0814   using map_z_tower_z_ID_t = std::map<double, int>;
0815   map_z_tower_z_ID_t map_z_tower_z_ID;
0816   double phi_min = NAN;
0817
0818   for (const auto &tower_pair : tower_map)
0819   {
0820     const int &tower_ID = tower_pair.first;
0821     const PHG4CylinderGeom_Spacalv3::geom_tower &tower = tower_pair.second;
0822
0823     // inspect index in sector 0
0824     std::pair<int, int> tower_z_phi_ID = layergeom->get_tower_z_phi_ID(tower_ID, 0);
0825
0826     const int &tower_ID_z = tower_z_phi_ID.first;
0827     const int &tower_ID_phi = tower_z_phi_ID.second;
0828
0829     if (tower_ID_phi == 0)
0830     {
0831       // assign phi min according phi bin 0
0832 // NOLINTNEXTLINE(bugprone-integer-division)
0833       phi_min = M_PI_2 - deltaphi * (layergeom->get_max_phi_bin_in_sec() * layergeom->get_n_subtower_phi() / 2)  // shift of first tower in sector
0834                 + sector_map.begin()->second;
0835     }
0836
0837     if (tower_ID_phi == layergeom->get_max_phi_bin_in_sec() / 2)
0838     {
0839       // assign eta min according phi bin 0
0840       map_z_tower_z_ID[tower.centralZ] = tower_ID_z;
0841     }
0842     // ...
0843   }  //       const auto &tower_pair: tower_map
0844
0845   assert(!std::isnan(phi_min));
0846   layerseggeo->set_phimin(phi_min);
0847   layerseggeo->set_phistep(deltaphi);
0848   layerseggeo->set_phibins(nphibin);
0849
0850   PHG4CylinderCellGeom_Spacalv1::tower_z_ID_eta_bin_map_t tower_z_ID_eta_bin_map;
0851   int eta_bin = 0;
0852   for (auto &z_tower_z_ID : map_z_tower_z_ID)
0853   {
0854     tower_z_ID_eta_bin_map[z_tower_z_ID.second] = eta_bin;
0855     eta_bin++;
0856   }
0857   layerseggeo->set_tower_z_ID_eta_bin_map(tower_z_ID_eta_bin_map);
0858   layerseggeo->set_etabins(eta_bin * layergeom->get_n_subtower_eta());
0859   layerseggeo->set_etamin(NAN);
0860   layerseggeo->set_etastep(NAN);
0861
0862   // build eta bin maps
0863   for (const auto &tower_pair : tower_map)
0864   {
0865     const int &tower_ID = tower_pair.first;
0866     const PHG4CylinderGeom_Spacalv3::geom_tower &tower = tower_pair.second;
0867
0868     // inspect index in sector 0
0869     std::pair<int, int> tower_z_phi_ID = layergeom->get_tower_z_phi_ID(tower_ID, 0);
0870     const int &tower_ID_z = tower_z_phi_ID.first;
0871     const int &tower_ID_phi = tower_z_phi_ID.second;
0872     const int &etabin = tower_z_ID_eta_bin_map[tower_ID_z];
0873
0874     if (tower_ID_phi == layergeom->get_max_phi_bin_in_sec() / 2)
0875     {
0876       // half z-range
0877       const double dz = fabs(0.5 * (tower.pDy1 + tower.pDy2) / sin(tower.pRotationAngleX));
0878       const double tower_radial = layergeom->get_tower_radial_position(tower);
0879
0880       auto z_to_eta = [&tower_radial](const double &z)
0881       { return -log(tan(0.5 * atan2(tower_radial, z))); };
0882
0883       const double eta_central = z_to_eta(tower.centralZ);
0884       // half eta-range
0885       const double deta = (z_to_eta(tower.centralZ + dz) - z_to_eta(tower.centralZ - dz)) / 2;
0886       assert(deta > 0);
0887
0888       for (int sub_tower_ID_y = 0; sub_tower_ID_y < tower.NSubtowerY;
0889            ++sub_tower_ID_y)
0890       {
0891         assert(tower.NSubtowerY <= layergeom->get_n_subtower_eta());
0892         // do not overlap to the next bin.
0893         const int sub_tower_etabin = etabin * layergeom->get_n_subtower_eta() + sub_tower_ID_y;
0894
0895         const std::pair<double, double> etabounds(eta_central - deta + sub_tower_ID_y * 2 * deta / tower.NSubtowerY,
0896                                                   eta_central - deta + (sub_tower_ID_y + 1) * 2 * deta / tower.NSubtowerY);
0897
0898         const std::pair<double, double> zbounds(tower.centralZ - dz + sub_tower_ID_y * 2 * dz / tower.NSubtowerY,
0899                                                 tower.centralZ - dz + (sub_tower_ID_y + 1) * 2 * dz / tower.NSubtowerY);
0900
0901         layerseggeo->set_etabounds(sub_tower_etabin, etabounds);
0902         layerseggeo->set_zbounds(sub_tower_etabin, zbounds);
0903       }
0904     }
0905     // ...
0906   }  //       const auto &tower_pair: tower_map
0907
0908   // add geo object filled by different binning methods
0909   seggeo->AddLayerCellGeom(layerseggeo);
0910   // save this to the run wise tree to store on DST
0911   PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0912   // PHCompositeNode *parNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "PAR"));
0913   PHNodeIterator runIter(runNode);
0914   PHCompositeNode *RunDetNode = dynamic_cast<PHCompositeNode *>(runIter.findFirst("PHCompositeNode", detector));
0915   if (!RunDetNode)
0916   {
0917     RunDetNode = new PHCompositeNode(detector);
0918     runNode->addNode(RunDetNode);
0919   }
0920   /*
0921   std::string paramnodename = "G4CELLPARAM_" + detector;
0922   SaveToNodeTree(RunDetNode, paramnodename);
0923   save this to the parNode for use
0924   PHNodeIterator parIter(parNode);
0925    PHCompositeNode *ParDetNode = dynamic_cast<PHCompositeNode *>(parIter.findFirst("PHCompositeNode", detector));
0926    if (!ParDetNode)
0927    {
0928      ParDetNode = new PHCompositeNode(detector);
0929      parNode->addNode(ParDetNode);
0930    }
0931    std::string cellgeonodename = "G4CELLGEO_" + detector;
0932    PutOnParNode(ParDetNode, cellgeonodename);
0933     */
0934   return;
0935 }