g4simulation/g4ihcal/PHG4IHCalDetector.cc

0001 #include "PHG4IHCalDetector.h"
0002
0003 #include "PHG4IHCalDisplayAction.h"
0004
0005 #include <g4detectors/PHG4DetectorSubsystem.h>
0006 #include <g4detectors/PHG4HcalDefs.h>
0007
0008 #include <phparameter/PHParameters.h>
0009
0010 #include <g4gdml/PHG4GDMLConfig.hh>
0011 #include <g4gdml/PHG4GDMLUtility.hh>
0012
0013 #include <calobase/RawTowerDefs.h>           // for convert_name_...
0014 #include <calobase/RawTowerGeom.h>           // for RawTowerGeom
0015 #include <calobase/RawTowerGeomContainer.h>  // for RawTowerGeomC...
0016 #include <calobase/RawTowerGeomContainer_Cylinderv1.h>
0017 #include <calobase/RawTowerGeomv1.h>
0018
0019 #include <g4main/PHG4Detector.h>
0020 #include <g4main/PHG4DisplayAction.h>
0021 #include <g4main/PHG4Subsystem.h>
0022 #include <g4main/PHG4Utils.h>
0023
0024 #include <ffamodules/CDBInterface.h>
0025
0026 #include <phool/PHCompositeNode.h>
0027 #include <phool/PHIODataNode.h>
0028 #include <phool/PHNode.h>  // for PHNode
0029 #include <phool/PHNodeIterator.h>
0030 #include <phool/PHObject.h>  // for PHObject
0031 #include <phool/getClass.h>
0032 #include <phool/phool.h>
0033 #include <phool/recoConsts.h>
0034
0035 #include <TSystem.h>
0036
0037 #include <Geant4/G4AssemblyVolume.hh>
0038 #include <Geant4/G4IonisParamMat.hh>
0039 #include <Geant4/G4LogicalVolume.hh>
0040 #include <Geant4/G4Material.hh>
0041 #include <Geant4/G4MaterialTable.hh>
0042 #include <Geant4/G4PVPlacement.hh>
0043 #include <Geant4/G4RotationMatrix.hh>
0044 #include <Geant4/G4String.hh>
0045 #include <Geant4/G4SystemOfUnits.hh>
0046 #include <Geant4/G4ThreeVector.hh>
0047 #include <Geant4/G4Transform3D.hh>
0048 #include <Geant4/G4Tubs.hh>
0049 #include <Geant4/G4VPhysicalVolume.hh>  // for G4VPhysicalVolume
0050 #include <Geant4/G4VSolid.hh>
0051
0052 #pragma GCC diagnostic push
0053 #pragma GCC diagnostic ignored "-Wshadow"
0054 #pragma GCC diagnostic ignored "-Wpedantic"
0055 #include <Geant4/G4GDMLParser.hh>
0056 #include <Geant4/G4GDMLReadStructure.hh>  // for G4GDMLReadStructure
0057 #pragma GCC diagnostic pop
0058
0059 #include <boost/lexical_cast.hpp>
0060 #include <boost/tokenizer.hpp>
0061
0062 #include <cassert>
0063 #include <cmath>
0064 #include <cstdlib>
0065 #include <filesystem>
0066 #include <iostream>
0067 #include <memory>   // for unique_ptr
0068 #include <utility>  // for pair, make_pair
0069 #include <vector>   // for vector, vector<>::iter...
0070
0071 PHG4IHCalDetector::PHG4IHCalDetector(PHG4Subsystem *subsys, PHCompositeNode *Node, PHParameters *parameters, const std::string &dnam)
0072   : PHG4Detector(subsys, Node, dnam)
0073   , m_DisplayAction(dynamic_cast<PHG4IHCalDisplayAction *>(subsys->GetDisplayAction()))
0074   , m_Params(parameters)
0075   , gdml_config(PHG4GDMLUtility::GetOrMakeConfigNode(Node))
0076   , m_InnerRadius(m_Params->get_double_param("inner_radius") * cm)
0077   , m_OuterRadius(m_Params->get_double_param("outer_radius") * cm)
0078   , m_SizeZ(m_Params->get_double_param("size_z") * cm)
0079   , m_NumScintiPlates(m_Params->get_int_param(PHG4HcalDefs::scipertwr) * m_Params->get_int_param("n_towers"))
0080   , m_Active(m_Params->get_int_param("active"))
0081   , m_AbsorberActive(m_Params->get_int_param("absorberactive"))
0082   , m_GDMPath(m_Params->get_string_param("GDMPath"))
0083 {
0084   assert(gdml_config);
0085   // changes in the parameters have to be made here
0086   // otherwise they will not be propagated to the node tree
0087   if (std::filesystem::path(m_GDMPath).extension() != ".gdml")
0088   {
0089     m_GDMPath = CDBInterface::instance()->getUrl(m_GDMPath);
0090     m_Params->set_string_param("GDMPath", m_GDMPath);
0091   }
0092 }
0093
0094 PHG4IHCalDetector::~PHG4IHCalDetector()
0095 {
0096   delete m_ScintiMotherAssembly;
0097 }
0098
0099 //_______________________________________________________________
0100 //_______________________________________________________________
0101 int PHG4IHCalDetector::IsInIHCal(G4VPhysicalVolume *volume) const
0102 {
0103   if (m_AbsorberActive)
0104   {
0105     if (m_SteelAbsorberLogVolSet.contains(volume->GetLogicalVolume()))
0106     {
0107       return -1;
0108     }
0109   }
0110   if (m_Active)
0111   {
0112     if (m_ScintiTileLogVolSet.contains(volume->GetLogicalVolume()))
0113     {
0114       return 1;
0115     }
0116   }
0117   return 0;
0118 }
0119
0120 // Construct the envelope and the call the
0121 // actual inner hcal construction
0122 void PHG4IHCalDetector::ConstructMe(G4LogicalVolume *logicWorld)
0123 {
0124   recoConsts *rc = recoConsts::instance();
0125   G4Material *worldmat = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
0126   G4VSolid *hcal_envelope_cylinder = new G4Tubs("IHCal_envelope_solid", m_InnerRadius, m_OuterRadius, m_SizeZ / 2., 0, 2 * M_PI);
0127   m_VolumeEnvelope = hcal_envelope_cylinder->GetCubicVolume();
0128   G4LogicalVolume *hcal_envelope_log = new G4LogicalVolume(hcal_envelope_cylinder, worldmat, "Hcal_envelope", nullptr, nullptr, nullptr);
0129
0130   G4RotationMatrix hcal_rotm;
0131   hcal_rotm.rotateX(m_Params->get_double_param("rot_x") * deg);
0132   hcal_rotm.rotateY(m_Params->get_double_param("rot_y") * deg);
0133   hcal_rotm.rotateZ(m_Params->get_double_param("rot_z") * deg);
0134   G4VPhysicalVolume *mothervol = new G4PVPlacement(G4Transform3D(hcal_rotm, G4ThreeVector(m_Params->get_double_param("place_x") * cm, m_Params->get_double_param("place_y") * cm, m_Params->get_double_param("place_z") * cm)), hcal_envelope_log, "IHCalEnvelope", logicWorld, false, false, OverlapCheck());
0135   m_DisplayAction->SetMyTopVolume(mothervol);
0136   ConstructIHCal(hcal_envelope_log);
0137
0138   // disable GDML export for HCal geometries for memory saving and compatibility issues
0139   assert(gdml_config);
0140   gdml_config->exclude_physical_vol(mothervol);
0141   gdml_config->exclude_logical_vol(hcal_envelope_log);
0142
0143   const G4MaterialTable *mtable = G4Material::GetMaterialTable();
0144   int nMaterials = G4Material::GetNumberOfMaterials();
0145   for (auto i = 0; i < nMaterials; ++i)
0146   {
0147     const G4Material *mat = (*mtable)[i];
0148     if (mat->GetName() == "Uniplast_scintillator")
0149     {
0150       if ((mat->GetIonisation()->GetBirksConstant()) == 0)
0151       {
0152         mat->GetIonisation()->SetBirksConstant(m_Params->get_double_param("Birk_const"));
0153       }
0154     }
0155   }
0156   if (!m_Params->get_int_param("saveg4hit"))
0157   {
0158     AddGeometryNode();
0159   }
0160   return;
0161 }
0162
0163 int PHG4IHCalDetector::ConstructIHCal(G4LogicalVolume *hcalenvelope)
0164 {
0165   // import the staves from the geometry file
0166   std::unique_ptr<G4GDMLReadStructure> reader(new G4GDMLReadStructure());
0167   G4GDMLParser gdmlParser(reader.get());
0168   gdmlParser.SetOverlapCheck(OverlapCheck());
0169   if (!std::filesystem::exists(m_GDMPath))
0170   {
0171     std::cout << PHWHERE << " Inner HCal gdml file " << m_GDMPath << " not found" << std::endl;
0172     gSystem->Exit(1);
0173     exit(1);
0174   }
0175   gdmlParser.Read(m_GDMPath, false);
0176
0177   G4AssemblyVolume *abs_asym = reader->GetAssembly("InnerSector");      // absorber
0178   m_ScintiMotherAssembly = reader->GetAssembly("InnerTileAssembly90");  // scintillator
0179   std::vector<G4VPhysicalVolume *>::iterator it = abs_asym->GetVolumesIterator();
0180   static const unsigned int tilepersec = 24 * 4 * 2;
0181   for (unsigned int isector = 0; isector < abs_asym->TotalImprintedVolumes(); isector++)
0182   {
0183     m_DisplayAction->AddSteelVolume((*it)->GetLogicalVolume());
0184     m_SteelAbsorberLogVolSet.insert((*it)->GetLogicalVolume());
0185     hcalenvelope->AddDaughter((*it));
0186     m_AbsorberPhysVolMap.insert(std::make_pair(*it, isector));
0187     m_VolumeSteel += (*it)->GetLogicalVolume()->GetSolid()->GetCubicVolume();
0188     std::vector<G4VPhysicalVolume *>::iterator its = m_ScintiMotherAssembly->GetVolumesIterator();
0189     unsigned int ioff = isector * tilepersec;
0190     for (unsigned int j = 0; j < ioff; j++)
0191     {
0192       ++its;
0193     }
0194     for (unsigned int j = ioff; j < ioff + tilepersec; j++)
0195     {
0196       m_DisplayAction->AddScintiVolume((*its)->GetLogicalVolume());
0197       m_ScintiTileLogVolSet.insert((*its)->GetLogicalVolume());
0198       hcalenvelope->AddDaughter((*its));
0199       m_ScintiTilePhysVolMap.insert(std::make_pair(*its, ExtractLayerTowerId(isector, *its)));
0200       m_VolumeScintillator += (*its)->GetLogicalVolume()->GetSolid()->GetCubicVolume();
0201       ++its;
0202     }
0203
0204     ++it;
0205   }
0206   return 0;
0207 }
0208
0209 int PHG4IHCalDetector::ConsistencyCheck() const
0210 {
0211   // just make sure the parameters make a bit of sense
0212   if (m_InnerRadius >= m_OuterRadius)
0213   {
0214     std::cout << PHWHERE << ": Inner Radius " << m_InnerRadius / cm
0215               << " cm larger than Outer Radius " << m_OuterRadius / cm
0216               << " cm" << std::endl;
0217     gSystem->Exit(1);
0218   }
0219   return 0;
0220 }
0221
0222 void PHG4IHCalDetector::Print(const std::string &what) const
0223 {
0224   std::cout << "Inner Hcal Detector:" << std::endl;
0225   if (what == "ALL" || what == "VOLUME")
0226   {
0227     std::cout << "Volume Envelope: " << m_VolumeEnvelope / cm3 << " cm^3" << std::endl;
0228     std::cout << "Volume Steel: " << m_VolumeSteel / cm3 << " cm^3" << std::endl;
0229     std::cout << "Volume Scintillator: " << m_VolumeScintillator / cm3 << " cm^3" << std::endl;
0230     std::cout << "Volume Air: " << (m_VolumeEnvelope - m_VolumeSteel - m_VolumeScintillator) / cm3 << " cm^3" << std::endl;
0231   }
0232   return;
0233 }
0234
0235 std::tuple<int, int, int> PHG4IHCalDetector::GetLayerTowerId(G4VPhysicalVolume *volume) const
0236 {
0237   auto it = m_ScintiTilePhysVolMap.find(volume);
0238   if (it != m_ScintiTilePhysVolMap.end())
0239   {
0240     return it->second;
0241   }
0242   std::cout << "could not locate volume " << volume->GetName()
0243             << " in Inner Hcal scintillator map" << std::endl;
0244   gSystem->Exit(1);
0245   // that's dumb but code checkers do not know that gSystem->Exit()
0246   // terminates, so using the standard exit() makes them happy
0247   exit(1);
0248 }
0249
0250 int PHG4IHCalDetector::GetSectorId(G4VPhysicalVolume *volume) const
0251 {
0252   auto it = m_AbsorberPhysVolMap.find(volume);
0253   if (it != m_AbsorberPhysVolMap.end())
0254   {
0255     return it->second;
0256   }
0257   std::cout << "could not locate volume " << volume->GetName()
0258             << " in Inner Hcal Absorber map" << std::endl;
0259   gSystem->Exit(1);
0260   // that's dumb but code checkers do not know that gSystem->Exit()
0261   // terminates, so using the standard exit() makes them happy
0262   exit(1);
0263 }
0264
0265 std::tuple<int, int, int> PHG4IHCalDetector::ExtractLayerTowerId(const unsigned int isector, G4VPhysicalVolume *volume)
0266 {
0267   boost::char_separator<char> sep("_");
0268   boost::tokenizer<boost::char_separator<char>> tok(volume->GetName(), sep);
0269   boost::tokenizer<boost::char_separator<char>>::const_iterator tokeniter;
0270   int layer_id = -1;
0271   int tower_id = -1;
0272   for (tokeniter = tok.begin(); tokeniter != tok.end(); ++tokeniter)
0273   {
0274     if (*tokeniter == "impr")
0275     {
0276       ++tokeniter;
0277       if (tokeniter != tok.end())
0278       {
0279         layer_id = boost::lexical_cast<int>(*tokeniter) / 2;
0280         layer_id--;
0281         if (layer_id < 0 || layer_id >= m_NumScintiPlates)
0282         {
0283           std::cout << "invalid scintillator row " << layer_id
0284                     << ", valid range 0 < row < " << m_NumScintiPlates << std::endl;
0285           gSystem->Exit(1);
0286         }
0287       }
0288       else
0289       {
0290         std::cout << PHWHERE << " Error parsing " << volume->GetName()
0291                   << " for mother volume number " << std::endl;
0292         gSystem->Exit(1);
0293       }
0294       break;
0295     }
0296   }
0297   for (tokeniter = tok.begin(); tokeniter != tok.end(); ++tokeniter)
0298   {
0299     if (*tokeniter == "pv")
0300     {
0301       ++tokeniter;
0302       if (tokeniter != tok.end())
0303       {
0304         tower_id = boost::lexical_cast<int>(*tokeniter);
0305       }
0306     }
0307   }
0308   int column = map_towerid(tower_id);
0309   int row = map_layerid(layer_id);
0310   return std::make_tuple(isector, row, column);
0311 }
0312
0313 // map gdml tower ids to our columns
0314 int PHG4IHCalDetector::map_towerid(const int tower_id)
0315 {
0316   // odd id's go in one direction, even id's in the other one
0317   // this is a shortcut to derive the observed dependency
0318   // commented out after this code
0319   int itwr = -1;
0320   int itmp = tower_id / 2;
0321   if (tower_id % 2)
0322   {
0323     itwr = 12 + itmp;
0324   }
0325   else
0326   {
0327     itwr = 11 - itmp;
0328   }
0329   return itwr;
0330   // here is the mapping in long form
0331   // if (tower_id == 23)
0332   // {
0333   //   itwr = 0;
0334   // }
0335   // else if (tower_id == 21)
0336   // {
0337   //   itwr = 1;
0338   // }
0339   // else if (tower_id ==19 )
0340   // {
0341   //   itwr = 2;
0342   // }
0343   // else if (tower_id == 17)
0344   // {
0345   //   itwr = 3;
0346   // }
0347   // else if (tower_id == 15)
0348   // {
0349   //   itwr = 4;
0350   // }
0351   // else if (tower_id == 13)
0352   // {
0353   //   itwr = 5;
0354   // }
0355   // else if (tower_id == 11)
0356   // {
0357   //   itwr = 6;
0358   // }
0359   // else if (tower_id == 9)
0360   // {
0361   //   itwr = 7;
0362   // }
0363   // else if (tower_id == 7)
0364   // {
0365   //   itwr = 8;
0366   // }
0367   // else if (tower_id == 5)
0368   // {
0369   //   itwr = 9;
0370   // }
0371   // else if (tower_id == 3)
0372   // {
0373   //   itwr = 10;
0374   // }
0375   // else if (tower_id == 1)
0376   // {
0377   //   itwr = 11;
0378   // }
0379   // else if (tower_id == 0)
0380   // {
0381   //   itwr = 12;
0382   // }
0383   // else if (tower_id == 2)
0384   // {
0385   //   itwr = 13;
0386   // }
0387   // else if (tower_id == 4)
0388   // {
0389   //   itwr = 14;
0390   // }
0391   // else if (tower_id == 6)
0392   // {
0393   //   itwr = 15;
0394   // }
0395   // else if (tower_id == 8)
0396   // {
0397   //   itwr = 16;
0398   // }
0399   // else if (tower_id == 10)
0400   // {
0401   //   itwr = 17;
0402   // }
0403   // else if (tower_id == 12)
0404   // {
0405   //   itwr = 18;
0406   // }
0407   // else if (tower_id == 14)
0408   // {
0409   //   itwr = 19;
0410   // }
0411   // else if (tower_id == 16)
0412   // {
0413   //   itwr = 20;
0414   // }
0415   // else if (tower_id == 18)
0416   // {
0417   //   itwr = 21;
0418   // }
0419   // else if (tower_id == 20)
0420   // {
0421   //   itwr = 22;
0422   // }
0423   // else if (tower_id == 22)
0424   // {
0425   //   itwr = 23;
0426   // }
0427   // else
0428   // {
0429   //   std::cout << PHWHERE << " cannot map tower " << tower_id << std::endl;
0430   //   gSystem->Exit(1);
0431   //   exit(1);
0432   // }
0433 }
0434
0435 int PHG4IHCalDetector::map_layerid(const int layer_id)
0436 {
0437   int rowid = -1;
0438
0439   if (layer_id < 188)
0440   {
0441     rowid = layer_id + 68;
0442   }
0443   else  // (layer_id >= 188)
0444   {
0445     rowid = layer_id - 188;
0446   }
0447   // shift the row index up by 4
0448   rowid += 4;
0449   if (rowid > 255)
0450   {
0451     rowid -= 256;
0452   }
0453
0454   if (rowid > 255 || rowid < 0)
0455   {
0456     std::cout << PHWHERE << " row id out of range: " << rowid << std::endl;
0457     gSystem->Exit(1);
0458   }
0459   return rowid;
0460 }
0461
0462 // This is dulplicated code, we can get rid of it when we have the code to make towergeom for real data reco.
0463 void PHG4IHCalDetector::AddGeometryNode()
0464 {
0465   PHNodeIterator iter(topNode());
0466   PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0467   if (!runNode)
0468   {
0469     std::cout << PHWHERE << "Run Node missing, exiting." << std::endl;
0470     gSystem->Exit(1);
0471     exit(1);
0472   }
0473   PHNodeIterator runIter(runNode);
0474   PHCompositeNode *RunDetNode = dynamic_cast<PHCompositeNode *>(runIter.findFirst("PHCompositeNode", m_SuperDetector));
0475   if (!RunDetNode)
0476   {
0477     RunDetNode = new PHCompositeNode(m_SuperDetector);
0478     runNode->addNode(RunDetNode);
0479   }
0480   m_TowerGeomNodeName = "TOWERGEOM_" + m_SuperDetector;
0481   m_RawTowerGeom = findNode::getClass<RawTowerGeomContainer>(topNode(), m_TowerGeomNodeName);
0482   if (!m_RawTowerGeom)
0483   {
0484     m_RawTowerGeom = new RawTowerGeomContainer_Cylinderv1(RawTowerDefs::convert_name_to_caloid(m_SuperDetector));
0485     PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(m_RawTowerGeom, m_TowerGeomNodeName, "PHObject");
0486     RunDetNode->addNode(newNode);
0487   }
0488   double innerrad = m_Params->get_double_param(PHG4HcalDefs::innerrad);
0489   double thickness = m_Params->get_double_param(PHG4HcalDefs::outerrad) - innerrad;
0490   m_RawTowerGeom->set_radius(innerrad);
0491   m_RawTowerGeom->set_thickness(thickness);
0492   m_RawTowerGeom->set_phibins(m_Params->get_int_param(PHG4HcalDefs::n_towers));
0493   m_RawTowerGeom->set_etabins(m_Params->get_int_param("etabins"));
0494   double geom_ref_radius = innerrad + thickness / 2.;
0495   double phistart = m_Params->get_double_param("phistart");
0496   if (!std::isfinite(phistart))
0497   {
0498     std::cout << PHWHERE << " phistart is not finite: " << phistart
0499               << ", exiting now (this will crash anyway)" << std::endl;
0500     gSystem->Exit(1);
0501   }
0502   for (int i = 0; i < m_Params->get_int_param(PHG4HcalDefs::n_towers); i++)
0503   {
0504     double phiend = phistart + 2. * M_PI / m_Params->get_int_param(PHG4HcalDefs::n_towers);
0505     std::pair<double, double> range = std::make_pair(phiend, phistart);
0506     phistart = phiend;
0507     int tempi = i + 1;
0508     if (tempi >= m_Params->get_int_param(PHG4HcalDefs::n_towers))
0509     {
0510       tempi -= m_Params->get_int_param(PHG4HcalDefs::n_towers);
0511     }
0512     m_RawTowerGeom->set_phibounds(tempi, range);
0513   }
0514   double etalowbound = -m_Params->get_double_param("scinti_eta_coverage_neg");
0515   for (int i = 0; i < m_Params->get_int_param("etabins"); i++)
0516   {
0517     // double etahibound = etalowbound + 2.2 / get_int_param("etabins");
0518     double etahibound = etalowbound +
0519                         (m_Params->get_double_param("scinti_eta_coverage_neg") + m_Params->get_double_param("scinti_eta_coverage_pos")) / m_Params->get_int_param("etabins");
0520     std::pair<double, double> range = std::make_pair(etalowbound, etahibound);
0521     m_RawTowerGeom->set_etabounds(i, range);
0522     etalowbound = etahibound;
0523   }
0524   for (int iphi = 0; iphi < m_RawTowerGeom->get_phibins(); iphi++)
0525   {
0526     for (int ieta = 0; ieta < m_RawTowerGeom->get_etabins(); ieta++)
0527     {
0528       const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::convert_name_to_caloid(m_SuperDetector), ieta, iphi);
0529
0530       const double x(geom_ref_radius * std::cos(m_RawTowerGeom->get_phicenter(iphi)));
0531       const double y(geom_ref_radius * std::sin(m_RawTowerGeom->get_phicenter(iphi)));
0532       const double z(geom_ref_radius / std::tan(PHG4Utils::get_theta(m_RawTowerGeom->get_etacenter(ieta))));
0533
0534       RawTowerGeom *tg = m_RawTowerGeom->get_tower_geometry(key);
0535       if (tg)
0536       {
0537         if (Verbosity() > 0)
0538         {
0539           std::cout << "IHCalDetector::InitRun - Tower geometry " << key << " already exists" << std::endl;
0540         }
0541
0542         if (std::fabs(tg->get_center_x() - x) > 1e-4)
0543         {
0544           std::cout << "IHCalDetector::InitRun - Fatal Error - duplicated Tower geometry " << key << " with existing x = " << tg->get_center_x() << " and expected x = " << x
0545                     << std::endl;
0546
0547           return;
0548         }
0549         if (std::fabs(tg->get_center_y() - y) > 1e-4)
0550         {
0551           std::cout << "IHCalDetector::InitRun - Fatal Error - duplicated Tower geometry " << key << " with existing y = " << tg->get_center_y() << " and expected y = " << y
0552                     << std::endl;
0553           return;
0554         }
0555         if (std::fabs(tg->get_center_z() - z) > 1e-4)
0556         {
0557           std::cout << "IHCalDetector::InitRun - Fatal Error - duplicated Tower geometry " << key << " with existing z= " << tg->get_center_z() << " and expected z = " << z
0558                     << std::endl;
0559           return;
0560         }
0561       }
0562       else
0563       {
0564         if (Verbosity() > 0)
0565         {
0566           std::cout << "IHCalDetector::InitRun - building tower geometry " << key << "" << std::endl;
0567         }
0568
0569         tg = new RawTowerGeomv1(key);
0570
0571         tg->set_center_x(x);
0572         tg->set_center_y(y);
0573         tg->set_center_z(z);
0574         m_RawTowerGeom->add_tower_geometry(tg);
0575       }
0576     }
0577   }
0578   if (Verbosity() > 0)
0579   {
0580     m_RawTowerGeom->identify();
0581   }
0582 }