Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:18:15 

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!