diff --git a/Simulation/G4Utilities/Geo2G4/Geo2G4/Geo2G4Builder.h b/Simulation/G4Utilities/Geo2G4/Geo2G4/Geo2G4Builder.h old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/Geo2G4/Geo2G4MaterialFactory.h b/Simulation/G4Utilities/Geo2G4/Geo2G4/Geo2G4MaterialFactory.h old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/Geo2G4/Geo2G4SvcAccessor.h b/Simulation/G4Utilities/Geo2G4/Geo2G4/Geo2G4SvcAccessor.h old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/Geo2G4/Geo2G4SvcBase.h b/Simulation/G4Utilities/Geo2G4/Geo2G4/Geo2G4SvcBase.h old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/Geo2G4/IGeo2G4Svc.h b/Simulation/G4Utilities/Geo2G4/Geo2G4/IGeo2G4Svc.h old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/Geo2G4/LogicalVolume.h b/Simulation/G4Utilities/Geo2G4/Geo2G4/LogicalVolume.h old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/Geo2G4/VolumeBuilder.h b/Simulation/G4Utilities/Geo2G4/Geo2G4/VolumeBuilder.h old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/cmt/requirements b/Simulation/G4Utilities/Geo2G4/cmt/requirements index c39c52afe9710147892b653e471504009055e6ff..f7230870a8250903944997b766fda029fffcc838 100755 --- a/Simulation/G4Utilities/Geo2G4/cmt/requirements +++ b/Simulation/G4Utilities/Geo2G4/cmt/requirements @@ -2,8 +2,7 @@ package Geo2G4 author ADA -branches src cmt - +public use AtlasPolicy AtlasPolicy-* use AthenaKernel AthenaKernel-* Control use GaudiInterface GaudiInterface-* External @@ -19,6 +18,7 @@ use GeoSpecialShapes GeoSpecialShapes-* DetectorDescription/GeoModel use GeoModelInterfaces GeoModelInterfaces-* DetectorDescription/GeoModel use StoreGate StoreGate-* Control use SimHelpers SimHelpers-* Simulation/G4Sim +use SGTools SGTools-* Control end_private include_dirs "$(Geo2G4_root)" "$(Geo2G4_root)/Geo2G4" @@ -26,14 +26,15 @@ include_dirs "$(Geo2G4_root)" "$(Geo2G4_root)/Geo2G4" # Specify the required ROOT components for cmake (transparent to CMT) apply_pattern cmake_add_command command="find_package(ROOT COMPONENTS MathCore RIO)" -# Specify to cmake that this package has a non-standard include path (transparent to CMT) -apply_pattern cmake_add_command command="include_directories(Geo2G4)" - -# Override the library type for cmake so that it is linkable by clients (transparent to CMT) -apply_pattern cmake_override_library_type library=Geo2G4 type=installed_library - -library Geo2G4 ../src/*.cxx -s=components *.cxx +# ZLM to CMake folks: I believe that now that this is dual use and cleaner, these are not needed? +## Specify to cmake that this package has a non-standard include path (transparent to CMT) +#apply_pattern cmake_add_command command="include_directories(Geo2G4)" +# +## Override the library type for cmake so that it is linkable by clients (transparent to CMT) +#apply_pattern cmake_override_library_type library=Geo2G4 type=installed_library -apply_pattern component_library +apply_pattern dual_use_library files="*.cxx" -apply_pattern install_runtime +# Make dict for LArWheelSolidChecker +use AtlasReflex AtlasReflex-* External -no_auto_imports +apply_pattern lcgdict dict=LArWheelSolidChecker selectionfile=../src/lcg_dict/selection.xml headerfiles="../src/LArWheelSolidDDProxy.h" diff --git a/Simulation/G4Utilities/Geo2G4/src/ExtParameterisedVolumeBuilder.cxx b/Simulation/G4Utilities/Geo2G4/src/ExtParameterisedVolumeBuilder.cxx old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/ExtParameterisedVolumeBuilder.h b/Simulation/G4Utilities/Geo2G4/src/ExtParameterisedVolumeBuilder.h old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/GenericVolumeBuilder.cxx b/Simulation/G4Utilities/Geo2G4/src/GenericVolumeBuilder.cxx deleted file mode 100755 index f69e6b0cc5f22bb5cce00f7f736b63437e4a043e..0000000000000000000000000000000000000000 --- a/Simulation/G4Utilities/Geo2G4/src/GenericVolumeBuilder.cxx +++ /dev/null @@ -1,63 +0,0 @@ -/* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ - -#include "GenericVolumeBuilder.h" -#include "Geo2G4LogicalVolumeFactory.h" -#include "SimHelpers/ServiceAccessor.h" -#include "GeoModelKernel/GeoLogVol.h" -#include "AthenaBaseComps/AthMsgStreamMacros.h" -#include "G4PVPlacement.hh" -#include "G4ReflectionFactory.hh" -#include "globals.hh" -#include <iostream> - -LogicalVolume* GenericVolumeBuilder::Build(const PVConstLink theGeoPhysVolume, OpticalVolumesMap* /*optical_volumes*/) const -{ - static Geo2G4LogicalVolumeFactory LVFactory; - - const GeoLogVol * geoLog = theGeoPhysVolume->getLogVol(); - ATH_MSG_DEBUG ( " Start converting volume "<<geoLog->getName() ); - - G4LogicalVolume * theG4LogVolume = LVFactory.Build(geoLog); - - if (theGeoPhysVolume->getNChildVols()==0) return theG4LogVolume; - if (theG4LogVolume->GetNoDaughters() ) return theG4LogVolume; - // - // Loop over the children of the GeoVolume and place them - // - for(size_t ii = 0; ii<theGeoPhysVolume->getNChildVols(); ii++) - { - std::string nameChild = theGeoPhysVolume->getNameOfChildVol(ii); - // - // Get the id from GeoModel - Query<int> Qint = theGeoPhysVolume->getIdOfChildVol(ii); - int id = 90999; - if(Qint.isValid() ) id = Qint; - // - // Get the child Phys volume ii - // - PVConstLink theGeoPhysChild = theGeoPhysVolume->getChildVol(ii); - // - // Build the child - // - G4LogicalVolume* theG4LogChild = Build(theGeoPhysChild); - - // Get its transform - const G4Transform3D theG4Position(theGeoPhysVolume->getXToChildVol(ii)); - - if (nameChild == "ANON") nameChild=theG4LogChild->GetName(); - // log<<MSG::VERBOSE<<"\t Positioning "<<theG4LogChild->GetName()<< - // " into "<<theG4LogVolume->GetName()<< " with name "<<nameChild - // << "and Id:" << id <<endreq; - G4ReflectionFactory::Instance()->Place(theG4Position, - nameChild, - theG4LogChild, - theG4LogVolume, - false, - id); - - } - - return theG4LogVolume; -} diff --git a/Simulation/G4Utilities/Geo2G4/src/GenericVolumeBuilder.h b/Simulation/G4Utilities/Geo2G4/src/GenericVolumeBuilder.h deleted file mode 100755 index c253423773d79bb5a2d843abcf0365ae152818e6..0000000000000000000000000000000000000000 --- a/Simulation/G4Utilities/Geo2G4/src/GenericVolumeBuilder.h +++ /dev/null @@ -1,27 +0,0 @@ -/* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ - -#ifndef GenericVolumeBuilder_H -#define GenericVolumeBuilder_H - -#include "Geo2G4/VolumeBuilder.h" -#include "AthenaKernel/MsgStreamMember.h" - -class GenericVolumeBuilder: public VolumeBuilder -{ -public: - GenericVolumeBuilder(std::string n):VolumeBuilder(n),m_msg(n) - {} - - LogicalVolume* Build(PVConstLink pv, OpticalVolumesMap* optical_volumes = 0) const; - /// Log a message using the Athena controlled logging system - MsgStream& msg( MSG::Level lvl ) const { return m_msg << lvl; } - /// Check whether the logging system is active at the provided verbosity level - bool msgLvl( MSG::Level lvl ) const { return m_msg.get().level() <= lvl; } -private: - /// Private message stream member - mutable Athena::MsgStreamMember m_msg; -}; - -#endif diff --git a/Simulation/G4Utilities/Geo2G4/src/Geo2G4AssemblyFactory.cxx b/Simulation/G4Utilities/Geo2G4/src/Geo2G4AssemblyFactory.cxx old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/Geo2G4AssemblyFactory.h b/Simulation/G4Utilities/Geo2G4/src/Geo2G4AssemblyFactory.h old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/Geo2G4AssemblyVolume.cxx b/Simulation/G4Utilities/Geo2G4/src/Geo2G4AssemblyVolume.cxx old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/Geo2G4AssemblyVolume.h b/Simulation/G4Utilities/Geo2G4/src/Geo2G4AssemblyVolume.h old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/Geo2G4Builder.cxx b/Simulation/G4Utilities/Geo2G4/src/Geo2G4Builder.cxx old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/Geo2G4ElementFactory.cxx b/Simulation/G4Utilities/Geo2G4/src/Geo2G4ElementFactory.cxx old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/Geo2G4ElementFactory.h b/Simulation/G4Utilities/Geo2G4/src/Geo2G4ElementFactory.h old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/Geo2G4LVFactory.cxx b/Simulation/G4Utilities/Geo2G4/src/Geo2G4LVFactory.cxx old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/Geo2G4LVFactory.h b/Simulation/G4Utilities/Geo2G4/src/Geo2G4LVFactory.h old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/Geo2G4LogicalVolumeFactory.cxx b/Simulation/G4Utilities/Geo2G4/src/Geo2G4LogicalVolumeFactory.cxx old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/Geo2G4LogicalVolumeFactory.h b/Simulation/G4Utilities/Geo2G4/src/Geo2G4LogicalVolumeFactory.h old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/Geo2G4MatPropTableFactory.cxx b/Simulation/G4Utilities/Geo2G4/src/Geo2G4MatPropTableFactory.cxx old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/Geo2G4MatPropTableFactory.h b/Simulation/G4Utilities/Geo2G4/src/Geo2G4MatPropTableFactory.h old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/Geo2G4MaterialFactory.cxx b/Simulation/G4Utilities/Geo2G4/src/Geo2G4MaterialFactory.cxx old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/Geo2G4OpticalSurfaceFactory.cxx b/Simulation/G4Utilities/Geo2G4/src/Geo2G4OpticalSurfaceFactory.cxx old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/Geo2G4OpticalSurfaceFactory.h b/Simulation/G4Utilities/Geo2G4/src/Geo2G4OpticalSurfaceFactory.h old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/Geo2G4STParameterisation.cxx b/Simulation/G4Utilities/Geo2G4/src/Geo2G4STParameterisation.cxx old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/Geo2G4STParameterisation.h b/Simulation/G4Utilities/Geo2G4/src/Geo2G4STParameterisation.h old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/Geo2G4SolidFactory.cxx b/Simulation/G4Utilities/Geo2G4/src/Geo2G4SolidFactory.cxx old mode 100755 new mode 100644 index 5ea379ec56ef09d7c8a1757fd7918548706e9cae..398fb3fac9d1c2e9a192be8b7244d67c0bd461a6 --- a/Simulation/G4Utilities/Geo2G4/src/Geo2G4SolidFactory.cxx +++ b/Simulation/G4Utilities/Geo2G4/src/Geo2G4SolidFactory.cxx @@ -4,6 +4,7 @@ #include "Geo2G4SolidFactory.h" #include "LArWheelSolid.h" +#include "LArWheelSolidDDProxy.h" #include "GeoSpecialShapes/LArCustomShape.h" @@ -53,7 +54,59 @@ typedef std::map<const GeoShape*, G4VSolid*, std::less<const GeoShape*> > shapesMap; typedef std::map<std::string, G4VSolid*,std::less<std::string> > customSolidMap; -Geo2G4SolidFactory::Geo2G4SolidFactory(): m_msg("Geo2G4SolidFactory") +const Geo2G4SolidFactory::LArWheelSolid_typemap Geo2G4SolidFactory::s_lwsTypes = { + { "LAr::EMEC::InnerWheel::Absorber", {InnerAbsorberWheel, 1} }, + { "LAr::EMEC::Pos::InnerWheel::Absorber", {InnerAbsorberWheel, 1} }, + + { "LAr::EMEC::OuterWheel::Absorber", {OuterAbsorberWheel, 1} }, + { "LAr::EMEC::Pos::OuterWheel::Absorber", {OuterAbsorberWheel, 1} }, + + { "LAr::EMEC::InnerWheel::Electrode", {InnerElectrodWheel, 1} }, + { "LAr::EMEC::Pos::InnerWheel::Electrode", {InnerElectrodWheel, 1} }, + + { "LAr::EMEC::OuterWheel::Electrode", {OuterElectrodWheel, 1} }, + { "LAr::EMEC::Pos::OuterWheel::Electrode", {OuterElectrodWheel, 1} }, + + { "LAr::EMEC::Neg::InnerWheel::Absorber", {InnerAbsorberWheel, -1} }, + + { "LAr::EMEC::Neg::OuterWheel::Absorber", {OuterAbsorberWheel, -1} }, + + { "LAr::EMEC::Neg::InnerWheel::Electrode", {InnerElectrodWheel, -1} }, + + { "LAr::EMEC::Neg::OuterWheel::Electrode", {OuterElectrodWheel, -1} }, + + { "LAr::EMEC::InnerModule::Absorber", {InnerAbsorberModule, 1} }, + + { "LAr::EMEC::OuterModule::Absorber", {OuterAbsorberModule, 1} }, + + { "LAr::EMEC::InnerModule::Electrode", {InnerElectrodModule, 1} }, + + { "LAr::EMEC::OuterModule::Electrode", {OuterElectrodModule, 1} }, + + { "LAr::EMEC::InnerWheel::Glue", {InnerGlueWheel, 1} }, + { "LAr::EMEC::Pos::InnerWheel::Glue", {InnerGlueWheel, 1} }, + + { "LAr::EMEC::InnerWheel::Lead", {InnerLeadWheel, 1} }, + { "LAr::EMEC::Pos::InnerWheel::Lead", {InnerLeadWheel, 1} }, + + { "LAr::EMEC::OuterWheel::Glue", {OuterGlueWheel, 1} }, + { "LAr::EMEC::Pos::OuterWheel::Glue", {OuterGlueWheel, 1} }, + + { "LAr::EMEC::OuterWheel::Lead", {OuterLeadWheel, 1} }, + { "LAr::EMEC::Pos::OuterWheel::Lead", {OuterLeadWheel, 1} }, + + { "LAr::EMEC::Neg::InnerWheel::Glue", {InnerGlueWheel, -1} }, + + { "LAr::EMEC::Neg::OuterWheel::Glue", {OuterGlueWheel, -1} }, + + { "LAr::EMEC::Neg::InnerWheel::Lead", {InnerLeadWheel, -1} }, + + { "LAr::EMEC::Neg::OuterWheel::Lead", {OuterLeadWheel, -1} } +}; + +Geo2G4SolidFactory::Geo2G4SolidFactory() : + m_msg("Geo2G4SolidFactory"), + m_detStore( "StoreGateSvc/DetectorStore", "Geo2G4SolidFactory" ) { } @@ -431,47 +484,54 @@ G4VSolid *Geo2G4SolidFactory::Build(const GeoShape* geoShape, std::string name) else { theSolid = 0; - if(customName == "LAr::EMEC::InnerWheel::Absorber" || customName == "LAr::EMEC::Pos::InnerWheel::Absorber"){ - theSolid = new LArWheelSolid(customName, InnerAbsorberWheel, 1); - } else if(customName == "LAr::EMEC::OuterWheel::Absorber" || customName == "LAr::EMEC::Pos::OuterWheel::Absorber"){ - theSolid = new LArWheelSolid(customName, OuterAbsorberWheel, 1); - } else if(customName == "LAr::EMEC::InnerWheel::Electrode" || customName == "LAr::EMEC::Pos::InnerWheel::Electrode"){ - theSolid = new LArWheelSolid(customName, InnerElectrodWheel, 1); - } else if(customName == "LAr::EMEC::OuterWheel::Electrode" || customName == "LAr::EMEC::Pos::OuterWheel::Electrode"){ - theSolid = new LArWheelSolid(customName, OuterElectrodWheel, 1); - } else if(customName == "LAr::EMEC::Neg::InnerWheel::Absorber"){ - theSolid = new LArWheelSolid(customName, InnerAbsorberWheel, -1); - } else if(customName == "LAr::EMEC::Neg::OuterWheel::Absorber"){ - theSolid = new LArWheelSolid(customName, OuterAbsorberWheel, -1); - } else if(customName == "LAr::EMEC::Neg::InnerWheel::Electrode"){ - theSolid = new LArWheelSolid(customName, InnerElectrodWheel, -1); - } else if(customName == "LAr::EMEC::Neg::OuterWheel::Electrode"){ - theSolid = new LArWheelSolid(customName, OuterElectrodWheel, -1); - } else if(customName == "LAr::EMEC::InnerModule::Absorber"){ - theSolid = new LArWheelSolid(customName, InnerAbsorberModule, 1); - } else if(customName == "LAr::EMEC::OuterModule::Absorber"){ - theSolid = new LArWheelSolid(customName, OuterAbsorberModule, 1); - } else if(customName == "LAr::EMEC::InnerModule::Electrode"){ - theSolid = new LArWheelSolid(customName, InnerElectrodModule, 1); - } else if(customName == "LAr::EMEC::OuterModule::Electrode"){ - theSolid = new LArWheelSolid(customName, OuterElectrodModule, 1); - } else if(customName == "LAr::EMEC::InnerWheel::Glue" || customName == "LAr::EMEC::Pos::InnerWheel::Glue"){ - theSolid = new LArWheelSolid(customName, InnerGlueWheel, 1); - } else if(customName == "LAr::EMEC::InnerWheel::Lead" || customName == "LAr::EMEC::Pos::InnerWheel::Lead"){ - theSolid = new LArWheelSolid(customName, InnerLeadWheel, 1); - } else if(customName == "LAr::EMEC::OuterWheel::Glue" || customName == "LAr::EMEC::Pos::OuterWheel::Glue"){ - theSolid = new LArWheelSolid(customName, OuterGlueWheel, 1); - } else if(customName == "LAr::EMEC::OuterWheel::Lead" || customName == "LAr::EMEC::Pos::OuterWheel::Lead"){ - theSolid = new LArWheelSolid(customName, OuterLeadWheel, 1); - } else if(customName == "LAr::EMEC::Neg::InnerWheel::Glue"){ - theSolid = new LArWheelSolid(customName, InnerGlueWheel, -1); - } else if(customName == "LAr::EMEC::Neg::InnerWheel::Lead"){ - theSolid = new LArWheelSolid(customName, InnerLeadWheel, -1); - } else if(customName == "LAr::EMEC::Neg::OuterWheel::Glue"){ - theSolid = new LArWheelSolid(customName, OuterGlueWheel, -1); - } else if(customName == "LAr::EMEC::Neg::OuterWheel::Lead"){ - theSolid = new LArWheelSolid(customName, OuterLeadWheel, -1); - } +// if(customName == "LAr::EMEC::InnerWheel::Absorber" || customName == "LAr::EMEC::Pos::InnerWheel::Absorber"){ +// theSolid = new LArWheelSolid(customName, InnerAbsorberWheel, 1); +// } else if(customName == "LAr::EMEC::OuterWheel::Absorber" || customName == "LAr::EMEC::Pos::OuterWheel::Absorber"){ +// theSolid = new LArWheelSolid(customName, OuterAbsorberWheel, 1); +// } else if(customName == "LAr::EMEC::InnerWheel::Electrode" || customName == "LAr::EMEC::Pos::InnerWheel::Electrode"){ +// theSolid = new LArWheelSolid(customName, InnerElectrodWheel, 1); +// } else if(customName == "LAr::EMEC::OuterWheel::Electrode" || customName == "LAr::EMEC::Pos::OuterWheel::Electrode"){ +// theSolid = new LArWheelSolid(customName, OuterElectrodWheel, 1); +// } else if(customName == "LAr::EMEC::Neg::InnerWheel::Absorber"){ +// theSolid = new LArWheelSolid(customName, InnerAbsorberWheel, -1); +// } else if(customName == "LAr::EMEC::Neg::OuterWheel::Absorber"){ +// theSolid = new LArWheelSolid(customName, OuterAbsorberWheel, -1); +// } else if(customName == "LAr::EMEC::Neg::InnerWheel::Electrode"){ +// theSolid = new LArWheelSolid(customName, InnerElectrodWheel, -1); +// } else if(customName == "LAr::EMEC::Neg::OuterWheel::Electrode"){ +// theSolid = new LArWheelSolid(customName, OuterElectrodWheel, -1); +// } else if(customName == "LAr::EMEC::InnerModule::Absorber"){ +// theSolid = new LArWheelSolid(customName, InnerAbsorberModule, 1); +// } else if(customName == "LAr::EMEC::OuterModule::Absorber"){ +// theSolid = new LArWheelSolid(customName, OuterAbsorberModule, 1); +// } else if(customName == "LAr::EMEC::InnerModule::Electrode"){ +// theSolid = new LArWheelSolid(customName, InnerElectrodModule, 1); +// } else if(customName == "LAr::EMEC::OuterModule::Electrode"){ +// theSolid = new LArWheelSolid(customName, OuterElectrodModule, 1); +// } else if(customName == "LAr::EMEC::InnerWheel::Glue" || customName == "LAr::EMEC::Pos::InnerWheel::Glue"){ +// theSolid = new LArWheelSolid(customName, InnerGlueWheel, 1); +// } else if(customName == "LAr::EMEC::InnerWheel::Lead" || customName == "LAr::EMEC::Pos::InnerWheel::Lead"){ +// theSolid = new LArWheelSolid(customName, InnerLeadWheel, 1); +// } else if(customName == "LAr::EMEC::OuterWheel::Glue" || customName == "LAr::EMEC::Pos::OuterWheel::Glue"){ +// theSolid = new LArWheelSolid(customName, OuterGlueWheel, 1); +// } else if(customName == "LAr::EMEC::OuterWheel::Lead" || customName == "LAr::EMEC::Pos::OuterWheel::Lead"){ +// theSolid = new LArWheelSolid(customName, OuterLeadWheel, 1); +// } else if(customName == "LAr::EMEC::Neg::InnerWheel::Glue"){ +// theSolid = new LArWheelSolid(customName, InnerGlueWheel, -1); +// } else if(customName == "LAr::EMEC::Neg::InnerWheel::Lead"){ +// theSolid = new LArWheelSolid(customName, InnerLeadWheel, -1); +// } else if(customName == "LAr::EMEC::Neg::OuterWheel::Glue"){ +// theSolid = new LArWheelSolid(customName, OuterGlueWheel, -1); +// } else if(customName == "LAr::EMEC::Neg::OuterWheel::Lead"){ +// theSolid = new LArWheelSolid(customName, OuterLeadWheel, -1); +// } + + theSolid = createLArWheelSolid(customName, s_lwsTypes.at(customName) ); // map.at throws std::out_of_range exception on unknown shape name + if ( 0 == theSolid ) { + std::string error = std::string("Can't create LArWheelSolid for name ") + customName + " in Geo2G4SolidFactory::Build"; + throw std::runtime_error(error); + } + if(theSolid != 0) customSolids[customName] = theSolid; } } @@ -489,3 +549,19 @@ G4VSolid *Geo2G4SolidFactory::Build(const GeoShape* geoShape, std::string name) sharedShapes[geoShape] = theSolid; return theSolid; } + +G4VSolid* Geo2G4SolidFactory::createLArWheelSolid(const std::string& name, const LArWheelSolidDef_t & lwsdef) const { // LArWheelSolid_t wheelType, int zside + LArWheelSolid_t wheelType = lwsdef.first; + int zside = lwsdef.second; + + LArWheelSolid * theLWS = new LArWheelSolid(name, wheelType, zside); + + LArWheelSolidDDProxy * theLWS_p = new LArWheelSolidDDProxy(theLWS); + // ownership is passed to detStore + if ( detStore()->record(theLWS_p, name).isFailure() ) { + msg(MSG::WARNING) << "Can't store proxy for LArWheelSolid to the DetectorStore" <<endreq; + delete theLWS_p; + } + return theLWS; +} + diff --git a/Simulation/G4Utilities/Geo2G4/src/Geo2G4SolidFactory.h b/Simulation/G4Utilities/Geo2G4/src/Geo2G4SolidFactory.h old mode 100755 new mode 100644 index 1b92869247c2374acafb33b6e4aa57cab33a5f34..5caa38685e787420176f24cc293cc76ff4e65678 --- a/Simulation/G4Utilities/Geo2G4/src/Geo2G4SolidFactory.h +++ b/Simulation/G4Utilities/Geo2G4/src/Geo2G4SolidFactory.h @@ -7,7 +7,13 @@ #include <map> #include <string> + +#include "GaudiKernel/ServiceHandle.h" +//#include "GaudiKernel/StatusCode.h" + #include "AthenaKernel/MsgStreamMember.h" +#include "StoreGate/StoreGateSvc.h" +#include "LArWheelSolid_type.h" class G4VSolid; class GeoShape; @@ -15,15 +21,34 @@ class GeoShape; class Geo2G4SolidFactory { public: + typedef ServiceHandle<StoreGateSvc> StoreGateSvc_t; + typedef std::pair<LArWheelSolid_t, int> LArWheelSolidDef_t; + typedef std::map<std::string, LArWheelSolidDef_t> LArWheelSolid_typemap; + Geo2G4SolidFactory(); G4VSolid* Build(const GeoShape*, std::string name=std::string("")) const; /// Log a message using the Athena controlled logging system MsgStream& msg( MSG::Level lvl ) const { return m_msg << lvl; } /// Check whether the logging system is active at the provided verbosity level bool msgLvl( MSG::Level lvl ) const { return m_msg.get().level() <= lvl; } + + /** @brief The standard @c StoreGateSvc/DetectorStore + * Returns (kind of) a pointer to the @c StoreGateSvc + */ + StoreGateSvc_t& detStore() const; private: + G4VSolid* createLArWheelSolid(const std::string& name, const LArWheelSolidDef_t & lwsdef) const; + + static const LArWheelSolid_typemap s_lwsTypes; + /// Private message stream member mutable Athena::MsgStreamMember m_msg; + /// Pointer to StoreGate (detector store by default) + mutable StoreGateSvc_t m_detStore; }; +inline ServiceHandle<StoreGateSvc>& Geo2G4SolidFactory::detStore() const { + return m_detStore; +} + #endif diff --git a/Simulation/G4Utilities/Geo2G4/src/Geo2G4Svc.cxx b/Simulation/G4Utilities/Geo2G4/src/Geo2G4Svc.cxx old mode 100755 new mode 100644 index d067192b69c1fb6fe8648dc9656916a71fe89b48..3f08d6f379eebbaca62c14f63d92fd592c180522 --- a/Simulation/G4Utilities/Geo2G4/src/Geo2G4Svc.cxx +++ b/Simulation/G4Utilities/Geo2G4/src/Geo2G4Svc.cxx @@ -10,7 +10,7 @@ #include "StoreGate/StoreGateSvc.h" #include "Geo2G4/Geo2G4SvcAccessor.h" -#include "GenericVolumeBuilder.h" +#include "Geo2G4/VolumeBuilder.h" void InitializeBuilders(); diff --git a/Simulation/G4Utilities/Geo2G4/src/Geo2G4Svc.h b/Simulation/G4Utilities/Geo2G4/src/Geo2G4Svc.h old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/InitializeBuilders.cxx b/Simulation/G4Utilities/Geo2G4/src/InitializeBuilders.cxx old mode 100755 new mode 100644 index 6893ba48142f86ac25ef1f6156f0ac5676dcc049..57770d6c46345c439e9631e91220d3be1ff57a45 --- a/Simulation/G4Utilities/Geo2G4/src/InitializeBuilders.cxx +++ b/Simulation/G4Utilities/Geo2G4/src/InitializeBuilders.cxx @@ -2,15 +2,9 @@ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration */ -#include "GenericVolumeBuilder.h" -#include "ParameterisedVolumeBuilder.h" #include "ExtParameterisedVolumeBuilder.h" -#include "SingleLVCopyBuilder.h" void InitializeBuilders() { - GenericVolumeBuilder *gn __attribute__ ((unused)) = new GenericVolumeBuilder("Generic_Volume_Builder"); - ParameterisedVolumeBuilder *pb __attribute__ ((unused)) = new ParameterisedVolumeBuilder("Parameterised_Volume_Builder"); ExtParameterisedVolumeBuilder *epb __attribute__ ((unused)) = new ExtParameterisedVolumeBuilder("Extended_Parameterised_Volume_Builder"); - SingleLVCopyBuilder *sb __attribute__ ((unused)) = new SingleLVCopyBuilder("Single_LV_Copy_Builder"); } diff --git a/Simulation/G4Utilities/Geo2G4/src/LArFanSection.cxx b/Simulation/G4Utilities/Geo2G4/src/LArFanSection.cxx index f911d861678bdfcb65f6fb3036333f7314ecc76c..6a2b1a8591faa7560948fa4ba4e6e67cb6bf187f 100644 --- a/Simulation/G4Utilities/Geo2G4/src/LArFanSection.cxx +++ b/Simulation/G4Utilities/Geo2G4/src/LArFanSection.cxx @@ -7,295 +7,6 @@ #include "GeoSpecialShapes/LArWheelCalculator.h" #include<iostream> -#ifdef LArWheelSolidDTI_NEW -G4double LArWheelSolid::distance_to_in -#else -G4double LArWheelSolid::distance_to_in_ref -#endif -(G4ThreeVector &p, const G4ThreeVector &v) const -{ -#ifdef DEBUG_LARWHEELSOLID - if(Verbose) std::cout << "dti new1 " << p << " " << v << std::endl; -#endif - - G4double distance = 0.; - // p cannot be ouside on R due to bounding polycone limitation - if(p.x() > m_fs->xmax){ - if(v.x() >= 0.){ -#ifdef DEBUG_LARWHEELSOLID - if(Verbose) std::cout << "DTI inf1" << std::endl; -#endif - return kInfinity; - } - const G4double b = (m_fs->xmax - p.x()) / v.x(); - const G4double y2 = p.y() + v.y() * b; - const G4double z2 = p.z() + v.z() * b; - p.set(m_fs->xmax, y2, z2); - distance += b; - } else if(p.x() < m_fs->xmin){ - if(v.x() <= 0.){ -#ifdef DEBUG_LARWHEELSOLID - if(Verbose) std::cout << "DTI inf2" << std::endl; -#endif - return kInfinity; - } - const G4double b = (m_fs->xmin - p.x()) / v.x(); - const G4double y2 = p.y() + v.y() * b; - const G4double z2 = p.z() + v.z() * b; - p.set(m_fs->xmin, y2, z2); - distance += b; - } - G4double dist_p = Calculator->DistanceToTheNeutralFibre(p); - if(FanHalfThickness - fabs(dist_p) > Tolerance) return distance; - - size_t i0 = m_fs->select_section(p.z()); - - G4int step, start, stop; - if(v.z() < 0.){ step = -1; start = i0; stop = -1; } - else { step = 1; start = i0 + 1; stop = m_fs->z.size(); } - - static G4ThreeVector out_point; - G4int i; - G4double b = kInfinity; - bool inner_escape = fs_inner_escape(b, p, v); -#ifdef DEBUG_LARWHEELSOLID - if(Verbose){ - std::cout << "inner escape: " << (inner_escape? "yes ": "no"); - if(inner_escape) std::cout << b << " " << (p + v * b); - std::cout << std::endl; - } -#endif - const G4double &x2 = v.x() > 0.? m_fs->xmax: m_fs->xmin; -#ifdef DEBUG_LARWHEELSOLID - if(Verbose){ - std::cout << "from " << start << " to " << stop - << " step " << step << std::endl; - } -#endif - for(i = start; i != stop; i += step){ -#ifdef DEBUG_LARWHEELSOLID - if(Verbose){ - std::cout << "in dti at " << i << " " << p - << " d = " << distance << std::endl; - } -#endif - const G4double a = (m_fs->z[i] - p.z()) / v.z(); - const G4double x1 = p.x() + v.x() * a; - if(x1 < m_fs->xmin || x1 > m_fs->xmax){ // escape through sides - if(v.x() == 0.){ -#ifdef DEBUG_LARWHEELSOLID - if(Verbose) std::cout << "DTI INF" << std::endl; -#endif - return kInfinity;// parallel to OYZ, outside of x range => no intersection - } - const G4double b1 = (x2 - p.x()) / v.x(); - if(b1 < b) b = b1; -#ifdef DEBUG_LARWHEELSOLID - if(Verbose) std::cout << "b = " << b << std::endl; -#endif - break; - } - const G4double y1 = p.y() + v.y() * a; - if(y1 < 0.){ -#ifdef DEBUG_LARWHEELSOLID - if(Verbose) std::cout << "y escape" << std::endl; -#endif - break; // escape trough inner cone or a side - } - const G4double r1 = x1*x1 + y1*y1; - if(r1 < m_fs->rmin2[i] || r1 > m_fs->rmax2[i]){ -#ifdef DEBUG_LARWHEELSOLID - if(Verbose) std::cout << "R escape" << std::endl; -#endif - break; // escape trough one of cones or a side - } - out_point.set(x1, y1, m_fs->z[i]); - - const G4double d1 = (out_point - p).mag(); - if(inner_escape){ - if(distance + d1 > b){ -#ifdef DEBUG_LARWHEELSOLID - if(Verbose) std::cout << "inner escape " << std::endl; -#endif - break; - } - } - - G4ThreeVector p1(out_point); - G4double dist_out = Calculator->DistanceToTheNeutralFibre(out_point); - G4double dd = kInfinity; -#ifdef DEBUG_LARWHEELSOLID - if(Verbose){ - std::cout << ">" << p << " " << dist_p << " " << out_point - << " " << dist_out << std::endl; - } -#endif - if(dist_p * dist_out < 0.){// it certanly cross current half-wave - dd = in_iteration_process(p, dist_p, out_point); - p1 = p + v * dd; - } - G4double d2 = search_for_nearest_point(p, p1); - if(d2 < kInfinity){ - return distance + d2; // this half-wave is intersected - } else if(dd < kInfinity){ - return distance + dd; - } - distance += d1; - p = out_point; - dist_p = dist_out; - } - if(i == stop){ // escape trough z edge -#ifdef DEBUG_LARWHEELSOLID - if(Verbose) std::cout << "escape through z edge " << out_point << std::endl; -#endif - return kInfinity; //return distance; - } - - b -= distance; -#ifdef DEBUG_LARWHEELSOLID - if(Verbose) std::cout << "before fsoe: " << i << " " << b << std::endl; -#endif - fs_outer_escape(i, b, p, v); - out_point = p + v * b; - - G4double dist_out = Calculator->DistanceToTheNeutralFibre(out_point); -#ifdef DEBUG_LARWHEELSOLID - if(Verbose){ - std::cout << p << " " << dist_p << " " << out_point - << " " << dist_out << std::endl; - } -#endif - G4double dd = kInfinity; - if(dist_p * dist_out < 0.){// it certanly cross current half-wave - dd = in_iteration_process(p, dist_p, out_point); - out_point = p + v * dd; - } - G4double d2 = search_for_nearest_point(p, out_point); - if(d2 < kInfinity){ - return distance + d2; // this half-wave is intersected - } else if(dd < kInfinity){ - return distance + dd; - } - return kInfinity; -} - -#ifdef LArWheelSolidDTO_NEW -G4double LArWheelSolid::distance_to_out -#else -G4double LArWheelSolid::distance_to_out_ref -#endif -(const G4ThreeVector &p0, const G4ThreeVector &v) const -{ - static G4ThreeVector p, C; - static G4ThreeVector out_section, out, CP; - p = p0; - G4int fan_section = m_fs->select_section(p.z()); - if(v.z() > 0.) fan_section ++; - - G4double distance = 0.; - G4int step = v.z() < 0.? (-1): 1; - G4double x1; - G4double b = kInfinity; - bool inner_escape = fs_inner_escape(b, p, v); -#ifdef DEBUG_LARWHEELSOLID - if(Verbose){ - std::cout << "inner escape: " << (inner_escape? "yes ": "no"); - if(inner_escape) std::cout << b << " " << (p + v * b); - std::cout << std::endl; - } -#endif - do { - const G4double a = (m_fs->z[fan_section] - p.z()) / v.z(); - x1 = p.x() + v.x() * a; - if(x1 < m_fs->xmin || x1 > m_fs->xmax) break; - const G4double y1 = p.y() + v.y() * a; - if(y1 < 0.) break; - const G4double r1 = x1*x1 + y1*y1; - if(r1 > m_fs->rmax2[fan_section] || r1 < m_fs->rmin2[fan_section]) break; - - out_section.set(x1, y1, m_fs->z[fan_section]); -#ifdef DEBUG_LARWHEELSOLID - if(Verbose){ - std::cout << "fs " << fan_section << " " << p - << " " << out_section << std::endl; - } -#endif - const G4double d1 = (out_section - p).mag(); - if(inner_escape){ - if(distance + d1 > b){ -#ifdef DEBUG_LARWHEELSOLID - if(Verbose) std::cout << "inner escape " << std::endl; -#endif - break; - } - } - - G4double d = (out_section - p).mag(); - if(fabs(Calculator->DistanceToTheNeutralFibre(out_section)) - FanHalfThickness - > Tolerance) - { - d = out_iteration_process(p, out_section); - out_section = p + v * d; - //while(search_for_most_remoted_point(p, out_section, C)){ - if(search_for_most_remoted_point(p, out_section, C)){ - d = out_iteration_process(p, C); - out_section = p + v * d; - } - return distance + d; - } - if(search_for_most_remoted_point(p, out_section, C)){ - d = out_iteration_process(p, C); - return distance + d; - } - distance += d; - - if(fan_section == 0 || fan_section == m_fs->last_fs){ - return distance; - } - p = out_section; - fan_section += step; - } while(true); - - // for current fan_section exit point is certainly out of z edge - - // it's either on cones or on x size planes - - b -= distance; - const G4double &x2 = v.x() > 0.? m_fs->xmax: m_fs->xmin; - if(x1 > m_fs->xmax || x1 < m_fs->xmin){ - if(v.x() == 0.){ -#ifdef DEBUG_LARWHEELSOLID - if(Verbose) std::cout << "fs DTO parallel" << std::endl; -#endif - return 0;// parallel to OYZ, outside of x range => no intersection - } - G4double b1 = (x2 - p.x()) / v.x(); - if(b1 < b) b = b1; -#ifdef DEBUG_LARWHEELSOLID - if(Verbose) std::cout << "b = " << b << std::endl; -#endif - } - -#ifdef DEBUG_LARWHEELSOLID - if(Verbose) std::cout << "before fsoe: " << fan_section << " " << b << std::endl; -#endif - fs_outer_escape(fan_section, b, p, v); - out_section = p + v * b; - - if(fabs(Calculator->DistanceToTheNeutralFibre(out_section)) - FanHalfThickness - > Tolerance) - { - b = out_iteration_process(p, out_section); - out_section = p + v * b; - } - //while(search_for_most_remoted_point(out, out1, C)){ - if(search_for_most_remoted_point(p, out_section, C)){ - b = out_iteration_process(p, C); - out_section = p + v * b; - } - distance += b; - - return distance; -} void LArFanSections::print(void) const { @@ -303,135 +14,217 @@ void LArFanSections::print(void) const << "Amin = " << Amin << ", Amax = " << Amax << std::endl << "Bmin = " << Bmin << ", Bmax = " << Bmax << std::endl - << "xmin = " << xmin << ", xmax = " << xmax << std::endl - << "first_flat = " << first_flat; - if(Cflat2 > 0.) std::cout << ", Cflat2 is set " << Cflat2; - std::cout << std::endl; - std::cout << "Limits: (" << z.size() << ")" << std::endl; - for(size_t i = 0; i < z.size(); ++ i){ - std::cout << i << " " << z[i] << " " - << rmin2[i] << " " << rmax2[i] << std::endl; - } + << "xmin = " << xmin << ", xmax = " << xmax + << "Cflat2 = " << Cflat2 << std::endl; } LArFanSections::LArFanSections( - G4double Ain, G4double Aout, - G4double Bin, G4double Bout, - G4double Xmax - ){ - Cflat2=0.; - last_fs=0; - first_flat = -1; - Amin = Ain; Amax = Aout; - Bmin = Bin; Bmax = Bout; - Amin2 = Amin*Amin; Amax2 = Amax*Amax; - Bmin2 = Bmin*Bmin; Bmax2 = Bmax*Bmax; - ABmin = Amin*Bmin; ABmax = Amax*Bmax; - - xmax = Xmax; - xmin = -xmax; + G4double ri1, G4double ri2, G4double ro1, G4double ro2, + G4double Xmax, G4double z1, G4double z2 +) +{ + const G4double dz = z2 - z1; + Amin = (ri2 - ri1) / dz; + Amax = (ro2 - ro1) / dz; + Bmin = ri1 - Amin * z1; + Bmax = ro1 - Amax * z1; + Cflat2 = ro2*ro2; + + Amin2 = Amin*Amin; Amax2 = Amax*Amax; + Bmin2 = Bmin*Bmin; Bmax2 = Bmax*Bmax; + ABmin = Amin*Bmin; ABmax = Amax*Bmax; + + xmax = Xmax; + xmin = -xmax; } -size_t LArFanSections::select_section(const G4double &Z) +G4bool LArWheelSolid::check_D( + G4double &t1, G4double A, G4double B, G4double C, G4bool out +) const { - for(size_t i0 = last_fs - 1; i0 > 0; -- i0){ - if(Z > z[i0]) return i0; - } - return 0; + // G4bool out is to be set true if the point is surface-outside + // then have to discard first intersection + + const G4double D = B*B - A*C; + LWSDBG(8, std::cout << "check D=" << D << " out=" << out << std::endl); + if(D < 0.) return false; + + const G4double D1 = sqrt(D); + t1 = (-B + D1) / A; + const G4double t2 = (-B - D1) / A; + LWSDBG(8, std::cout << "t1=" << t1 << " t2=" << t2 << std::endl); + if(t1 > 0.){ + if(t2 > 0.){ + if(out){ + if(t2 > t1) t1 = t2; + } else { + if(t2 < t1) t1 = t2; + } + } else { + if(out) return false; + } + } else { + if(t2 > 0.){ + if(out) return false; + t1 = t2; + } else { + return false; + } + } + return true; } -bool LArWheelSolid::fs_inner_escape(G4double &b, - const G4ThreeVector &p, const G4ThreeVector &v - ) const +// p must be not outside of the "FanBound" +// if track crosses inner cone in valid (z, x) interval, +// returns true, sets q to the cross point +bool LArWheelSolid::fs_cross_lower( + const G4ThreeVector &p, const G4ThreeVector &v, + G4ThreeVector &q) const { - const G4double A = v.perp2() - m_fs->Amin2*v.z()*v.z(); - const G4double B = p.x()*v.x() + p.y()*v.y() - - m_fs->Amin2*p.z()*v.z() - m_fs->ABmin*v.z(); - const G4double C = p.perp2() - - m_fs->Amin2*p.z()*p.z() - - 2.*m_fs->ABmin*p.z() - m_fs->Bmin2; - const G4double D = B*B - A*C; - if(D > 0.){ - const G4double &zmin = m_fs->z.front(); - const G4double &zmax = m_fs->z.back(); - const G4double D1 = sqrt(D); - G4double t1 = (-B + D1) / A; - const G4double zz1 = p.z() + v.z() * t1; - if(zz1 < zmin || zz1 > zmax) t1 = kInfinity; - if(t1 > -Tolerance){ - const G4double t2 = (-B - D1) / A; - const G4double zz2 = p.z() + v.z() * t2; - if(!(zz2 < zmin || zz2 > zmax) && t2 > -Tolerance) - { - if(t1 > Tolerance){ - if(t2 > Tolerance){ - b = t2 < t1? t2: t1; - return true; - } else { - b = t2; - return true; - } - } else { - if(t2 > Tolerance){ - b = t1; - return true; - }// else point on surface, v just touching the cone, let's keep b - } - } // else keep b - } // else keep b - } - return false; + LWSDBG(7, std::cout << "fcl" << std::endl); + const G4double A = v.perp2() - m_fs->Amin2*v.z()*v.z(); + const G4double B = p.x()*v.x() + p.y()*v.y() + - m_fs->Amin2*p.z()*v.z() - m_fs->ABmin*v.z(); + const G4double C = p.perp2() - m_fs->Amin2*p.z()*p.z() + - 2.*m_fs->ABmin*p.z() - m_fs->Bmin2; + G4double t1(0.0); + const G4double out_dist = m_fs->Amin*p.z() + m_fs->Bmin - p.perp(); + LWSDBG(8, std::cout << "fcl out_dist(p)=" << out_dist << " Tolerance=" << Tolerance << std::endl); + const G4bool out = out_dist >= 0.0; + if(check_D(t1, A, B, C, out)){ + const G4double zz1 = p.z() + v.z() * t1; + if(zz1 < Zsect.front() || zz1 > Zsect.back()){ + LWSDBG(8, std::cout << "fcl out on Z " << zz1 << std::endl); + return false; + } + const G4double xx1 = p.x() + v.x() * t1; + if(xx1 < m_fs->xmin || xx1 > m_fs->xmax){ + LWSDBG(8, std::cout << "fcl out on X " << xx1 << std::endl); + return false; + } + q.setX(xx1); + q.setY(p.y() + v.y() * t1); + q.setZ(zz1); + LWSDBG(8, std::cout << "fcl got " << t1 << std::endl); + return true; + } + LWSDBG(8, std::cout << "fcl no intersection" << std::endl); + return false; } -void LArWheelSolid::fs_outer_escape( - const G4int &fan_section, G4double &b, - const G4ThreeVector &p, const G4ThreeVector &v) const +// p must be not outside of the "FanBound" +// if track crosses outer cone in valid (z, x) interval, +// returns true, adds to b the distance to the cross point, +// sets q to the cross point +bool LArWheelSolid::fs_cross_upper( + const G4ThreeVector &p, const G4ThreeVector &v, + G4ThreeVector &q) const { - G4double A = v.perp2(); - G4double B = p.x()*v.x() + p.y()*v.y(); - G4double C = p.perp2(); - if(fan_section < m_fs->first_flat){ // cone part - A -= m_fs->Amax2*v.z()*v.z(); - B -= m_fs->Amax2*p.z()*v.z() + m_fs->ABmax*v.z(); - C -= m_fs->Amax2*p.z()*p.z() + 2.*m_fs->ABmax*p.z() + m_fs->Bmax2; - } else { // cyliner part - C -= m_fs->Cflat2; - } - const G4double D = B*B - A*C; - if(D > 0.){ - const G4double D1 = sqrt(D); - const G4double t1 = (-B + D1) / A; - const G4double t2 = (-B - D1) / A; - if(t1 < -Tolerance){ - if(t2 >= -Tolerance && b > t2){ - b = t2; - } - } else if(t1 < Tolerance){ - if(t2 < -Tolerance){ - if(b > t1) b = t1; - } else if(t2 >= Tolerance && b > t2) b = t2; - } else {//if t1 > Tolerance - if(t2 < Tolerance){ - if(b > t1) b = t1; - } else { - const G4double t = t1 > t2? t2: t1; - if(b > t) b = t; - } - } - } + LWSDBG(7, std::cout << "fcu" << std::endl); + G4double A = v.perp2(); + G4double B = p.x()*v.x() + p.y()*v.y(); + G4double C = p.perp2(); + + if(IsOuter){ + const G4double &Af = A, &Bf = B; + const G4double Cf = C - m_fs->Cflat2; + G4double b1; + if(check_D(b1, Af, Bf, Cf, Cf >= 0.)){ + const G4double zz1 = p.z() + v.z() * b1; + if(zz1 >= Zmid && zz1 <= Zsect.back()){ + const G4double xx1 = p.x() + v.x() * b1; + if(xx1 < m_fs->xmin || xx1 > m_fs->xmax) return false; + q.setX(xx1); + q.setY(p.y() + v.y() * b1); + q.setZ(zz1);; + return true; + } + } + LWSDBG(8, std::cout << "fcu no cyl intersection" << std::endl); + } + + A -= m_fs->Amax2*v.z()*v.z(); + B -= m_fs->Amax2*p.z()*v.z() + m_fs->ABmax*v.z(); + C -= m_fs->Amax2*p.z()*p.z() + 2.*m_fs->ABmax*p.z() + m_fs->Bmax2; + + G4double t1; + const G4bool out = m_fs->Amax*p.z() + m_fs->Bmax <= p.perp(); + if(check_D(t1, A, B, C, out)){ + const G4double zz1 = p.z() + v.z() * t1; + if(zz1 < Zsect.front() || zz1 > Zmid) return false; + const G4double xx1 = p.x() + v.x() * t1; + if(xx1 < m_fs->xmin || xx1 > m_fs->xmax) return false; + q.setX(xx1); + q.setY(p.y() + v.y() * t1); + q.setZ(zz1); + return true; + } + LWSDBG(8, std::cout << "fcu no cone intersection" << std::endl); + return false; } -void LArFanSections::prepare(void) +/* p must be not outside "FanBound" */ +LArWheelSolid::FanBoundExit_t LArWheelSolid::find_exit_point( + const G4ThreeVector &p, const G4ThreeVector &v, + G4ThreeVector &q) const { - Cflat2 = 0.; - for(first_flat = 1; first_flat < G4int(z.size()); ++ first_flat){ - if(fabs(rmax2[first_flat] - rmax2[first_flat - 1]) < 0.00001){ - Cflat2 = rmax2[first_flat]; - break; - } - } - last_fs = z.size() - 1; - Amin2 = Amin*Amin; Amax2 = Amax*Amax; - Bmin2 = Bmin*Bmin; Bmax2 = Bmax*Bmax; - ABmin = Amin*Bmin; ABmax = Amax*Bmax; + LWSDBG(6, std::cout << "in fep p" << MSG_VECTOR(p)<< ", v"<< MSG_VECTOR(v) << ", q" << MSG_VECTOR(q) << std::endl); + +/* by construction, cannot have true from both upper and lower */ +/* the only problem is the points on surface but "slightly outside" */ +/* fs_cross_* account for (x, z) range */ +// lower has to be checked first, since outer might find more distant +// intersection in the acceptable (x, z) range + if(fs_cross_lower(p, v, q)) return ExitAtInner; + LWSDBG(6, std::cout << "after fs_cross_lower q" << MSG_VECTOR(q) << std::endl); + if(fs_cross_upper(p, v, q)) return ExitAtOuter; + LWSDBG(6, std::cout << "after fs_cross_upper q" << MSG_VECTOR(q) << std::endl); + + FanBoundExit_t result = ExitAtSide; + G4double d; + if(v.x() > 0.) d = (m_fs->xmax - p.x()) / v.x(); + else if(v.x() < 0.) d = (m_fs->xmin - p.x()) / v.x(); + else d = kInfinity; + + G4double dz; + FanBoundExit_t resultz = NoCross; + if(v.z() > 0.){ + dz = (Zsect.back() - p.z()) / v.z(); + resultz = ExitAtBack; + } else if(v.z() < 0.){ + dz = (Zsect.front() - p.z()) / v.z(); + resultz = ExitAtFront; + } else { + dz = kInfinity; + } + if(d > dz){ + d = dz; + result = resultz; + } + q = p + v * d; + LWSDBG(7, std::cout << "fep side " << d << " " << result << " q" << MSG_VECTOR(q) << std::endl); + const G4double out_distlower = m_fs->Amin*q.z() + m_fs->Bmin - q.perp(); // > 0 - below lower cone + LWSDBG(7, std::cout << "fep out_distlower(q)=" << out_distlower << " Tolerance=" << Tolerance << std::endl); + if (out_distlower >= 0.0) { + // side intersection point is below lower cone + // initial point p was at exit boundary + q = p; + return NoCross; + } + + if (IsOuter && q.z() >= Zmid && q.z() <= Zsect.back()+Tolerance && q.perp2() >= m_fs->Cflat2) { + // outside of upper cylinder + q = p; + return NoCross; + } + const G4double out_distupper = m_fs->Amax*q.z() + m_fs->Bmax - q.perp(); // < 0 - above upper cone + LWSDBG(7, std::cout << "fep out_distupper(q)=" << out_distupper << " Tolerance=" << Tolerance << std::endl); + if (out_distupper <= 0.0) { + // side intersection point is above upper cone + // initial point p was at exit boundary + q = p; + return NoCross; + } + assert((q - p).mag() < kInfinity); + return result; } diff --git a/Simulation/G4Utilities/Geo2G4/src/LArFanSection.h b/Simulation/G4Utilities/Geo2G4/src/LArFanSection.h index 738f5fe74839c67f2faf48dafbe1f2960d5d5ad7..48b0f645f466590bfb56bd8843f9698d8ff26baf 100644 --- a/Simulation/G4Utilities/Geo2G4/src/LArFanSection.h +++ b/Simulation/G4Utilities/Geo2G4/src/LArFanSection.h @@ -8,41 +8,22 @@ // helper class to replace G4Polycone // in certain LArWheelSolid operations -#include<vector> - -#include "G4VSolid.hh" - -#define LAR_FAN_SECTIONS_MULT 1 - struct LArFanSections { - G4double Amin, Amax; - G4double Bmin, Bmax; - G4double Amin2, Amax2; - G4double Bmin2, Bmax2; - G4double xmin, xmax; - G4double Cflat2, ABmax, ABmin; - G4int last_fs, first_flat; - std::vector<G4double> z; - std::vector<G4double> rmin2; - std::vector<G4double> rmax2; - - LArFanSections( - G4double Ain, G4double Aout, - G4double Bin, G4double Bout, - G4double Xmax - ); - - ~LArFanSections() - { - z.clear(); - rmin2.clear(); - rmax2.clear(); - } - - void prepare(void); - void print(void) const; - size_t select_section(const G4double &); + G4double Amin, Amax; + G4double Bmin, Bmax; + G4double Amin2, Amax2; + G4double Bmin2, Bmax2; + G4double xmin, xmax; + G4double Cflat2, ABmax, ABmin; + + LArFanSections( + G4double ri1, G4double ri2, + G4double ro1, G4double ro2, + G4double Xmax, G4double z1, G4double z2 + ); + + void print(void) const; }; #endif // __LARFANSECTION_H__ diff --git a/Simulation/G4Utilities/Geo2G4/src/LArWheelSolid.cxx b/Simulation/G4Utilities/Geo2G4/src/LArWheelSolid.cxx old mode 100755 new mode 100644 index 168c96c9282692877d297bc5ec8c9cd42f9d1c10..5f783272a56df03f4ab7717edbaae364cdc87e5c --- a/Simulation/G4Utilities/Geo2G4/src/LArWheelSolid.cxx +++ b/Simulation/G4Utilities/Geo2G4/src/LArWheelSolid.cxx @@ -2,72 +2,64 @@ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration */ -// LArWheelSolid.cc (see also LArWheelSolid*.cc) -// Authors: A. M. Soukharev, with revisions by William G. Seligman -// 16-Oct-2001 WGS: Incorporated corrected for first and last wave -// from Jozsef Toth. -// Nov 2001 AMS: bugfixes and updates of algorithms -// also in files LArWheelSolid*.cc -// 01-Apr-2002 AMS: Added GetPhiGap routine. -// August 2002: speed improvements -// 15-Apr-2003 AMS: moved geometry code to LArWheelCalculator - #include "G4VGraphicsScene.hh" #include "G4VisExtent.hh" -#include "G4Polycone.hh" -#include "G4ios.hh" #include "GeoSpecialShapes/LArWheelCalculator.h" #include "LArWheelSolid.h" -class G4Polyhedron; class G4NURBS; class G4VoxelLimits; class G4AffineTransform; EInside LArWheelSolid::Inside(const G4ThreeVector &inputP) const { - static G4ThreeVector p; - EInside inside_bp = BoundingPolycone->Inside(inputP); - if(inside_bp == kOutside) return kOutside; - p = inputP; - G4double d = FanHalfThickness - fabs(Calculator->DistanceToTheNearestFan(p)); - if(d > Tolerance) return inside_bp; - if(d > -Tolerance) return kSurface; - return kOutside; + LWSDBG(10, std::cout << std::setprecision(25)); + LWSDBG(1, std::cout << TypeStr() << " Inside " << MSG_VECTOR(inputP) << std::endl); + static G4ThreeVector p; + const EInside inside_BS = BoundingShape->Inside(inputP); + if(inside_BS == kOutside){ + LWSDBG(2, std::cout << "outside BS" << std::endl); + return kOutside; + } + p = inputP; + const G4double d = fabs(Calculator->DistanceToTheNearestFan(p)); + if(d > FHTplusT){ + LWSDBG(2, std::cout << "outside fan d=" << d << ", FHTplusT=" << FHTplusT << std::endl); + return kOutside; + } + if(d < FHTminusT){ + LWSDBG(2, std::cout << "inside fan d=" << d << ", FHTminusT=" << FHTminusT << ", inside_BS=" << inside(inside_BS) << std::endl); + return inside_BS; + } + LWSDBG(2, std::cout << "surface" << std::endl); + return kSurface; } -G4ThreeVector LArWheelSolid::SurfaceNormal (const G4ThreeVector &inputP) const +G4ThreeVector LArWheelSolid::SurfaceNormal(const G4ThreeVector &inputP) const { - static G4ThreeVector p, d; - EInside inside_bp = BoundingPolycone->Inside(inputP); - if(inside_bp == kOutside - || (inside_bp == kSurface && Inside(inputP) == kSurface)) - { - return BoundingPolycone->SurfaceNormal(inputP); - } - p = inputP; - Calculator->DistanceToTheNearestFan(p); - d = Calculator->NearestPointOnNeutralFibre(p); - d.rotateZ(inputP.phi() - p.phi()); // rotate back to initial position - p = inputP - d; - return(p.unit()); -} - -G4int LArWheelSolid::select_fan_section(G4double z) const -{ - G4int i; - for(i = 1; i <= MaxFanSectionLimits; i ++){ - if(FanSectionLimits[i] >= z) return i - 1; - } - return i - 2; + LWSDBG(1, std::cout << TypeStr() << " SurfaceNormal" << MSG_VECTOR(inputP) << std::endl); + static G4ThreeVector p, d; + EInside inside_BS = BoundingShape->Inside(inputP); + if(inside_BS != kInside){ + LWSDBG(2, std::cout << "not inside BS" << std::endl); + return BoundingShape->SurfaceNormal(inputP); + } + p = inputP; + Calculator->DistanceToTheNearestFan(p); + d = Calculator->NearestPointOnNeutralFibre(p); + d.rotateZ(inputP.phi() - p.phi()); // rotate back to initial position + LWSDBG(4, std::cout << "npnf" << MSG_VECTOR(d) << std::endl); + p = inputP - d; + LWSDBG(2, std::cout << "sn " << MSG_VECTOR(p.unit()) << std::endl); + return(p.unit()); } G4bool LArWheelSolid::CalculateExtent(const EAxis a, const G4VoxelLimits &vl, const G4AffineTransform &t, G4double &p, G4double &q) const { - return BoundingPolycone->CalculateExtent(a, vl, t, p, q); + return BoundingShape->CalculateExtent(a, vl, t, p, q); } G4GeometryType LArWheelSolid::GetEntityType() const @@ -121,10 +113,21 @@ void LArWheelSolid::DescribeYourselfTo(G4VGraphicsScene &scene) const G4VisExtent LArWheelSolid::GetExtent() const { - return BoundingPolycone->GetExtent(); + return BoundingShape->GetExtent(); } G4Polyhedron* LArWheelSolid::CreatePolyhedron() const { - return BoundingPolycone->CreatePolyhedron(); + return BoundingShape->CreatePolyhedron(); +} + +/* + * returns the number of lower z boundary of z-section containing Z + */ +G4int LArWheelSolid::select_section(const G4double &Z) const +{ + for(G4int i = Zsect_start_search; i > 0; -- i){ + if(Z > Zsect[i]) return i; + } + return 0; } diff --git a/Simulation/G4Utilities/Geo2G4/src/LArWheelSolid.h b/Simulation/G4Utilities/Geo2G4/src/LArWheelSolid.h old mode 100755 new mode 100644 index 21d72c757af275b09d8ab499c9a09d65258b23cc..4e27149ade5de5a06403d794dad19c05ce3b6f9d --- a/Simulation/G4Utilities/Geo2G4/src/LArWheelSolid.h +++ b/Simulation/G4Utilities/Geo2G4/src/LArWheelSolid.h @@ -2,16 +2,8 @@ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration */ -// LArWheelSolid -// Authors: A. M. Soukharev, with revisions by William Seligman - -// Revised 11-May-2001 WGS to integrate it with LArEMECConstruction -// November 2001 updates and bugfixes by A. Soukharev -// August 2002 - A. Soukharev - speed improvements -// 15 Apr 2003 AMS move geometry calculations to Calculator - -#ifndef __LArWheelSolid_HH__ -#define __LArWheelSolid_HH__ +#ifndef __LArWheelSolid_H__ +#define __LArWheelSolid_H__ #include "AthenaKernel/MsgStreamMember.h" #include "G4VSolid.hh" @@ -20,15 +12,30 @@ // disabled by default to avoid any performance degradation //#define DEBUG_LARWHEELSOLID -// set this to use speed-improved version of DistanceToOut -#define LArWheelSolidDTO_NEW +// set this to use native G4 FanBound's methods for DisToIn +// instead of local calculations +//#define LARWHEELSOLID_USE_FANBOUND + +// set this to use BoundingShape's methods for DisToOut +// instead of local calculations +//#define LARWHEELSOLID_USE_BS_DTO + +// change this to have more z sections +#define LARWHEELSOLID_ZSECT_MULT 1 + -// set this to use speed-improved version of DistanceToIn -#define LArWheelSolidDTI_NEW +// set this to check in dti and dto functions if particle direction +// pointing inside or outside of volume to return zero fast when it is required by spec. +// currently at development stage, requires accurate surface normal calculations +//#define CHECK_DIRTONORM_ANGLE_ON_SURFACE #ifdef DEBUG_LARWHEELSOLID -extern G4bool Verbose; +#define LWSDBG(a, b) if(Verbose > a) b #define MSG_VECTOR(v) "(" << v.x() << ", " << v.y() << ", " << v.z() << ")" +#define LWS_HARD_TEST_DTI +#define LWS_HARD_TEST_DTO +#else +#define LWSDBG(a, b) #endif // Forward declarations. @@ -42,21 +49,9 @@ class G4Polycone; class LArWheelCalculator; class TF1; class LArFanSections; +class G4Polyhedra; -typedef enum { - InnerAbsorberWheel, - OuterAbsorberWheel, - InnerElectrodWheel, - OuterElectrodWheel, - InnerAbsorberModule, - OuterAbsorberModule, - InnerElectrodModule, - OuterElectrodModule, - InnerGlueWheel, - OuterGlueWheel, - InnerLeadWheel, - OuterLeadWheel -} LArWheelSolid_t; +#include "LArWheelSolid_type.h" inline const char *LArWheelSolidTypeString(LArWheelSolid_t type) { @@ -119,25 +114,34 @@ public: private: static const G4double Tolerance; + static const G4double AngularTolerance; static const G4double IterationPrecision; static const G4double IterationPrecision2; static const unsigned int IterationsLimit; - LArWheelSolid_t Type; + G4bool IsOuter; + const LArWheelSolid_t Type; LArWheelCalculator *Calculator; - G4double FanHalfThickness; + G4double FanHalfThickness, FHTplusT, FHTminusT; G4double FanPhiAmplitude; G4double MinPhi; G4double MaxPhi; - G4double PhiPosition; + const G4double PhiPosition; G4Polycone* BoundingPolycone; - G4Polycone** FanSection; - G4int MaxFanSection; - G4double *FanSectionLimits; - G4int MaxFanSectionLimits; + G4VSolid* BoundingShape; +#ifdef LARWHEELSOLID_USE_FANBOUND + G4VSolid* FanBound; +#endif + + std::vector<G4double> Zsect; + G4int Zsect_start_search; LArFanSections *m_fs; + // z at outer wheel "bend" + G4double Zmid; + // Special limit, used in dto + G4double Ymin; // limits for use in service functions G4double Zmin, Zmax, Rmin, Rmax; @@ -145,30 +149,34 @@ private: void outer_solid_init(const G4String &); void set_phi_size(void); - virtual G4double distance_to_in(const G4ThreeVector &, - const G4double, - const G4ThreeVector &) const; virtual G4double distance_to_in(G4ThreeVector &, const G4ThreeVector &) const; - virtual G4double distance_to_in_ref(G4ThreeVector &, const G4ThreeVector &) const; G4double in_iteration_process(const G4ThreeVector &, - G4double, - const G4ThreeVector &) const; - G4double search_for_nearest_point(const G4ThreeVector &, - const G4ThreeVector &) const; + G4double, G4ThreeVector &) const; + G4double search_for_nearest_point( + const G4ThreeVector &, const G4double, + const G4ThreeVector & + ) const; G4bool search_for_most_remoted_point(const G4ThreeVector &, const G4ThreeVector &, G4ThreeVector &) const; - G4double distance_to_out(const G4ThreeVector &, - const G4ThreeVector &) const; - G4double distance_to_out_ref(const G4ThreeVector &, const G4ThreeVector &) const; G4double out_iteration_process(const G4ThreeVector &, - const G4ThreeVector &) const; - G4int select_fan_section(G4double) const; - - bool fs_inner_escape(G4double &b, - const G4ThreeVector &p, const G4ThreeVector &v) const; - void fs_outer_escape(const G4int &fan_section, G4double &b, - const G4ThreeVector &p, const G4ThreeVector &v) const; + G4ThreeVector &) const; + + G4bool fs_cross_lower(const G4ThreeVector &p, const G4ThreeVector &v, + G4ThreeVector &q) const; + G4bool fs_cross_upper(const G4ThreeVector &p, const G4ThreeVector &v, + G4ThreeVector &q) const; + typedef enum { + NoCross, ExitAtInner, ExitAtOuter, + ExitAtFront, ExitAtBack, ExitAtSide + } FanBoundExit_t; + FanBoundExit_t find_exit_point(const G4ThreeVector &p, + const G4ThreeVector &v, + G4ThreeVector &q) const; + G4bool check_D(G4double &b, + G4double A, G4double B, G4double C, G4bool) const; + + G4int select_section(const G4double &Z) const; EInside Inside_accordion(const G4ThreeVector&) const; void get_point_on_accordion_surface(G4ThreeVector &) const; @@ -205,10 +213,6 @@ protected: friend double LArWheelSolid_fcn_side_area(double *, double *); #ifdef DEBUG_LARWHEELSOLID - G4double in_chord_method( - const G4ThreeVector &p0, const G4ThreeVector &p1, - const G4ThreeVector &v) const; - static const char* inside(EInside i) { switch(i){ @@ -218,7 +222,17 @@ protected: } return "unknown"; } + + public: + static G4int Verbose; + void SetVerbose(G4int v){ Verbose = v; } + G4bool test_dti(const G4ThreeVector &p, + const G4ThreeVector &v, const G4double distance) const; + G4bool test_dto(const G4ThreeVector &p, + const G4ThreeVector &v, const G4double distance) const; + private: + const char *TypeStr(void) const { return LArWheelSolidTypeString(Type); } #endif }; -#endif // __LArWheelSolid_HH__ +#endif // __LArWheelSolid_H__ diff --git a/Simulation/G4Utilities/Geo2G4/src/LArWheelSolidDDProxy.cxx b/Simulation/G4Utilities/Geo2G4/src/LArWheelSolidDDProxy.cxx new file mode 100644 index 0000000000000000000000000000000000000000..450291f4bb2f29c2f3f0cd5fef311b6389ad6336 --- /dev/null +++ b/Simulation/G4Utilities/Geo2G4/src/LArWheelSolidDDProxy.cxx @@ -0,0 +1,65 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +// LArWheelSolidDDProxy +// proxy for LArWheelSolid to store in DetectorDescription +// Author: D. A. Maximov + +#include "LArWheelSolidDDProxy.h" +#include "LArWheelSolid.h" +#include<iostream> + +LArWheelSolidDDProxy::LArWheelSolidDDProxy(LArWheelSolid* plws) : + m_plws(plws) + {} + +LArWheelSolidDDProxy::~LArWheelSolidDDProxy() {} + +int LArWheelSolidDDProxy::Inside(const CLHEP::Hep3Vector& p) const { + return m_plws->Inside(p); +} + +double LArWheelSolidDDProxy::DistanceToIn(const CLHEP::Hep3Vector& p, const CLHEP::Hep3Vector& v) const { + return m_plws->DistanceToIn(p, v); +} + +double LArWheelSolidDDProxy::DistanceToIn(const CLHEP::Hep3Vector& p) const { + return m_plws->DistanceToIn(p); +} + +double LArWheelSolidDDProxy::DistanceToOut(const CLHEP::Hep3Vector& p, const CLHEP::Hep3Vector& v) const { + return m_plws->DistanceToOut(p, v); +} + +double LArWheelSolidDDProxy::DistanceToOut(const CLHEP::Hep3Vector& p) const { + return m_plws->DistanceToOut(p); +} + +CLHEP::Hep3Vector LArWheelSolidDDProxy::SurfaceNormal(const CLHEP::Hep3Vector& p) const { + return m_plws->SurfaceNormal(p); +} + +CLHEP::Hep3Vector LArWheelSolidDDProxy::GetPointOnSurface() const { + return m_plws->GetPointOnSurface(); +} + +double LArWheelSolidDDProxy::GetCubicVolume() { + return m_plws->GetCubicVolume(); +} + +double LArWheelSolidDDProxy::GetSurfaceArea() { + return m_plws->GetSurfaceArea(); +} + +#ifdef DEBUG_LARWHEELSOLID +void LArWheelSolidDDProxy::SetVerbose(int v) const +{ + m_plws->SetVerbose(v); +} +#else +void LArWheelSolidDDProxy::SetVerbose(int) const +{ + std::cerr << "DEBUG_LARWHEELSOLID is off" << std::endl; +} +#endif diff --git a/Simulation/G4Utilities/Geo2G4/src/LArWheelSolidDDProxy.h b/Simulation/G4Utilities/Geo2G4/src/LArWheelSolidDDProxy.h new file mode 100644 index 0000000000000000000000000000000000000000..52a8ba2830ac93efb32cc3a571afa1ec60cfd655 --- /dev/null +++ b/Simulation/G4Utilities/Geo2G4/src/LArWheelSolidDDProxy.h @@ -0,0 +1,68 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +// LArWheelSolidDDProxy +// proxy for LArWheelSolid to store in DetectorDescription +// Author: D. A. Maximov + +#ifndef __LArWheelSolidDDProxy_HH__ +#define __LArWheelSolidDDProxy_HH__ + +// #include "AthenaKernel/MsgStreamMember.h" +#include "CLHEP/Vector/ThreeVector.h" +#include "SGTools/CLASS_DEF.h" + + +class LArWheelSolid; + +class LArWheelSolidDDProxy { +public: + + LArWheelSolidDDProxy(LArWheelSolid* plws); + virtual ~LArWheelSolidDDProxy(); + + // Mandatory for custom solid Geant4 functions +/* EInside Inside(const G4ThreeVector&) const; + G4double DistanceToIn(const G4ThreeVector&, + const G4ThreeVector&) const; + G4double DistanceToIn(const G4ThreeVector&) const; + G4double DistanceToOut(const G4ThreeVector&, + const G4ThreeVector&, + const G4bool calcNorm = false, + G4bool* validNorm = 0, + G4ThreeVector* n = 0) const; */ + int Inside(const CLHEP::Hep3Vector&) const; + + double DistanceToIn(const CLHEP::Hep3Vector&, const CLHEP::Hep3Vector&) const; + double DistanceToIn(const CLHEP::Hep3Vector&) const; + + double DistanceToOut(const CLHEP::Hep3Vector&, const CLHEP::Hep3Vector&) const; + + double DistanceToOut(const CLHEP::Hep3Vector&) const; + CLHEP::Hep3Vector SurfaceNormal (const CLHEP::Hep3Vector&) const; + +// G4bool CalculateExtent(const EAxis, +// const G4VoxelLimits&, +// const G4AffineTransform&, +// G4double&, +// G4double&) const; + + CLHEP::Hep3Vector GetPointOnSurface() const; + double GetCubicVolume(); + double GetSurfaceArea(); + + void SetVerbose(int v) const; + +private: + + LArWheelSolid * m_plws; + +}; + + +//using the macro below we can assign an identifier (and a version) +//This is required and checked at compile time when you try to record/retrieve +CLASS_DEF(LArWheelSolidDDProxy, 900345679 , 1) + +#endif // __LArWheelSolidDDProxy_HH__ diff --git a/Simulation/G4Utilities/Geo2G4/src/LArWheelSolidDisToIn.cxx b/Simulation/G4Utilities/Geo2G4/src/LArWheelSolidDisToIn.cxx old mode 100755 new mode 100644 index 7e9aa75e08b70a340247a49680219ec6b06fec68..d5f74a9c981328f6b59506aa1d44dde6a8004224 --- a/Simulation/G4Utilities/Geo2G4/src/LArWheelSolidDisToIn.cxx +++ b/Simulation/G4Utilities/Geo2G4/src/LArWheelSolidDisToIn.cxx @@ -5,417 +5,303 @@ // DistanceToIn stuff for LArWheelSolid #include <cassert> #include "AthenaBaseComps/AthMsgStreamMacros.h" -#include "G4Polycone.hh" #include "CLHEP/Units/PhysicalConstants.h" #include "GeoSpecialShapes/LArWheelCalculator.h" #include "LArWheelSolid.h" +#include "LArFanSection.h" G4double LArWheelSolid::DistanceToIn(const G4ThreeVector &inputP) const { - static G4ThreeVector p; - if(BoundingPolycone->Inside(inputP) == kOutside){ - // here is an approximation - for the point being outside BoundingPolycone - // the solid looks like solid BoundingPolycone - return BoundingPolycone->DistanceToIn(inputP); - } - p = inputP; - G4double d = fabs(Calculator->DistanceToTheNearestFan(p)) - FanHalfThickness; - if(d <= Tolerance) return 0.; - return d; + LWSDBG(1, std::cout << TypeStr() << " DisToIn" << MSG_VECTOR(inputP) << std::endl); + if(BoundingShape->Inside(inputP) == kOutside) { + // here is an approximation - for the point outside BoundingShape + // the solid looks like a BoundingShape + // it's okay since the result could be underestimated + LWSDBG(2, std::cout << "Outside BS" << std::endl); + return BoundingShape->DistanceToIn(inputP); + } + G4ThreeVector p(inputP); + const G4double d = fabs(Calculator->DistanceToTheNearestFan(p)); + if(d > FHTplusT) return d - FanHalfThickness; + LWSDBG(2, std::cout << "already inside" << MSG_VECTOR(p) << std::endl); + return 0.; } -// inputV should be unit G4double LArWheelSolid::DistanceToIn(const G4ThreeVector &inputP, const G4ThreeVector &inputV) const { - static G4ThreeVector p, v; - G4double distance = 0.; - if(BoundingPolycone->Inside(inputP) == kOutside){ - distance = BoundingPolycone->DistanceToIn(inputP, inputV); - if(distance == kInfinity) return kInfinity; - p = inputP + inputV * distance; - assert(BoundingPolycone->Inside(p) != kOutside); - } else p = inputP; - G4double phi0 = p.phi(); - G4double d = Calculator->DistanceToTheNearestFan(p); - if(FanHalfThickness - fabs(d) > Tolerance) return distance; - v = inputV; - v.rotateZ(p.phi() - phi0); -#ifndef DEBUG_LARWHEELSOLID - return distance + distance_to_in(p, v); -#else - G4ThreeVector p1(p), v1(v), p0(p), v0(v); - G4double dd1 = distance_to_in_ref(p1, v1); - G4double dd = distance_to_in(p, v); - if(fabs(dd - dd1) > IterationPrecision){ - static int cnt = 0; - std::cout << "DTI " << p0 << " " << v0 << " " - << LArWheelSolidTypeString(Type) << std::endl - << "DTI_ref " << dd1 << " DTI " << dd << std::endl; - std::cout << "DTI MISMATCH " << (dd - dd1) << std::endl; - Verbose = true; - p1 = p0; p = p0; v1 = v0; v = v0; - G4double dd2 = distance_to_in_ref(p1, v1); - G4double dd3 = distance_to_in(p, v); - std::cout << " === " << dd2 << " === " << dd3 << " === " << std::endl; - cnt ++; - // if(cnt > 10) exit(0); - Verbose = false; - if(cnt < 10){ - FILE *F = fopen("test_dti", "a"); - fwrite(&Type, sizeof(Type), 1, F); - fwrite(&p0, sizeof(p0), 1, F); - fwrite(&v0, sizeof(v0), 1, F); - fclose(F); - } - } - bool inf = false; - if(dd > 10.*m){ - dd = BoundingPolycone->DistanceToOut(p0, v); - inf = true; - } - size_t i = 1; - for(G4double t = 0.; t < dd; t += IterationPrecision * 10.){ - G4ThreeVector b = p0 + v * t; - if(FanHalfThickness - fabs(Calculator->DistanceToTheNeutralFibre(b)) > Tolerance){ - if(fabs(t - dd) <= IterationPrecision) break; - std::cout << "@# " << Type << " " << p << " " << v << " "; - std::cout << i << " " << (inf? kInfinity: dd) << " " << t << std::endl; - break; - } - i ++; - } - if(inf) return kInfinity; - return(distance + dd); -#endif + LWSDBG(1, std::cout << TypeStr() << " DisToIn" << MSG_VECTOR(inputP) + << MSG_VECTOR(inputV) << std::endl); + + G4double distance = 0.; + const EInside inside_BS = BoundingShape->Inside(inputP); + G4ThreeVector p(inputP); + if(inside_BS == kOutside) { + distance = BoundingShape->DistanceToIn(inputP, inputV); + if(distance == kInfinity) { + LWSDBG(2, std::cout << "Infinity distance to BoundingShape" << MSG_VECTOR(inputP) << MSG_VECTOR(inputV) << std::endl); + return kInfinity; + } + p += inputV * distance; + assert(BoundingShape->Inside(p) != kOutside); + LWSDBG(2, std::cout << "shift" << MSG_VECTOR(inputP) << std::endl); + } + + const G4double phi0 = p.phi(); + const G4double d = Calculator->DistanceToTheNearestFan(p); + if(fabs(d) < FHTminusT){ + LWSDBG(2, std::cout << "already inside fan" << MSG_VECTOR(p) << std::endl); + // if initial point is on BS surface and inside fan volume it is a real surface + if(inside_BS == kSurface) { + LWSDBG(2, std::cout << "On BS surface" << std::endl); + return BoundingShape->DistanceToIn(inputP, inputV); + } + return distance; + } + G4ThreeVector v(inputV); + v.rotateZ(p.phi() - phi0); + const G4double d0 = distance_to_in(p, v); + +#ifdef DEBUG_LARWHEELSOLID + if(Verbose > 2){ + if(Verbose > 3){ + std::cout << MSG_VECTOR(inputP) + << " " << MSG_VECTOR(inputV) << std::endl; + } + std::cout << "DTI: " << d0 << std::endl; + } + //distance += d0; + if(Verbose > 3){ + if(distance < kInfinity){ + G4ThreeVector p2 = inputP + inputV * (distance+d0); + EInside i = Inside(p2); + std::cout << "DTI hit at " << MSG_VECTOR(p2) << ", " + << inside(i) << std::endl; + } + } +#ifdef LWS_HARD_TEST_DTI + if(test_dti(inputP, inputV, distance+d0 )){ + if(Verbose == 1){ + std::cout << TypeStr() << " DisToIn" << MSG_VECTOR(inputP) + << MSG_VECTOR(inputV) << std::endl; + } + } + if(Verbose == 1){ + std::cout << TypeStr() << " DisToIn" << MSG_VECTOR(inputP) + << MSG_VECTOR(inputV) << " " << distance+d0 << std::endl; + } +#endif // ifdef LWS_HARD_TEST_DTI + +#endif // ifdef DEBUG_LARWHEELSOLID + + return distance+d0; +//#else +// return distance + distance_to_in(p, v); } // This functions should be called in the case when we are sure that // points p (which sould be OUTSIDE of vertical fan) and out_point have // the surface of the vertical fan between them. // returns distance from point p to absorber surface -G4double LArWheelSolid::in_iteration_process(const G4ThreeVector &p, - G4double dist_p, - const G4ThreeVector &out_point) - const +// sets last parameter to the founded point +G4double LArWheelSolid::in_iteration_process( + const G4ThreeVector &p, G4double dist_p, G4ThreeVector &B +) const { -#ifdef DEBUG_LARWHEELSOLID - if(Verbose){ - std::cout << "iip from " << p << " to " << out_point - << " dir " << (out_point - p).unit() << std::endl; - } -#endif - static G4ThreeVector A, B, C, diff; - A = p; - B = out_point; - G4double dist_c; - unsigned int niter = 0; - // assert(fabs(Calculator->DistanceToTheNeutralFibre(A)) - FanHalfThickness > Tolerance); + LWSDBG(6, std::cout << "iip from " << p << " to " << B + << " dir " << (B - p).unit() + << std::endl); + + static G4ThreeVector A, C, diff; + A = p; + G4double dist_c; + unsigned int niter = 0; + // assert(fabs(Calculator->DistanceToTheNeutralFibre(A)) > FHTplusT); // assert(Calculator->DistanceToTheNeutralFibre(A) == dist_p); - do { - C = A + B; - C *= 0.5; - dist_c = Calculator->DistanceToTheNeutralFibre(C); - if(dist_c * dist_p < 0 - || (FanHalfThickness - fabs(dist_c) > Tolerance)){ - B = C; - } else { - A = C; - } - niter ++; - diff = A - B; - } while(diff.mag2() > IterationPrecision2 && niter < IterationsLimit); - assert(niter < IterationsLimit); - assert(fabs(Calculator->DistanceToTheNeutralFibre(B)) - FanHalfThickness < Tolerance); - /*std::cout << "diff^2 " << diff.mag2() << " " << IterationPrecision2; - std::cout << " dist_c " << dist_c << " " << (FanHalfThickness - fabs(dist_c)) << " tol " << Tolerance - << std::endl;*/ - diff = p - B; -#ifdef DEBUG_LARWHEELSOLID - if(Verbose){ - std::cout << "iip result in " << niter << " = " << B - << " " << diff.mag() << std::endl; - } -#endif - return diff.mag(); + do { + C = A + B; + C *= 0.5; + dist_c = Calculator->DistanceToTheNeutralFibre(C); + if(dist_c * dist_p < 0. || fabs(dist_c) < FHTminusT){ + B = C; + } else { + A = C; + } + niter ++; + diff = A - B; + } while(diff.mag2() > IterationPrecision2 && niter < IterationsLimit); + assert(niter < IterationsLimit); + assert(fabs(Calculator->DistanceToTheNeutralFibre(B)) < FHTplusT); + diff = p - B; + LWSDBG(7, std::cout << "iip result in " << niter << " = " << B + << " " << diff.mag() << std::endl); + return diff.mag(); } // search for the nearest to the neutral fibre of the vertical fan point // on the segment between p_in and p_out -G4double LArWheelSolid::search_for_nearest_point(const G4ThreeVector &p_in, - const G4ThreeVector &p_out) const +G4double LArWheelSolid::search_for_nearest_point( + const G4ThreeVector &p_in, const G4double dist_p_in, + const G4ThreeVector &p_out +) const { - static G4ThreeVector A, B, C, l, diff; - A = p_in; - B = p_out; - diff = B - A; - l = diff.unit() * IterationPrecision; + LWSDBG(6, std::cout << "sfnp " << MSG_VECTOR(p_in) << " " + << MSG_VECTOR(p_out) << std::endl); + + static G4ThreeVector A, B, C, l, diff; + A = p_in; + B = p_out; + diff = B - A; + l = diff.unit() * IterationPrecision; // this is to correctly take the sign of the distance into account - G4double dist_p_in = Calculator->DistanceToTheNeutralFibre(p_in); - G4double sign = dist_p_in < 0.? -1. : 1.; - G4double d_prime; - G4double dist_c; - unsigned long niter = 0; - do{ - C = A + B; - C *= 0.5; - dist_c = Calculator->DistanceToTheNeutralFibre(C); - if(dist_c * sign <= 0.){ // we are in coditions for in_iteration_process - return in_iteration_process(p_in, dist_p_in, C); - } + G4double sign = dist_p_in < 0.? -1. : 1.; + G4double d_prime; + G4double dist_c; + unsigned long niter = 0; + do { + C = A + B; + C *= 0.5; + dist_c = Calculator->DistanceToTheNeutralFibre(C); + if(dist_c * sign <= 0.){ // we are in coditions for in_iteration_process + LWSDBG(7, std::cout << "sfnp0 " << dist_c << std::endl); + return in_iteration_process(p_in, dist_p_in, C); + } // calculate sign of derivative of distance to the neutral fibre - // d_prime = (Calculator->DistanceToTheNeutralFibre(C + l) - - // Calculator->DistanceToTheNeutralFibre(C - l)) * sign; // hope this substitution is acceptable - diff = C - l; - d_prime = (dist_c - Calculator->DistanceToTheNeutralFibre(diff)) * sign; - if(d_prime < 0.) A = C; - else B = C; - niter ++; - diff = A - B; - } while(diff.mag2() > IterationPrecision2 && niter < IterationsLimit); - /* if(niter >= IterationsLimit){ - G4cout << p_in << " " << p_out << G4endl; - G4cout << A << " " << B << " " << diff.mag() << " " - << dist_p_in << " " << dist_c << " " << sign << " " - << G4endl; - } - */ assert(niter < IterationsLimit); - if(FanHalfThickness - fabs(dist_c) > Tolerance){ - return in_iteration_process(p_in, dist_p_in, C); - } + diff = C - l; + d_prime = (dist_c - Calculator->DistanceToTheNeutralFibre(diff)) * sign; + if(d_prime < 0.) A = C; + else B = C; + niter ++; + diff = A - B; + } while(diff.mag2() > IterationPrecision2 && niter < IterationsLimit); + assert(niter < IterationsLimit); + if(fabs(dist_c) < FHTminusT){ + LWSDBG(7, std::cout << "sfnp1 " << dist_c << std::endl); + return in_iteration_process(p_in, dist_p_in, C); + } // let's check at p_in and p_out - if(dist_p_in * sign < dist_c * sign){ - C = p_in; - dist_c = dist_p_in; - } - G4double dist_p_out = Calculator->DistanceToTheNeutralFibre(p_out); - if(dist_p_out *sign < dist_c * sign) C = p_out; - if(FanHalfThickness - fabs(dist_p_out) > Tolerance){ - return in_iteration_process(p_in, dist_p_in, C); - } - return kInfinity; + if(dist_p_in * sign < dist_c * sign){ + C = p_in; + dist_c = dist_p_in; + } + G4double dist_p_out = Calculator->DistanceToTheNeutralFibre(p_out); + if(dist_p_out *sign < dist_c * sign) C = p_out; + if(fabs(dist_p_out) < FHTminusT){ + LWSDBG(7, std::cout << "sfnp2 " << dist_p_out << std::endl); + return in_iteration_process(p_in, dist_p_in, C); + } + return kInfinity; } -#ifdef LArWheelSolidDTI_NEW -G4double LArWheelSolid::distance_to_in_ref -#else -G4double LArWheelSolid::distance_to_in -#endif -(G4ThreeVector &p, const G4ThreeVector &v) const +G4double LArWheelSolid::distance_to_in(G4ThreeVector &p, const G4ThreeVector &v) const { - static G4ThreeVector out_point; -#ifdef DEBUG_LARWHEELSOLID - G4ThreeVector p0 = p; - if(Verbose) std::cout << "dti old" << std::endl; -#endif - G4int sign, start, stop; - if(v[2] < 0.){ sign = -1; start = MaxFanSection; stop = -1; } - else { sign = 1; start = 0; stop = MaxFanSection + 1; } - G4double distance = 0., result = kInfinity; - G4double dist_p = 0.; - for(G4int i = start; i != stop; i += sign){ -#ifdef DEBUG_LARWHEELSOLID - if(Verbose){ - std::cout << "dti old at " << i << " " << p << " d = " - << distance << std::endl; - } -#endif - if(FanSection[i]->Inside(p) == kOutside){ - G4double d1 = FanSection[i]->DistanceToIn(p, v); - //if(d1 < kInfinity){ // crossed this FS - // I guess here is again bug in G4Polycone. Assertion below failed sometimes. - // Let's try to work it around. - if(d1 < 10.*CLHEP::m){ // crossed this FS - distance += d1; - p += v * d1; - } else { -#ifdef DEBUG_LARWHEELSOLID - if(d1 < kInfinity){ - ATH_MSG_ERROR ( "DistanceToIn(p, v):" << G4endl - << "FanSection[" << i << "]->DistanceToIn(" << MSG_VECTOR(p) - << ", " << MSG_VECTOR(v) << ") = " << d1 ); - } -#endif - continue; - } - } - assert(FanSection[i]->Inside(p) != kOutside); - dist_p = Calculator->DistanceToTheNeutralFibre(p); - if(FanHalfThickness - fabs(dist_p) > Tolerance){ - result = distance; - break; - } + LWSDBG(4, std::cout << "dti: " << MSG_VECTOR(p) << " " + << MSG_VECTOR(v) << std::endl); - G4double dist_from_p_to_out = FanSection[i]->DistanceToOut(p, v); - - /* workaround for bug in G4Polycone */ - if(dist_from_p_to_out > 10 * CLHEP::m){ -#ifdef DEBUG_LARWHEELSOLID - ATH_MSG_WARNING ("distance_to_in: workaround for a 'bug' in G4Polycone " - << " at " << MSG_VECTOR(p) << ": " - << dist_from_p_to_out ); -#endif - return kInfinity; - } - - out_point = p + v * dist_from_p_to_out; - - G4double dist_out = Calculator->DistanceToTheNeutralFibre(out_point); -#ifdef DEBUG_LARWHEELSOLID - if(Verbose){ - std::cout << ">" << p << " " << dist_p << " " << out_point - << " " << dist_out << std::endl - << "out_point r " << out_point.perp() << " phi " - << out_point.phi() << std::endl; - } -#endif - if(dist_p * dist_out < 0.){// it certanly cross current half-wave - result = distance + in_iteration_process(p, dist_p, out_point); - break; - } - G4double d2 = search_for_nearest_point(p, out_point); - if(d2 < kInfinity){ - result = distance + d2; // this half-wave is intersected - break; - } - p = out_point; - distance += dist_from_p_to_out; - } -#ifdef DEBUG_LARWHEELSOLID - if(Verbose) std::cout << "dti old end " << result << std::endl; + G4double distance = 0.; +#ifdef LARWHEELSOLID_USE_FANBOUND + if(FanBound->Inside(p) == kOutside) { + const G4double d = FanBound->DistanceToIn(p, v); + p += v * d; + distance += d; + } +#else + if(p.x() > m_fs->xmax) { + if(v.x() >= 0.) return kInfinity; + const G4double b = (m_fs->xmax - p.x()) / v.x(); + const G4double y2 = p.y() + v.y() * b; + const G4double z2 = p.z() + v.z() * b; + p.set(m_fs->xmax, y2, z2); + distance += b; + } else if(p.x() < m_fs->xmin) { + if(v.x() <= 0.) return kInfinity; + const G4double b = (m_fs->xmin - p.x()) / v.x(); + const G4double y2 = p.y() + v.y() * b; + const G4double z2 = p.z() + v.z() * b; + p.set(m_fs->xmin, y2, z2); + distance += b; + } #endif - return result; -} +// here p is on surface of or inside the "FanBound", +// distance corrected, misses are accounted for + LWSDBG(5, std::cout << "dti corrected: " << MSG_VECTOR(p) << std::endl); -// this is obsolete version of distance_to_in -// it is kept for information -// calculates distance from point p along vector v to the nearest surface of -// the vertical fan -G4double LArWheelSolid::distance_to_in(const G4ThreeVector &p, - const G4double dist_p, - const G4ThreeVector &v) const -{ - static G4ThreeVector out_point, A, B; - static G4int warning = 0; + G4double dist_p = Calculator->DistanceToTheNeutralFibre(p); + if(fabs(dist_p) < FHTminusT) { + LWSDBG(5, std::cout << "hit fan dist_p=" << dist_p << ", FHTminusT=" << FHTminusT << std::endl); + return distance; + } - G4int fan_section = select_fan_section(p.z()); - assert(fan_section >= 0); - assert(fan_section <= MaxFanSection); - G4Polycone *current_section = FanSection[fan_section]; - // assert(current_section->Inside(p) != kOutside); - - if(current_section->Inside(p) == kOutside){ -#ifdef DEBUG_LARWHEELSOLID - ATH_MSG_WARNING ( "#" << (warning ++) << " from istance_to_in:" << G4endl - << "current_section->Inside(p) == kOutside!!!" << G4endl - << " point: " << MSG_VECTOR(p) << ", phi = " << p.phi() - << ", r = " << p.r() - << " DTI: " << current_section->DistanceToIn(p) - << ", DTNeaFan: " << dist_p - << " fs parametes:" << G4endl - << "\tphi borders: " << current_section->GetStartPhi() - << ", " << current_section->GetEndPhi() - ); +#ifdef CHECK_DIRTONORM_ANGLE_ON_SURFACE + if(fabs(dist_p) > FHTplusT) { + LWSDBG(5, std::cout << "outside fan dist_p=" << dist_p << ", FHTplusT=" << FHTplusT << std::endl); + } else { + LWSDBG(5, std::cout << "on fan surface dist_p=" << dist_p << ", FHTplusT=" << FHTplusT << ", FHTminusT=" << FHTminusT << std::endl); + + const G4ThreeVector d = Calculator->NearestPointOnNeutralFibre(p); + // check dot product between normal and v + if ( (p-d).cosTheta(v) < -AngularTolerance ) { + // direction "v" definitely pointing inside + // return 0.0, it should be in "distance" + return distance; + } + } #endif - //FIXME: do we need to increment warning even if DEBUG_LARWHEELSOLID is not defined? - if(warning >= 5){ - ATH_MSG_WARNING ( "iterations are stopped" ); - return kInfinity; - } - G4double d = current_section->DistanceToIn(p, v); - ATH_MSG_WARNING ( " DTI(p, v) = " << d ); - if(d < kInfinity){ - G4ThreeVector p1 = p + d * v; - d += distance_to_in(p1, dist_p, v); - } - return d; - } - warning = 0; - /**/ - // determine where vector v cross current polycone section - G4double dist_from_p_to_out = current_section->DistanceToOut(p, v); - out_point = p + v * dist_from_p_to_out; - G4double dist_out = Calculator->DistanceToTheNeutralFibre(out_point); - if(dist_p * dist_out < 0){// it certanly cross current half-wave - return in_iteration_process(p, dist_p, out_point); - } - G4double distance = search_for_nearest_point(p, out_point); - if(distance < kInfinity) return distance; // this half-wave is intersected - // another section or no intersection at all - G4double z_min = FanSectionLimits[fan_section] + Tolerance; - G4double z_max = FanSectionLimits[fan_section + 1] - Tolerance; - if((out_point.z() >= z_max && fan_section < MaxFanSection) - || (out_point.z() <= z_min && fan_section != 0)) - { - // another section is crossed, no escape from bounding - distance = dist_from_p_to_out; - // select direction - if(v.z() < 0.) current_section = FanSection[fan_section - 1]; - else current_section = FanSection[fan_section + 1]; - A = out_point + v * current_section->DistanceToOut(out_point, v); - G4double side_a = Calculator->DistanceToTheNeutralFibre(A); - if(dist_p * side_a < 0.){ - distance += in_iteration_process(out_point, dist_out, A); - } else { - G4double dd = search_for_nearest_point(out_point, A); - if(( (fan_section == 0 && A.z() >= FanSectionLimits[2]) - || (fan_section == MaxFanSection && A.z() <= - FanSectionLimits[MaxFanSection - 1])) - && !(dd < kInfinity)) - { - if(fan_section == 0) B = A + v * FanSection[2]->DistanceToOut(A, v); - else B = A + v * FanSection[fan_section - 2]->DistanceToOut(A, v); - dd += search_for_nearest_point(A, B); - } - distance += dd; - } - } - return distance; -} -#ifdef DEBUG_LARWHEELSOLID + static G4ThreeVector q; +#ifdef LARWHEELSOLID_USE_FANBOUND + q = p + v * FanBound->DistanceToOut(p, v); +#else + find_exit_point(p, v, q); +#endif -G4double LArWheelSolid::in_chord_method( - const G4ThreeVector &p0, const G4ThreeVector &p1, - const G4ThreeVector &v) const -{ -#if 0 - if(Verbose) std::cout << "icm: " << p0 << " " << v << std::endl; - const G4double d0 = Calculator->DistanceToTheNeutralFibre(p0); - const G4double ht = d0 < 0.? FanHalfThickness: (-FanHalfThickness); - G4double f = d0 + ht; - if(Verbose) std::cout << "at p0: " << f << std::endl; + G4int start = select_section(p.z()); + G4int stop = select_section(q.z()); + G4int step = -1; + if(stop > start) { step = 1; start ++; stop ++; } + LWSDBG(5, std::cout << "dti sections " << start << " " << stop + << " " << step << std::endl); + static G4ThreeVector p1; + for(G4int i = start; i != stop; i += step){ +// v.z() can't be 0, otherwise start == stop, so the exit point could be only +// at z border of the fan section + const G4double d1 = (Zsect[i] - p.z()) / v.z(); + const G4double x1 = p.x() + v.x() * d1, y1 = p.y() + v.y() * d1; + p1.set(x1, y1, Zsect[i]); + G4double dist_p1 = Calculator->DistanceToTheNeutralFibre(p1); + LWSDBG(5, std::cout << i << ">" << p << " " << dist_p << " " + << p1 << " " << dist_p1 << std::endl); + G4double dd = kInfinity; + if(dist_p * dist_p1 < 0.){// it certanly cross current half-wave + dd = in_iteration_process(p, dist_p, p1); + } + G4double d2 = search_for_nearest_point(p, dist_p, p1); + LWSDBG(6, std::cout << i << "> dd=" << dd << ", d2=" << d2 << ", distance=" << distance << std::endl); + if(d2 < kInfinity){ + return distance + d2; // this half-wave is intersected + } else if(dd < kInfinity){ + return distance + dd; + } + distance += d1; + p.set(x1, y1, Zsect[i]); + dist_p = dist_p1; + } - { - FILE *F = fopen("test_icm", "r"); - if(F == 0){ - F = fopen("test_icm", "w"); - for(double x = 0.; x < (p1 - p0).mag(); x += IterationPrecision){ - G4ThreeVector pi = p0 + v * x; - fprintf(F, "%e %e\n", x, Calcuator->DistanceToTheNeutralFibre(pi) + ht); - } - std::cout << "Vz = " << v.z() << std::endl; - fclose(F); - } else fclose(F); - } - G4double F = Calculator->DistanceToTheNeutralFibre(p1) + ht; - G4double T = (p1 - p0).mag(); - if(Verbose) std::cout << "at p1: T " << T << " " << " F " << F << std::endl; - G4double t = 0.; - for(size_t i = 0; i < IterationsLimit; ++ i){ - const G4double dt = f * (t - T) / (f - F); - if(Verbose) std::cout << "at " << i << ": t " << t << " dt " << dt << " f " << f << std::endl; - if(fabs(dt) < IterationPrecision){ - return t; - } - if(std::signbit(f) != std::signbit - t -= dt; - f = Calculator->DistanceToTheNeutralFibre(p0 + v * t) + ht; - } -#endif - return 0.; - } -#endif + G4double dist_q = Calculator->DistanceToTheNeutralFibre(q); + LWSDBG(5, std::cout << "dti exit point: " << MSG_VECTOR(q) << " " + << dist_q << std::endl); + G4double dd = kInfinity; + if(dist_p * dist_q < 0.){// it certanly cross current half-wave + dd = in_iteration_process(p, dist_p, q); + } + G4double d2 = search_for_nearest_point(p, dist_p, q); + if(d2 < kInfinity){ + return distance + d2; // this half-wave is intersected + } else if(dd < kInfinity){ + return distance + dd; + } + return kInfinity; +} diff --git a/Simulation/G4Utilities/Geo2G4/src/LArWheelSolidDisToOut.cxx b/Simulation/G4Utilities/Geo2G4/src/LArWheelSolidDisToOut.cxx old mode 100755 new mode 100644 index 1dd42f8fe1844d2c9fc80855b7a1b19dda6b6527..a19893c4b52872d8bf8617f5902563130a40ab06 --- a/Simulation/G4Utilities/Geo2G4/src/LArWheelSolidDisToOut.cxx +++ b/Simulation/G4Utilities/Geo2G4/src/LArWheelSolidDisToOut.cxx @@ -2,273 +2,302 @@ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration */ -// November 2001, A. Soukharev -// LArWheelSolidDisToOut.cc (see also LArWheelSolid*.cc) -// DistanceToOut functions and stuff for them - #include <cassert> -#include "AthenaBaseComps/AthMsgStreamMacros.h" -#include "G4Polycone.hh" -#include "CLHEP/Units/PhysicalConstants.h" - #include "GeoSpecialShapes/LArWheelCalculator.h" #include "LArWheelSolid.h" - -#include<stdio.h> +#include "LArFanSection.h" G4double LArWheelSolid::DistanceToOut(const G4ThreeVector &inputP) const { - static G4ThreeVector p; - if(BoundingPolycone->Inside(inputP) == kOutside){ -#ifdef DEBUG_LARWHEELSOLID - ATH_MSG_VERBOSE ( "DistanceToOut(p):" - << " point " << MSG_VECTOR(inputP) - << " is outside of BoundingPolycone."); -#endif // DEBUG_LARWHEELSOLID - return 0.; - } - G4double d0 = BoundingPolycone->DistanceToOut(inputP); - p = inputP; - G4double d = FanHalfThickness - fabs(Calculator->DistanceToTheNearestFan(p)); - if(d < Tolerance) return 0.; - if(d > d0) return d0; - else return d; + LWSDBG(1, std::cout << TypeStr() << " DisToOut" << MSG_VECTOR(inputP) << std::endl); + //static G4ThreeVector p; + if(BoundingShape->Inside(inputP) != kInside){ + LWSDBG(2, std::cout << "DistanceToOut(p):" + << " point " << MSG_VECTOR(inputP) + << " is not inside of the BoundingShape." + << std::endl); + return 0.; + } + G4ThreeVector p( inputP ); + const G4double d = FanHalfThickness - fabs(Calculator->DistanceToTheNearestFan(p)); + if(d < Tolerance){ + LWSDBG(2, std::cout << "already not inside " << MSG_VECTOR(p) << std::endl); + return 0.; + } + const G4double d0 = BoundingShape->DistanceToOut(inputP); + LWSDBG(2, std::cout << "dto " << d << " " << d0 << std::endl); + if(d > d0) return d0; + else return d; } G4double LArWheelSolid::DistanceToOut(const G4ThreeVector &inputP, const G4ThreeVector &inputV, const G4bool calcNorm, G4bool* validNorm, - G4ThreeVector* ) const + G4ThreeVector* sn) const { - if(calcNorm && validNorm != 0) *validNorm = false; - if(BoundingPolycone->Inside(inputP) == kOutside){ -#ifdef DEBUG_LARWHEELSOLID - ATH_MSG_VERBOSE ( "DistanceToOut(p, v):" - << " point " << MSG_VECTOR(inputP) - << " is outside of BoundingPolycone." - ); -#endif // DEBUG_LARWHEELSOLID - return 0.; - } - static G4ThreeVector p, v; - p = inputP; - if(fabs(Calculator->DistanceToTheNearestFan(p)) - FanHalfThickness < Tolerance){ - v = inputV; - v.rotateZ(p.phi() - inputP.phi()); -#ifndef DEBUG_LARWHEELSOLID - return distance_to_out(p, v); -#else - G4double old = distance_to_out(p, v); - G4double dd1 = distance_to_out_ref(p, v); - if(fabs(dd1 - old) > IterationPrecision){ - static int cnt = 0; - std::cout << "DTO " << p << "(" << p.perp() << ") " << v - << " dto: " << old << " ref: " << dd1 << std::endl; - std::cout << "DTO MISMATCH " << (old - dd1) << " " - << LArWheelSolidTypeString(Type) << std::endl; - Verbose = true; - distance_to_out(p, v); - std::cout << "-----" << std::endl; - distance_to_out_ref(p, v); - Verbose = false; - //if(cnt ++ > 10) exit(0); - if(cnt ++ < 10){ - FILE *F = fopen("test_dto", "a"); - fwrite(&Type, sizeof(Type), 1, F); - fwrite(&p, sizeof(p), 1, F); - fwrite(&v, sizeof(v), 1, F); - fclose(F); - } - } - size_t i = 1; - for(G4double t = 0.; t < old; t += IterationPrecision * 10.){ - G4ThreeVector b = p + v * t; - if(fabs(Calculator->DistanceToTheNeutralFibre(b)) - FanHalfThickness > Tolerance){ - if(fabs(t - old) <= IterationPrecision) break; - std::cout << "@@ " << Type << " " << p << " " << v << " "; - std::cout << i << " " << old << " " << t << std::endl; - break; - } - i ++; - } - return old; + LWSDBG(1, std::cout << TypeStr() << " DisToOut" << MSG_VECTOR(inputP) + << MSG_VECTOR(inputV) << std::endl); + + const EInside inside_BS = BoundingShape->Inside(inputP); + if(inside_BS == kOutside){ + LWSDBG(2, std::cout << "DistanceToOut(p):" + << " point " << MSG_VECTOR(inputP) + << " is outside of BoundingShape." << std::endl); + if(calcNorm) *validNorm = false; + return 0.; + } + + // If here inside or on surface of BS + G4ThreeVector p(inputP); + + const G4double adtnf_p = fabs(Calculator->DistanceToTheNearestFan(p)); + if(adtnf_p >= FHTplusT) { + LWSDBG(2, std::cout << "DistanceToOut(p, v): point " + << MSG_VECTOR(inputP) + << " is outside of solid." << std::endl); + if(calcNorm) *validNorm = false; + return 0.; + } + + G4ThreeVector v(inputV); + const G4double phi0 = p.phi() - inputP.phi(); + v.rotateZ(phi0); + +#ifdef CHECK_DIRTONORM_ANGLE_ON_SURFACE + if(adtnf_p < FHTminusT) { + LWSDBG(5, std::cout << "inside fan point " << MSG_VECTOR(inputP) << ", FHTminusT=" << FHTminusT << std::endl); + } else { + LWSDBG(5, std::cout << "on fan surface adtnf_p=" << adtnf_p << ", FHTplusT=" << FHTplusT << ", FHTminusT=" << FHTminusT << std::endl); + + const G4ThreeVector d = Calculator->NearestPointOnNeutralFibre(p); + // check dot product between normal and v + if ( (p-d).cosTheta(v) > AngularTolerance ) { + // direction "v" definitely pointing outside + // return 0.0 + return 0.; + } + } #endif - } -#ifdef DEBUG_LARWHEELSOLID - ATH_MSG_VERBOSE ( "DistanceToOut(p, v):" - << " point " << MSG_VECTOR(inputP) - << " is outside of solid." ); -#endif // DEBUG_LARWHEELSOLID - return 0.; -} -// calculates the distance from point p along vector v to the nearest -// surface of vertical fan -#ifdef LArWheelSolidDTO_NEW -G4double LArWheelSolid::distance_to_out_ref +// former distance_to_out starts here + LWSDBG(4, std::cout << "dto: " << MSG_VECTOR(p) << " " + << MSG_VECTOR(v) << std::endl); + + G4ThreeVector q(p); +#ifdef LARWHEELSOLID_USE_BS_DTO + const G4double dto_bs = BoundingShape->DistanceToOut( + inputP, inputV, calcNorm, validNorm, sn + ); + q = p + v * dto_bs; + if(q.y() < Ymin){ + LWSDBG(5, std::cout << "dto exit point too low " << MSG_VECTOR(q) << std::endl); + const G4double dy = (Ymin - p.y()) / v.y(); + q.setX(p.x() + v.x() * dy); + q.setY(Ymin); + q.setZ(p.z() + v.z() * dy); + } #else -G4double LArWheelSolid::distance_to_out + FanBoundExit_t exit = find_exit_point(p, v, q); + LWSDBG(5, std::cout << "dto exit " << exit << std::endl); #endif -(const G4ThreeVector &p, const G4ThreeVector &v) const -{ - static G4ThreeVector out_section, out, out1, diff, C; + LWSDBG(5, std::cout << "dto exit point " << MSG_VECTOR(q) << std::endl); - G4int fan_section = select_fan_section(p.z()); - assert(FanSection[fan_section]->Inside(p) != kOutside); - G4double dist_from_p_to_out_section = - FanSection[fan_section]->DistanceToOut(p, v); - /* workaround for bug in G4Polycone */ - if(dist_from_p_to_out_section > 10 * CLHEP::m){ -#ifdef DEBUG_LARWHEELSOLID - ATH_MSG_WARNING ( "distance_to_out: workaround for a 'bug' in G4Polycone " - << " at " << MSG_VECTOR(p) << ": " - << dist_from_p_to_out_section ); -#endif - return 0; - } + G4double distance = 0.; + G4int start = select_section(p.z()); + G4int stop = select_section(q.z()); + G4int step = -1; + if(stop > start){ step = 1; start ++; stop ++; } + LWSDBG(5, std::cout << "dto sections " << start << " " << stop << " " << step << std::endl); - out_section = p + v * dist_from_p_to_out_section; - G4bool outside = (fabs(Calculator->DistanceToTheNeutralFibre(out_section)) - - FanHalfThickness > Tolerance); - G4double dist_from_p_to_out; - if(outside){ - // we are in the conditions for out_iteration_process, - dist_from_p_to_out = out_iteration_process(p, out_section); - out = p + v * dist_from_p_to_out; - // but we are not sure the answer being right (see description) - // so make fruther calculations - } else { - out = out_section; - dist_from_p_to_out = dist_from_p_to_out_section; - } + G4double tmp; + G4ThreeVector p1, C; - if(search_for_most_remoted_point(p, out, C)){ - return out_iteration_process(p, C); - } - // if initial exiting (from FanSection) point was outside fan, - // out_iteration_process had given correct result in the very beginning of the - // procedure. Let's remember and return it. - G4double distance = dist_from_p_to_out; - if(outside) return distance; - // otherwise move search into neighbours - G4double z_min = FanSectionLimits[fan_section] + Tolerance; - G4double z_max = FanSectionLimits[fan_section + 1] - Tolerance; -#ifdef DEBUG_LARWHEELSOLID - if(Verbose) std::cout << z_min << " " << z_max << " " << out << " " << fan_section << std::endl; + for(G4int i = start; i != stop; i += step){ + const G4double d1 = (Zsect[i] - p.z()) / v.z(); +// v.z() can't be 0, otherwise start == stop, so the exit point could be only +// at z border of the fan section + LWSDBG(5, std::cout << "at " << i << " dist to zsec = " << d1 << std::endl); + const G4double x1 = p.x() + v.x() * d1, y1 = p.y() + v.y() * d1; + p1.set(x1, y1, Zsect[i]); + const G4double dd = fabs(Calculator->DistanceToTheNeutralFibre(p1)); + if(dd > FHTplusT){ + tmp = out_iteration_process(p, p1); + //while(search_for_most_remoted_point(p, out_section, C)){ + if(search_for_most_remoted_point(p, p1, C)){ + tmp = out_iteration_process(p, C); + } + distance += tmp; +#ifndef LARWHEELSOLID_USE_BS_DTO + exit = NoCross; #endif - if((out.z() >= z_max && fan_section < MaxFanSection) - || (out.z() <= z_min && fan_section > 0)){ // another section - // select direction - G4int next_fan_section = (v.z() < 0.)? fan_section - 1: fan_section + 1; - G4double d = FanSection[next_fan_section]->DistanceToOut(out, v); - static G4ThreeVector out_next_fs; - out_next_fs = out + v * d; -#ifdef DEBUG_LARWHEELSOLID - if(Verbose) std::cout << "next fs " << out_next_fs << std::endl; + goto end_dto; + } + if(search_for_most_remoted_point(p, p1, C)){ + distance += out_iteration_process(p, C); +#ifndef LARWHEELSOLID_USE_BS_DTO + exit = NoCross; +#endif + goto end_dto; + } + distance += d1; + p.set(x1, y1, Zsect[i]); + } + + if(fabs(Calculator->DistanceToTheNeutralFibre(q)) > FHTplusT){ + LWSDBG(5, std::cout << "q=" << MSG_VECTOR(q) << " outside fan cur distance=" << distance << ", FHTplusT=" << FHTplusT << std::endl); + tmp = out_iteration_process(p, q); +#ifndef LARWHEELSOLID_USE_BS_DTO + exit = NoCross; +#endif + } else { + tmp = (q - p).mag(); + } + //while(search_for_most_remoted_point(out, out1, C)){ + if(search_for_most_remoted_point(p, q, C)){ + tmp = out_iteration_process(p, C); +#ifndef LARWHEELSOLID_USE_BS_DTO + exit = NoCross; #endif - if(fabs(Calculator->DistanceToTheNeutralFibre(out_next_fs)) - FanHalfThickness - < Tolerance) - { - if(search_for_most_remoted_point(out, out_next_fs, C)){ - distance += out_iteration_process(out, C); - } else distance += d; - } else { - G4double d1 = out_iteration_process(out, out_next_fs); - static G4ThreeVector out1; - out1 = out + v * d1; - if(search_for_most_remoted_point(out, out1, C)){ - distance += out_iteration_process(out, C); - } else distance += d1; - } - } + } + distance += tmp; +// former distance_to_out ends here + end_dto: + LWSDBG(5, std::cout << "At end_dto distance=" << distance << std::endl); +#ifdef LARWHEELSOLID_USE_BS_DTO + if(calcNorm && distance < dto_bs) *validNorm = false; +#else + if(calcNorm){ + LWSDBG(5, std::cout << "dto exit1 " << exit << std::endl); + switch(exit){ + case ExitAtBack: + sn->set(0., 0., 1.); + *validNorm = true; + break; + case ExitAtFront: + sn->set(0., 0., -1.); + *validNorm = true; + break; + case ExitAtOuter: + q.rotateZ(-phi0); + sn->set(q.x(), q.y(), 0.); + if(q.z() <= Zmid) sn->setZ(- q.perp() * m_fs->Amax); + *validNorm = true; + break; + default: + *validNorm = false; + break; + } + } +#endif + #ifdef DEBUG_LARWHEELSOLID - if(Verbose) std::cout << "dto result " << distance << std::endl; + if(Verbose > 2){ + if(Verbose > 3){ + std::cout << MSG_VECTOR(inputP) + << " " << MSG_VECTOR(inputV) << std::endl; + } + std::cout << "DTO: " << distance << " "; + if (validNorm) { + std::cout << *validNorm << " " << MSG_VECTOR((*sn)); + } else { + std::cout << "Norm is not valid"; + } + std::cout << std::endl; + + if(Verbose > 3){ + G4ThreeVector p2 = inputP + inputV * distance; + EInside i = Inside(p2); + std::cout << "DTO hit at " << MSG_VECTOR(p2) << ", " + << inside(i) << std::endl; + } + } +#ifdef LWS_HARD_TEST_DTO + if(test_dto(inputP, inputV, distance)){ + if(Verbose == 1){ + std::cout << TypeStr() << " DisToOut" << MSG_VECTOR(inputP) + << MSG_VECTOR(inputV) << std::endl; + } + } #endif - return distance; +#endif + return distance; } // This functions should be called in the case when we are sure that // point p is NOT OUTSIDE of vertical fan and point out is NOT INSIDE vertical fan -// returns distance from point p to fan surface +// returns distance from point p to fan surface, sets last +// parameter to the found point // may give wrong answer - see description G4double LArWheelSolid::out_iteration_process(const G4ThreeVector &p, - const G4ThreeVector &out) const + G4ThreeVector &B) const { -#ifdef DEBUG_LARWHEELSOLID - if(Verbose) std::cout << "oip: " << p << " " << out; -#endif - assert(fabs(Calculator->DistanceToTheNeutralFibre(p)) - FanHalfThickness < Tolerance); - assert(fabs(Calculator->DistanceToTheNeutralFibre(out)) - FanHalfThickness > -Tolerance); - static G4ThreeVector A, B, C, diff; - A = p; - B = out; - unsigned int niter = 0; - do{ - C = A + B; - C *= 0.5; - if(fabs(Calculator->DistanceToTheNeutralFibre(C)) - FanHalfThickness - < Tolerance) A = C; - else B = C; - niter ++; - diff = A - B; - } while(diff.mag2() > IterationPrecision2 && niter < IterationsLimit); - assert(fabs(Calculator->DistanceToTheNeutralFibre(B)) - FanHalfThickness > Tolerance); - assert(niter < IterationsLimit); - diff = p - B; -#ifdef DEBUG_LARWHEELSOLID - if(Verbose) std::cout << " -> " << B << " " << diff.mag() << std::endl; -#endif - return diff.mag(); + LWSDBG(6, std::cout << "oip: " << p << " " << B); + assert(fabs(Calculator->DistanceToTheNeutralFibre(p)) < FHTplusT); + assert(fabs(Calculator->DistanceToTheNeutralFibre(B)) > FHTminusT); + static G4ThreeVector A, C, diff; + A = p; + unsigned int niter = 0; + do { + C = A + B; + C *= 0.5; + if(fabs(Calculator->DistanceToTheNeutralFibre(C)) < FHTplusT){ + A = C; + } else { + B = C; + } + niter ++; + diff = A - B; + } while(diff.mag2() > IterationPrecision2 && niter < IterationsLimit); + assert(fabs(Calculator->DistanceToTheNeutralFibre(B)) > FHTplusT); + assert(niter < IterationsLimit); + diff = p - B; + LWSDBG(7, std::cout << " -> " << B << " " << diff.mag()); + LWSDBG(6, std::cout << std::endl); + return diff.mag(); } // returns true if the point being outside vert. fan is found, also set C to // that point in this case // returns false if the whole track looks like being inside vert. fan -G4bool LArWheelSolid::search_for_most_remoted_point(const G4ThreeVector &a, - const G4ThreeVector &b, - G4ThreeVector &C) const +G4bool LArWheelSolid::search_for_most_remoted_point( + const G4ThreeVector &a, const G4ThreeVector &b, G4ThreeVector &C +) const { -#ifdef DEBUG_LARWHEELSOLID - if(Verbose) std::cout << "sfmrp " << a << " " << b << std::endl; -#endif - static G4ThreeVector A, B, diff, l; - diff = b - a; - if(diff.mag2() <= IterationPrecision2) return false; - A = a; - B = b; - l = diff.unit() * IterationPrecision; + LWSDBG(6, std::cout << "sfmrp " << a << " " << b << std::endl); + static G4ThreeVector A, B, diff, l; + diff = b - a; + if(diff.mag2() <= IterationPrecision2) return false; + A = a; + B = b; + l = diff.unit() * IterationPrecision; // find the most remoted point on the line AB // and check if it is outside vertical fan // small vector along the segment AB - G4double d1, d2; - unsigned int niter = 0; + G4double d1, d2; + unsigned int niter = 0; // searching for maximum of (Calculator->DistanceToTheNeutralFibre)^2 along AB - do{ - C = A + B; - C *= 0.5; - d1 = Calculator->DistanceToTheNeutralFibre(C); - if(fabs(d1) - FanHalfThickness > Tolerance){ + do { + C = A + B; + C *= 0.5; + d1 = Calculator->DistanceToTheNeutralFibre(C); + if(fabs(d1) > FHTplusT){ // here out_iteration_process gives the answer -#ifdef DEBUG_LARWHEELSOLID - if(Verbose){ - std::cout << "sfmrp -> " << C << " " << fabs(d1) - << " " << (C - a).unit() << " " << (C - a).mag() - << std::endl; - } -#endif - return true; - } + LWSDBG(7, std::cout << "sfmrp -> " << C << " " << fabs(d1) + << " " << (C - a).unit() << " " + << (C - a).mag() << std::endl); + return true; + } // sign of derivative //d1 = Calculator->DistanceToTheNeutralFibre(C + l); - d2 = Calculator->DistanceToTheNeutralFibre(C - l); - if(d1 * d1 - d2 * d2 > 0.) A = C; - else B = C; - niter ++; - diff = A - B; - } while(diff.mag2() > IterationPrecision2 && niter < IterationsLimit); + d2 = Calculator->DistanceToTheNeutralFibre(C - l); + if(d1 * d1 - d2 * d2 > 0.) A = C; + else B = C; + niter ++; + diff = A - B; + } while(diff.mag2() > IterationPrecision2 && niter < IterationsLimit); // the line entirely lies inside fan - assert(niter < IterationsLimit); - return false; + assert(niter < IterationsLimit); + return false; } diff --git a/Simulation/G4Utilities/Geo2G4/src/LArWheelSolidInit.cxx b/Simulation/G4Utilities/Geo2G4/src/LArWheelSolidInit.cxx old mode 100755 new mode 100644 index 737b7f3685f87c84cf84dc0063e5db88690ddcff..ab9265e7e8841d604b4e9b0eeb6f7ab50c111c46 --- a/Simulation/G4Utilities/Geo2G4/src/LArWheelSolidInit.cxx +++ b/Simulation/G4Utilities/Geo2G4/src/LArWheelSolidInit.cxx @@ -2,45 +2,42 @@ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration */ -// LArWheelSolid -#include "LArWheelSolid.h" -#include "LArFanSection.h" -#include "GeoSpecialShapes/LArWheelCalculator.h" -#include "AthenaBaseComps/AthMsgStreamMacros.h" -#include "G4Polycone.hh" -#include "G4GeometryTolerance.hh" - #include <cassert> #include <stdexcept> +#include <iostream> #include "CLHEP/Units/PhysicalConstants.h" +#include "AthenaBaseComps/AthMsgStreamMacros.h" +#include "G4GeometryTolerance.hh" +#include "G4Polycone.hh" + +#include "GeoSpecialShapes/LArWheelCalculator.h" +#include "LArWheelSolid.h" +#include "LArFanSection.h" #ifdef DEBUG_LARWHEELSOLID -G4bool Verbose; +G4int LArWheelSolid::Verbose = 0; #endif - // these are internal technical constants, should not go in DB const unsigned int LArWheelSolid::IterationsLimit = 50; // That's enough even for 10e-15 IterationPrecision const G4double LArWheelSolid::Tolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance() / 2; +const G4double LArWheelSolid::AngularTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance() / 2; const G4double LArWheelSolid::IterationPrecision = 0.001*CLHEP::mm; const G4double LArWheelSolid::IterationPrecision2 = IterationPrecision * IterationPrecision; LArWheelSolid::LArWheelSolid(const G4String& name, LArWheelSolid_t type, G4int zside) - : G4VSolid(name),m_fs(0),m_msg("LArWheelSolid"),f_area(0),f_vol(0),f_area_on_pc(0),f_length(0),f_side_area(0) + : G4VSolid(name), Type(type), PhiPosition(CLHEP::halfpi), m_fs(0), + m_msg("LArWheelSolid"), + f_area(0), f_vol(0), f_area_on_pc(0), f_length(0), f_side_area(0) { - ATH_MSG_DEBUG ( "constructor started" ); - -#ifdef LArWheelSolidDTI_NEW - ATH_MSG_INFO ( "compiled with new DTI" ); -#endif -#ifdef LArWheelSolidDTO_NEW - ATH_MSG_INFO ( "compiled with new DTO" ); +#ifdef LARWHEELSOLID_USE_FANBOUND + ATH_MSG_INFO ( "compiled with G4 FanBound" ); +#else + ATH_MSG_INFO ( "compiled with private find_exit_point" ); #endif - Type = type; - PhiPosition = CLHEP::halfpi; LArWheelCalculator::LArWheelCalculator_t calc_type = LArWheelCalculator::LArWheelCalculator_t(0); switch(Type){ case InnerAbsorberWheel: @@ -84,13 +81,16 @@ LArWheelSolid::LArWheelSolid(const G4String& name, LArWheelSolid_t type, "Constructor: unknown LArWheelSolid_t"); } Calculator = new LArWheelCalculator(calc_type, zside); - G4String bp_name = name + "-BoundingPolycone"; + const G4String bs_name = name + "-Bounding"; #ifdef DEBUG_LARWHEELSOLID - Verbose = false; + const char *venv = getenv("LARWHEELSOLID_VERBOSE"); + if(venv) Verbose = atoi(venv); #endif // Initialize code that depends on wheel type: FanHalfThickness = Calculator->GetFanHalfThickness(); + FHTplusT = FanHalfThickness + Tolerance; + FHTminusT = FanHalfThickness - Tolerance; switch(Type){ case InnerAbsorberWheel: case InnerElectrodWheel: @@ -98,7 +98,7 @@ LArWheelSolid::LArWheelSolid(const G4String& name, LArWheelSolid_t type, case InnerElectrodModule: case InnerGlueWheel: case InnerLeadWheel: - inner_solid_init(bp_name); + inner_solid_init(bs_name); break; case OuterAbsorberWheel: case OuterElectrodWheel: @@ -106,17 +106,26 @@ LArWheelSolid::LArWheelSolid(const G4String& name, LArWheelSolid_t type, case OuterElectrodModule: case OuterGlueWheel: case OuterLeadWheel: - outer_solid_init(bp_name); + outer_solid_init(bs_name); break; default: G4Exception("LArWheelSolid", "UnknownSolidType", FatalException, "Constructor: unknown LArWheelSolid_t"); } - init_tests(); + Zsect_start_search = (Zsect.size() - 1) - 1; + init_tests(); test(); // activated by env. variable clean_tests(); + +#ifdef DEBUG_LARWHEELSOLID + m_fs->print(); + std::cout << "Limits: (" << Zsect.size() << ")" << std::endl; + for(size_t i = 0; i < Zsect.size(); ++ i){ + std::cout << i << " " << Zsect[i] << std::endl; + } +#endif ATH_MSG_DEBUG ( "solid of type " << LArWheelCalculator::LArWheelCalculatorTypeString(calc_type) << " initialized" ); @@ -124,302 +133,134 @@ LArWheelSolid::LArWheelSolid(const G4String& name, LArWheelSolid_t type, LArWheelSolid::~LArWheelSolid() { - delete [] FanSection; - if(m_fs) delete m_fs; + if(m_fs) delete m_fs; } // initialization of inner Absorber or Electrod wheels -void LArWheelSolid::inner_solid_init(const G4String &bp_name) +void LArWheelSolid::inner_solid_init(const G4String &bs_name) { - FanPhiAmplitude = 0.065; // internal technical constant, should not go in DB - set_phi_size(); - - G4double zPlane[2], rInner[3], rOuter[2]; - zPlane[0] = 0; - zPlane[1] = Calculator->GetWheelThickness(); - G4double wheel_thickness = zPlane[1] - zPlane[0]; - Calculator->GetWheelInnerRadius(rInner); // This can give us three radii - only want two... - Calculator->GetWheelOuterRadius(rOuter); - - Zmin = zPlane[0]; Zmax = zPlane[1]; - Rmin = rInner[0]; Rmax = rOuter[1]; - - BoundingPolycone = new G4Polycone(bp_name, MinPhi, MaxPhi - MinPhi, - 2, zPlane, rInner, rOuter); - - G4int number_of_half_waves = Calculator->GetNumberOfHalfWaves(); - FanSection = new G4Polycone* [number_of_half_waves + 2]; - MaxFanSection = number_of_half_waves + 1; - // We heed two additional fan sections to handle begin and end folds correctly - MaxFanSectionLimits = MaxFanSection + 1; - FanSectionLimits = new G4double[MaxFanSectionLimits + 1]; - - G4double sss = Calculator->GetStraightStartSection(); - G4double half_wave_length = Calculator->GetHalfWaveLength(); - - G4double Ain = (rInner[1] - rInner[0]) / wheel_thickness; - G4double Aout = (rOuter[1] - rOuter[0]) / wheel_thickness; - G4double Bin = rInner[0] - Ain * zPlane[0]; - G4double Bout = rOuter[0] - Aout * zPlane[0]; - - G4double phi_min = CLHEP::halfpi - FanPhiAmplitude - Calculator->GetFanStepOnPhi() * 2; - G4double phi_max = CLHEP::halfpi + FanPhiAmplitude + Calculator->GetFanStepOnPhi() * 2; - - G4int i; - for(i = 2; i < MaxFanSection; i ++){ - FanSectionLimits[i] = half_wave_length * (i - 1) + sss; - } - FanSectionLimits[0] = 0.; - FanSectionLimits[1] = sss + half_wave_length * 0.25; - FanSectionLimits[MaxFanSectionLimits] = wheel_thickness - FanSectionLimits[0]; - FanSectionLimits[MaxFanSectionLimits - 1] = wheel_thickness - FanSectionLimits[1]; - - m_fs = new LArFanSections(Ain, Aout, Bin, Bout, Rmax*cos(phi_min)); - - m_fs->z.push_back(FanSectionLimits[0]); - m_fs->rmin2.push_back(rInner[0]*rInner[0]); - m_fs->rmax2.push_back(rOuter[0]*rOuter[0]); - - char tmp_str[10]; - for(i = 0; i <= MaxFanSection; i ++){ - zPlane[0] = FanSectionLimits[i]; - zPlane[1] = FanSectionLimits[i + 1]; - rInner[0] = Bin + zPlane[0] * Ain; - rOuter[0] = Bout + zPlane[0] * Aout; - rInner[1] = Bin + zPlane[1] * Ain; - rOuter[1] = Bout + zPlane[1] * Aout; -#if LAR_FAN_SECTIONS_MULT > 1 - if(i != 0 && i != MaxFanSection){ - const G4double dzi = (FanSectionLimits[i + 1] - FanSectionLimits[i]) - / LAR_FAN_SECTIONS_MULT; - for(G4int di = 1; di < LAR_FAN_SECTIONS_MULT; ++ di){ - const G4double zi = FanSectionLimits[i] + dzi * di; - m_fs->z.push_back(zi); - const G4double ri = Bin + zi * Ain; - m_fs->rmin2.push_back(ri*ri); - const G4double ro = Bout + zi * Aout; - m_fs->rmax2.push_back(ro*ro); - } - } -#endif - m_fs->z.push_back(FanSectionLimits[i + 1]); - m_fs->rmin2.push_back(rInner[1]*rInner[1]); - m_fs->rmax2.push_back(rOuter[1]*rOuter[1]); - - sprintf(tmp_str, "%i", i); - G4String fs_name = bp_name + "-FanSection-" + tmp_str; -#ifdef DEBUG_LARWHEELSOLID - ATH_MSG_INFO ( fs_name << G4endl - << " z = (" << zPlane[0] << "," << zPlane[1] << ")" << G4endl - << " rIn = (" << rInner[0] << "," << rInner[1] << ")" << G4endl - << " rOut = (" << rOuter[0] << "," << rOuter[1] << ")" ); -#endif - FanSection[i] = new G4Polycone(fs_name, phi_min, phi_max - phi_min, - 2, zPlane, rInner, rOuter); -#ifdef DEBUG_LARWHEELSOLID - ATH_MSG_INFO ( " phi = (" << FanSection[i]->GetStartPhi() - << "," << FanSection[i]->GetEndPhi() << ")" ); + IsOuter = false; + FanPhiAmplitude = 0.065; // internal technical constant, should not go in DB + set_phi_size(); + + G4double zPlane[2], rInner[2], rOuter[2]; + zPlane[0] = 0.; + zPlane[1] = Calculator->GetWheelThickness(); + G4double wheel_thickness = zPlane[1] - zPlane[0]; + Calculator->GetWheelInnerRadius(rInner); + Calculator->GetWheelOuterRadius(rOuter); + const G4double phi_min = PhiPosition - FanPhiAmplitude + - Calculator->GetFanStepOnPhi() * 2; + + Zmin = zPlane[0]; Zmax = zPlane[1]; + Rmin = rInner[0]; Rmax = rOuter[1]; + Ymin = Rmin * 0.9; + Zmid = zPlane[1]; + + BoundingPolycone = new G4Polycone(bs_name + "Polycone", MinPhi, MaxPhi - MinPhi, + 2, zPlane, rInner, rOuter); + + BoundingShape = BoundingPolycone; +#ifdef LARWHEELSOLID_USE_FANBOUND + const G4double phi_size = (FanPhiAmplitude + Calculator->GetFanStepOnPhi() * 2) * 2; + FanBound = new G4Polycone(bs_name + "ofFan", phi_min, phi_size, + 2, zPlane, rInner, rOuter); #endif - } - m_fs->prepare(); -#ifdef DEBUG_LARWHEELSOLID - m_fs->print(); + ATH_MSG_INFO(BoundingShape->GetName() + " is the BoundingShape"); + + const G4double half_wave_length = Calculator->GetHalfWaveLength(); + const G4double sss = Calculator->GetStraightStartSection(); + Zsect.push_back(0.); + Zsect.push_back(sss + half_wave_length * 0.25); + const G4int num_fs = Calculator->GetNumberOfHalfWaves() + 1; + for(G4int i = 2; i < num_fs; i ++){ + const G4double zi = half_wave_length * (i - 1) + sss; +#if LARWHEELSOLID_ZSECT_MULT > 1 + for(G4int j = LARWHEELSOLID_ZSECT_MULT - 1; j > 0; -- j){ + Zsect.push_back(zi - half_wave_length * j / LARWHEELSOLID_ZSECT_MULT); + } #endif + Zsect.push_back(zi); + } + Zsect.push_back(wheel_thickness - Zsect[1]); + Zsect.push_back(wheel_thickness - Zsect[0]); + + m_fs = new LArFanSections( + rInner[0], rInner[1], rOuter[0], rOuter[1], + Rmax*cos(phi_min), Zsect.front(), Zsect.back() + ); } // initialization of outer Absorber or Electrod wheels -void LArWheelSolid::outer_solid_init(const G4String &bp_name) +void LArWheelSolid::outer_solid_init(const G4String &bs_name) { - FanPhiAmplitude = 0.02; // internal technical constant, should not go in DB - set_phi_size(); - - G4double zPlane[3], rInner[3], rOuter[3]; - zPlane[0] = 0; - zPlane[2] = Calculator->GetWheelThickness(); - G4double wheel_thickness = zPlane[2] - zPlane[0]; - zPlane[1] = Calculator->GetWheelInnerRadius(rInner); - Calculator->GetWheelOuterRadius(rOuter); - - Zmin = zPlane[0]; Zmax = zPlane[2]; - Rmin = rInner[0]; Rmax = rOuter[2]; - - BoundingPolycone = new G4Polycone(bp_name, MinPhi, MaxPhi - MinPhi, - 3, zPlane, rInner, rOuter); - - G4int number_of_half_waves = Calculator->GetNumberOfHalfWaves(); - FanSection = new G4Polycone* [number_of_half_waves + 2]; - MaxFanSection = number_of_half_waves + 1; - // We heed two additional fan sections to handle begin and end folds correctly - MaxFanSectionLimits = MaxFanSection + 1; - FanSectionLimits = new G4double[MaxFanSectionLimits + 1]; - - G4double sss = Calculator->GetStraightStartSection(); - G4double half_wave_length = Calculator->GetHalfWaveLength(); - - G4double phi_min = CLHEP::halfpi - FanPhiAmplitude - Calculator->GetFanStepOnPhi() * 2; - G4double phi_max = CLHEP::halfpi + FanPhiAmplitude + Calculator->GetFanStepOnPhi() * 2; - - G4int i; - for(i = 2; i < MaxFanSection; i ++){ - FanSectionLimits[i] = half_wave_length * (i - 1) + sss; - } - FanSectionLimits[0] = 0.; - FanSectionLimits[1] = sss + half_wave_length * 0.25; - FanSectionLimits[MaxFanSectionLimits] = wheel_thickness - FanSectionLimits[0]; - FanSectionLimits[MaxFanSectionLimits - 1] = wheel_thickness - FanSectionLimits[1]; - - // Note that as we go through the following loop, the [0] and [1] - // elements of the arrays change, but the [2] index remains - // constant except at the "bend." - - // Save the values at the "bend." - G4double zMid = zPlane[1]; - G4double rMidInner = rInner[1]; - G4double rMidOuter = rOuter[1]; - - char tmp_str[10]; - - // Calculate slopes and intercepts for beginning part - G4double Ain = (rInner[1] - rInner[0]) / (zPlane[1] - zPlane[0]); - G4double Aout = (rOuter[1] - rOuter[0]) / (zPlane[1] - zPlane[0]); - G4double Bin = rInner[0] - Ain * zPlane[0]; - G4double Bout = rOuter[0] - Aout * zPlane[0]; - - m_fs = new LArFanSections(Ain, Aout, Bin, Bout, Rmax*cos(phi_min)); - - m_fs->z.push_back(FanSectionLimits[0]); - m_fs->rmin2.push_back(rInner[0]*rInner[0]); - m_fs->rmax2.push_back(rOuter[0]*rOuter[0]); - - for(i = 0; i <= MaxFanSection; i ++){ - zPlane[0] = FanSectionLimits[i]; - zPlane[1] = FanSectionLimits[i + 1]; -#ifdef DEBUG_LARWHEELSOLID - ATH_MSG_INFO ( "zPlane[0]=" << zPlane[0] - << ", zMid=" << zMid << ", zPlane[1]=" - << zPlane[1] ); + IsOuter = true; + FanPhiAmplitude = 0.02; // internal technical constant, should not go in DB + set_phi_size(); + + G4double zPlane[3], rInner[3], rOuter[3]; + zPlane[0] = 0.; + zPlane[2] = Calculator->GetWheelThickness(); + G4double wheel_thickness = zPlane[2] - zPlane[0]; + zPlane[1] = Calculator->GetWheelInnerRadius(rInner); + Calculator->GetWheelOuterRadius(rOuter); + const G4double phi_min = PhiPosition - FanPhiAmplitude + - Calculator->GetFanStepOnPhi() * 2; + + Zmin = zPlane[0]; Zmax = zPlane[2]; + Rmin = rInner[0]; Rmax = rOuter[2]; + Ymin = Rmin * 0.9; + Zmid = zPlane[1]; + + BoundingPolycone = new G4Polycone(bs_name + "Polycone", MinPhi, MaxPhi - MinPhi, + 3, zPlane, rInner, rOuter); + BoundingShape = BoundingPolycone; +#ifdef LARWHEELSOLID_USE_FANBOUND + const G4double phi_size = (FanPhiAmplitude + Calculator->GetFanStepOnPhi() * 2) * 2; + FanBound = new G4Polycone(bs_name + "ofFan", phi_min, phi_size, + 3, zPlane, rInner, rOuter); #endif - if(zPlane[0] <= zMid && zMid < zPlane[1]){ - // We're creating a fan section with three z-planes (that is, - // this fan section includes the "bend" in the outer wheel). - - rInner[0] = Bin + zPlane[0] * Ain; - rOuter[0] = Bout + zPlane[0] * Aout; - // Recalculate the slopes and intercepts. - Ain = (rInner[2] - rMidInner) / (zPlane[2] - zMid); - Aout = (rOuter[2] - rMidOuter) / (zPlane[2] - zMid); - Bin = (zPlane[2] * rMidInner - rInner[2] * zMid) / - (zPlane[2] - zMid); - Bout = (zPlane[2] * rMidOuter - rOuter[2] * zMid) / - (zPlane[2] - zMid); - - zPlane[2] = zPlane[1]; - rInner[2] = Bin + zPlane[2] * Ain; - rOuter[2] = Bout + zPlane[2] * Aout; - zPlane[1] = zMid; - rInner[1] = rMidInner; - rOuter[1] = rMidOuter; - -#if LAR_FAN_SECTIONS_MULT > 1 - const G4double dzi = (FanSectionLimits[i + 1] - FanSectionLimits[i]) - / LAR_FAN_SECTIONS_MULT; - for(G4int di = 1; di < LAR_FAN_SECTIONS_MULT; ++ di){ - const G4double zi = FanSectionLimits[i] + dzi * di; - if(zi > zMid && m_fs->z.back() < zMid){ - m_fs->z.push_back(zMid); - m_fs->rmin2.push_back(rMidInner*rMidInner); - m_fs->rmax2.push_back(rMidOuter*rMidOuter); - } - m_fs->z.push_back(zi); - const G4double ri = Bin + zi * Ain; - m_fs->rmin2.push_back(ri*ri); - const G4double ro = Bout + zi * Aout; - m_fs->rmax2.push_back(ro*ro); - } -#else - m_fs->z.push_back(zMid); - m_fs->rmin2.push_back(rMidInner*rMidInner); - m_fs->rmax2.push_back(rMidOuter*rMidOuter); -#endif - m_fs->z.push_back(FanSectionLimits[i + 1]); - m_fs->rmin2.push_back(rInner[2]*rInner[2]); - m_fs->rmax2.push_back(rOuter[2]*rOuter[2]); - - sprintf(tmp_str, "%i", i); - G4String fs_name = bp_name + "-FanSection-" + tmp_str; -#ifdef DEBUG_LARWHEELSOLID - ATH_MSG_INFO ( fs_name << G4endl - << " z = (" << zPlane[0] << "," << zPlane[1] << "," << zPlane[2] << ")" << G4endl - << " rIn = (" << rInner[0] << "," << rInner[1] << "," << rInner[2] << ")" << G4endl - << " rOut = (" << rOuter[0] << "," << rOuter[1] << "," << rOuter[2] << ")" ); -#endif - FanSection[i] = new G4Polycone(fs_name, phi_min, phi_max - phi_min, - 3, zPlane, rInner, rOuter); -#ifdef DEBUG_LARWHEELSOLID - ATH_MSG_INFO ( " phi = (" << FanSection[i]->GetStartPhi() - << "," << FanSection[i]->GetEndPhi() << ")" ); -#endif - } else { - // We're creating a fan section with two z-planes (the more - // typical case). - rInner[0] = Bin + zPlane[0] * Ain; - rOuter[0] = Bout + zPlane[0] * Aout; - rInner[1] = Bin + zPlane[1] * Ain; - rOuter[1] = Bout + zPlane[1] * Aout; - -#if LAR_FAN_SECTIONS_MULT > 1 - if(i != 0 && i != MaxFanSection){ - const G4double dzi = (FanSectionLimits[i + 1] - FanSectionLimits[i]) - / LAR_FAN_SECTIONS_MULT; - for(G4int di = 1; di < LAR_FAN_SECTIONS_MULT; ++ di){ - const G4double zi = FanSectionLimits[i] + dzi * di; - m_fs->z.push_back(zi); - const G4double ri = Bin + zi * Ain; - m_fs->rmin2.push_back(ri*ri); - const G4double ro = Bout + zi * Aout; - m_fs->rmax2.push_back(ro*ro); - } - } -#endif - m_fs->z.push_back(FanSectionLimits[i + 1]); - m_fs->rmin2.push_back(rInner[1]*rInner[1]); - m_fs->rmax2.push_back(rOuter[1]*rOuter[1]); - - sprintf(tmp_str, "%i", i); - G4String fs_name = bp_name + "-FanSection-" + tmp_str; -#ifdef DEBUG_LARWHEELSOLID - ATH_MSG_INFO ( fs_name << G4endl - << " z = (" << zPlane[0] << "," << zPlane[1] << ")" << G4endl - << " rIn = (" << rInner[0] << "," << rInner[1] << ")" << G4endl - << " rOut = (" << rOuter[0] << "," << rOuter[1] << ")" ); -#endif - FanSection[i] = new G4Polycone(fs_name, phi_min, phi_max - phi_min, - 2, zPlane, rInner, rOuter); -#ifdef DEBUG_LARWHEELSOLID - ATH_MSG_INFO ( " phi = (" << FanSection[i]->GetStartPhi() - << "," << FanSection[i]->GetEndPhi() << ")" ); -#endif - } - } - m_fs->prepare(); -#ifdef DEBUG_LARWHEELSOLID - m_fs->print(); + ATH_MSG_INFO(BoundingShape->GetName() + "is the BoundingShape"); + + const G4double half_wave_length = Calculator->GetHalfWaveLength(); + const G4double sss = Calculator->GetStraightStartSection(); + Zsect.push_back(0.); + Zsect.push_back(sss + half_wave_length * 0.25); + const G4int num_fs = Calculator->GetNumberOfHalfWaves() + 1; + for(G4int i = 2; i < num_fs; i ++){ + const G4double zi = half_wave_length * (i - 1) + sss; +#if LARWHEELSOLID_ZSECT_MULT > 1 + for(G4int j = LARWHEELSOLID_ZSECT_MULT - 1; j > 0; -- j){ + G4double zj = zi - half_wave_length * j / LARWHEELSOLID_ZSECT_MULT; + Zsect.push_back(zj); + if(Zsect.back() <= Zmid && Zmid < zj){ + Zsect.push_back(Zmid); + } + } #endif + if(Zsect.back() <= Zmid && Zmid < zi){ + Zsect.push_back(Zmid); + } + Zsect.push_back(zi); + } + Zsect.push_back(wheel_thickness - Zsect[1]); + Zsect.push_back(wheel_thickness - Zsect[0]); + + m_fs = new LArFanSections( + rInner[0], rInner[1], rOuter[0], rOuter[1], + Rmax*cos(phi_min), Zsect.front(), Zmid + ); } // it should be called after FanPhiAmplitude has been set -// and before BoundingPolycone creation +// and before BoundingShape creation void LArWheelSolid::set_phi_size(void) { - if(Type == InnerAbsorberModule - || Type == InnerElectrodModule - || Type == OuterAbsorberModule - || Type == OuterElectrodModule) - { - MinPhi = PhiPosition - M_PI * (1. / 8.) - FanPhiAmplitude; - MaxPhi = PhiPosition + M_PI * (1. / 8.) + FanPhiAmplitude; - } else { - MinPhi = 0.; - MaxPhi = CLHEP::twopi; - } + if(Calculator->GetisModule()){ + MinPhi = PhiPosition - CLHEP::pi * (1. / 8.) - FanPhiAmplitude; + MaxPhi = PhiPosition + CLHEP::pi * (1. / 8.) + FanPhiAmplitude; + } else { + MinPhi = 0.; + MaxPhi = CLHEP::twopi; + } } diff --git a/Simulation/G4Utilities/Geo2G4/src/LArWheelSolidTests.cxx b/Simulation/G4Utilities/Geo2G4/src/LArWheelSolidTests.cxx index 159026029d12420c2d2c5377c07e5b5983fcb263..d0ac035ca54b8add0f2358a0812d900e09964fb6 100644 --- a/Simulation/G4Utilities/Geo2G4/src/LArWheelSolidTests.cxx +++ b/Simulation/G4Utilities/Geo2G4/src/LArWheelSolidTests.cxx @@ -2,23 +2,22 @@ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration */ -#include<iostream> -#include<stdexcept> -#include"boost/io/ios_state.hpp" -#include<map> - -#include"TRandom3.h" -#include"TF1.h" -#include"TNtupleD.h" -#include"TFile.h" - +#include <iostream> +#include <stdexcept> +#include "boost/io/ios_state.hpp" +#include <map> + +#include "TRandom3.h" +#include "TF1.h" +#include "TNtupleD.h" +#include "TFile.h" // For root version ifdef #include "TROOT.h" -#include"G4Polycone.hh" +#include "G4Polycone.hh" -#include"LArWheelSolid.h" #include "GeoSpecialShapes/LArWheelCalculator.h" +#include "LArWheelSolid.h" #include "CLHEP/Units/SystemOfUnits.h" #include "CLHEP/Units/PhysicalConstants.h" //#define LOCAL_DEBUG 1 @@ -78,7 +77,7 @@ G4ThreeVector LArWheelSolid::GetPointOnSurface(void) const const char *v1 = getenv("LARWHEELSOLID_TEST_MODE_LEVEL2"); if(v1) level2 = atof(v1); -#if DEBUG > 1 +#if LOCAL_DEBUG > 1 std::cout << "LWS::GPOS " << r << std::endl; #endif @@ -318,7 +317,7 @@ G4double LArWheelSolid::GetCubicVolume(void) f_vol->Integral(Rmin, Rmax, IntPrecision) #endif -#ifndef DEBUG +#ifndef LOCAL_DEBUG * Calculator->GetNumberOfFans() #endif ; @@ -391,7 +390,7 @@ G4double LArWheelSolid::GetSurfaceArea(void) // sagging ignored, effect should be negligible return result -#ifndef DEBUG +#ifndef LOCAL_DEBUG * Calculator->GetNumberOfFans() #endif ; @@ -399,37 +398,6 @@ G4double LArWheelSolid::GetSurfaceArea(void) void LArWheelSolid::test(void) { -#if 0 - FILE *FF = fopen("test_input", "r"); - if(FF){ - LArWheelSolid_t type; - fread(&type, sizeof(type), 1, FF); - if(type == Type){ - G4ThreeVector p0, v0; - fread(&p0, sizeof(p0), 1, FF); - fread(&v0, sizeof(v0), 1, FF); - Verbose = true; - std::cout << "AT TEST DTI" << p0 << " Rt = " << p0.perp() - << " phi = " << p0.phi() << " " << v0 << std::endl; - G4ThreeVector p1(p0), v1(v0), p(p0), v(v0); - G4double dd1 = distance_to_in_ref(p1, v1); - G4double dd = distance_to_in(p, v); - std::cout << "== DTI new " << dd1 << " old " << dd << " == " << std::endl; - EInside i1 = Inside(p0 + v0 * dd1); - EInside i = Inside(p0 + v0 * dd); - std::cout << "new " << inside(i1) << std::endl; - std::cout << "old " << inside(i) << std::endl; - - std::cout << std::endl << "AT TEST DTO" << p0 << " " << v0 << std::endl; - dd1 = distance_to_out_ref(p0, v0); - dd = distance_to_out(p0, v0); - std::cout << "== DTO new " << dd1 << " old " << dd << " == " << std::endl; - exit(0); - } - fclose(FF); - } -#endif - boost::io::ios_all_saver ias(std::cout); const char *on = getenv("LARWHEELSOLID_TEST"); if(on == 0) return; @@ -518,13 +486,29 @@ void LArWheelSolid::test(void) ias.restore(); } -void LArWheelSolid::clean_tests(void) -{ - if(f_area) delete f_area; - if(f_vol) delete f_vol; - if(f_area_on_pc) delete f_area_on_pc; - if(f_length) delete f_length; - if(f_side_area) delete f_side_area; +void LArWheelSolid::clean_tests(void) { + if(f_area) { + delete f_area; + f_area = 0; + } + if(f_vol) { + delete f_vol; + f_vol = 0; + } + + if(f_area_on_pc) { + delete f_area_on_pc; + f_area_on_pc = 0; + } + + if(f_length) { + delete f_length; + f_length = 0; + } + if(f_side_area) { + delete f_side_area; + f_side_area = 0; + } } G4double LArWheelSolid::get_area_at_r(G4double r) const @@ -631,11 +615,10 @@ void LArWheelSolid::init_tests(void) test_index = double(solid.size()); solid[test_index] = this; - std::cout << "LArWheelSolid::init_tests:" << std::endl; - std::cout << "Put " << this << " with index " << test_index - << " into list of solids" << std::endl; - std::cout << "Calculator: "<< Calculator - << " BP: " << BoundingPolycone << std::endl; +#ifdef DEBUG_LARWHEELSOLID + std::cout << "LArWheelSolid::init_tests: put " << this + << " with index " << test_index << std::endl; +#endif f_area = new TF1( (GetName() + "_f_area").c_str(), &LArWheelSolid_fcn_area, @@ -667,3 +650,147 @@ void LArWheelSolid::init_tests(void) ); f_side_area->FixParameter(0, test_index); } + +#ifdef DEBUG_LARWHEELSOLID +G4bool LArWheelSolid::test_dti( + const G4ThreeVector &inputP, const G4ThreeVector &inputV, + const G4double distance +) const +{ + if(distance > 10.*m){ + LWSDBG(1, std::cout << "DTI test skipped, distance > 10 m" + << std::endl); + return false; + } + static unsigned long counter = 0; + counter ++; + G4ThreeVector p, v; + if(BoundingShape->Inside(inputP) == kOutside){ + p = inputP + inputV * BoundingShape->DistanceToIn(inputP, inputV); + } else p = inputP; + const G4double phi0 = p.phi(); + const G4double d = Calculator->DistanceToTheNearestFan(p); + if(fabs(d) < FHTminusT){ + std::cout << "DTI test already inside" << MSG_VECTOR(p) << std::endl; + return false; + } + v = inputV; + v.rotateZ(p.phi() - phi0); + const G4double dd = IterationPrecision; + LWSDBG(1, std::cout << "Start DTI test, expect " + << long(distance / dd) << " iterations" + << std::endl); + + G4int V = Verbose; + Verbose = 0; + + const G4double d1 = distance - IterationPrecision; + static bool first = true; + for(G4double t = IterationPrecision; t < d1; t += dd){ + G4ThreeVector p1 = p + v * t; + if(fabs(Calculator->DistanceToTheNeutralFibre(p1)) < FHTminusT){ + std::cout << "DTI test at " << MSG_VECTOR(inputP) << " -> " + << MSG_VECTOR(inputV) << ", translated to " + << MSG_VECTOR(p) << " - > " << MSG_VECTOR(v) + << " in range of " + << distance << ": found nearer intersection at local point" + << MSG_VECTOR(p1) << ", distance " << d + << ", call " << counter + << std::endl; + Verbose = V; + + if(first){ + first = false; + FILE *F = fopen("dti_error.dat", "w"); + if(F){ + fprintf(F, "%10e %10e %10e\n", p.x(), p.y(), p.z()); + fprintf(F, "%10e %10e %10e\n", v.x(), v.y(), v.z()); + for(G4double e = IterationPrecision; e < d1; e += dd){ + p1 = p + v * e; + G4double f = fabs(Calculator->DistanceToTheNeutralFibre(p1)) - FanHalfThickness; + fprintf(F, "%10e %10e\n", e, f); + } + fclose(F); + } + } + + return true; + } + } + Verbose = V; + LWSDBG(1, std::cout << "DTI test at " << MSG_VECTOR(p) << " -> " + << MSG_VECTOR(v) << " in range of " + << distance << ": allright" << std::endl); + return false; +} + +G4bool LArWheelSolid::test_dto( + const G4ThreeVector &inputP, const G4ThreeVector &inputV, + const G4double distance +) const +{ + if(distance > 10.*m){ + LWSDBG(1, std::cout << "DTO test skipped, distance > 10 m" + << std::endl); + return false; + } + static unsigned long counter = 0; + counter ++; + G4ThreeVector p, v; + p = inputP; + const G4double phi0 = p.phi(); + const G4double d = Calculator->DistanceToTheNearestFan(p); + if(fabs(d) > FHTplusT){ + std::cout << "DTO test already outside" << MSG_VECTOR(p) << std::endl; + return false; + } + v = inputV; + v.rotateZ(p.phi() - phi0); + const G4double dd = IterationPrecision; + LWSDBG(1, std::cout << "Start DTO test, expect " + << long(distance / dd) << " iterations" + << std::endl); + + G4int V = Verbose; + Verbose = 0; + + const G4double d1 = distance - IterationPrecision; + static bool first = true; + for(G4double t = IterationPrecision; t < d1; t += dd){ + G4ThreeVector p1 = p + v * t; + if(fabs(Calculator->DistanceToTheNeutralFibre(p1)) > FHTplusT){ + std::cout << "DTO test at " << MSG_VECTOR(inputP) << " -> " + << MSG_VECTOR(inputV) << ", translated to " + << MSG_VECTOR(p) << " - > " << MSG_VECTOR(v) + << " in range of " + << distance << ": found nearer intersection at local point" + << MSG_VECTOR(p1) << ", distance " << d + << ", call " << counter + << std::endl; + Verbose = V; + + if(first){ + first = false; + FILE *F = fopen("dto_error.dat", "w"); + if(F){ + fprintf(F, "%10e %10e %10e\n", p.x(), p.y(), p.z()); + fprintf(F, "%10e %10e %10e\n", v.x(), v.y(), v.z()); + for(G4double e = IterationPrecision; e < d1; e += dd){ + p1 = p + v * e; + G4double f = fabs(Calculator->DistanceToTheNeutralFibre(p1)) - FanHalfThickness; + fprintf(F, "%10e %10e\n", e, f); + } + fclose(F); + } + } + + return true; + } + } + Verbose = V; + LWSDBG(1, std::cout << "DTO test at " << MSG_VECTOR(p) << " -> " + << MSG_VECTOR(v) << " in range of " + << distance << ": allright" << std::endl); + return false; +} +#endif diff --git a/Simulation/G4Utilities/Geo2G4/src/LArWheelSolid_type.h b/Simulation/G4Utilities/Geo2G4/src/LArWheelSolid_type.h new file mode 100644 index 0000000000000000000000000000000000000000..b30d4dc6c51bd19d9fc144b13241fee4711a6c29 --- /dev/null +++ b/Simulation/G4Utilities/Geo2G4/src/LArWheelSolid_type.h @@ -0,0 +1,28 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +// LArWheelSolid_type +// Author: D. A. Maximov +// typedefs for LArWheelSolid + +#ifndef __LArWheelSolid_type_HH__ +#define __LArWheelSolid_type_HH__ + +typedef enum { + InnerAbsorberWheel, + OuterAbsorberWheel, + InnerElectrodWheel, + OuterElectrodWheel, + InnerAbsorberModule, + OuterAbsorberModule, + InnerElectrodModule, + OuterElectrodModule, + InnerGlueWheel, + OuterGlueWheel, + InnerLeadWheel, + OuterLeadWheel +} LArWheelSolid_t; + + +#endif // __LArWheelSolid_type_HH__ diff --git a/Simulation/G4Utilities/Geo2G4/src/ParameterisedVolumeBuilder.cxx b/Simulation/G4Utilities/Geo2G4/src/ParameterisedVolumeBuilder.cxx deleted file mode 100755 index f2d5316ba8a4c65421df0b115c1a6e54d0a7386e..0000000000000000000000000000000000000000 --- a/Simulation/G4Utilities/Geo2G4/src/ParameterisedVolumeBuilder.cxx +++ /dev/null @@ -1,111 +0,0 @@ -/* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ - -#include "ParameterisedVolumeBuilder.h" -#include "Geo2G4LVFactory.h" -#include "Geo2G4STParameterisation.h" -#include "Geo2G4/LogicalVolume.h" -#include "GeoModelKernel/GeoAccessVolAndSTAction.h" -#include "GeoModelKernel/GeoLogVol.h" -#include "GeoModelKernel/GeoSerialTransformer.h" -#include "G4PVPlacement.hh" -#include "G4ReflectionFactory.hh" -#include "G4VPVParameterisation.hh" -#include "G4PVParameterised.hh" -#include "globals.hh" -#include <iostream> - - -// static ParameterisedVolumeBuilder b("Parameterised_Volume_Builder"); - -LogicalVolume* ParameterisedVolumeBuilder::Build(const PVConstLink theGeoPhysVolume, OpticalVolumesMap* /*optical_volumes*/) const -{ - PVConstLink theGeoPhysChild; - const GeoSerialTransformer* serialTransformerChild; - G4LogicalVolume* theG4LogChild; - unsigned int numChild; // number of child nodes (PV and ST) - bool descend; // flag to continue geo tree navigation - - // static SingleLogicalVolumeFactory LVFactory; - static Geo2G4LVFactory LVFactory; - - // const GeoLogVol* geoLog = theGeoPhysVolume->getLogVol(); - G4LogicalVolume* theG4LogVolume = LVFactory.Build(theGeoPhysVolume,descend); - - - if(!descend) return theG4LogVolume; - - // if (numChild==0) return theG4LogVolume; - // if (theG4LogVolume->GetNoDaughters()) return theG4LogVolume; - - // Loop over the children of the GeoVolume and place them - numChild = theGeoPhysVolume->getNChildVolAndST(); - - for(size_t ii=0; ii<numChild; ii++) - { - // Get the child Phys volume or SerialTransformer with ii index - GeoAccessVolAndSTAction action(ii); - theGeoPhysVolume->exec(&action); - - std::string nameChild = action.getName(); - int id=0; - - if(theGeoPhysChild=action.getVolume()) - { - // The ii node is a GeoPhysVol - // Build the child - theG4LogChild = Build(theGeoPhysChild); - if(!theG4LogChild) return 0; - - // Get its transform - const G4Transform3D theG4Position(action.getTransform()); - if (nameChild == "ANON") nameChild=theG4LogChild->GetName(); - - G4ReflectionFactory::Instance()->Place(theG4Position, - nameChild, - theG4LogChild, - theG4LogVolume, - false, - id); - } - else if((serialTransformerChild=action.getSerialTransformer())) - { - // The ii node is a GeoSerialTransformer - if(numChild > 1) - { - std::cerr << "\n\nParameterisedVolumeBuilder::Build ERROR building " << theGeoPhysVolume->getLogVol()->getName() - << ". In GEANT4 the parameterisation MUST be the mother's single daughter.\n"; - return 0; - } - else - { - theGeoPhysChild = serialTransformerChild->getVolume(); - - // Build the child - theG4LogChild = Build(theGeoPhysChild); - if(!theG4LogChild) return 0; - - if (nameChild == "ANON") nameChild=theG4LogChild->GetName(); - nameChild += "_Param"; - - Geo2G4STParameterisation* stParameterisation = new Geo2G4STParameterisation(serialTransformerChild->getFunction(), - serialTransformerChild->getNCopies()); - - G4VPhysicalVolume* pvParametrised __attribute__ ((unused)) = new G4PVParameterised(nameChild, - theG4LogChild, - theG4LogVolume, - kUndefined, - serialTransformerChild->getNCopies(), - stParameterisation); - } - } - else // node type is neither PV nor ST - { - std::cerr << "\n\nParameterisedVolumeBuilder::Build UNEXPECTED ERROR: action is returning 0.\n"; - return 0; - } - } - - return theG4LogVolume; -} diff --git a/Simulation/G4Utilities/Geo2G4/src/ParameterisedVolumeBuilder.h b/Simulation/G4Utilities/Geo2G4/src/ParameterisedVolumeBuilder.h deleted file mode 100755 index 37e32dd02541652aeb38ee91a05584438b8e3cc4..0000000000000000000000000000000000000000 --- a/Simulation/G4Utilities/Geo2G4/src/ParameterisedVolumeBuilder.h +++ /dev/null @@ -1,18 +0,0 @@ -/* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ - -#ifndef ParameterisedVolumeBuilder_H -#define ParameterisedVolumeBuilder_H - -#include "Geo2G4/VolumeBuilder.h" - -class ParameterisedVolumeBuilder: public VolumeBuilder -{ - public: - ParameterisedVolumeBuilder(std::string n):VolumeBuilder(n){} - - LogicalVolume* Build(PVConstLink pv, OpticalVolumesMap* optical_volumes = 0) const; -}; - -#endif diff --git a/Simulation/G4Utilities/Geo2G4/src/SingleLVCopyBuilder.cxx b/Simulation/G4Utilities/Geo2G4/src/SingleLVCopyBuilder.cxx deleted file mode 100755 index d946d3a4ae175a5ed2bc49ac77f94c1749efd205..0000000000000000000000000000000000000000 --- a/Simulation/G4Utilities/Geo2G4/src/SingleLVCopyBuilder.cxx +++ /dev/null @@ -1,63 +0,0 @@ -/* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ - -#include "SingleLVCopyBuilder.h" -#include "Geo2G4/LogicalVolume.h" -#include "SingleLogicalVolumeFactory.h" -#include "SimHelpers/ServiceAccessor.h" -#include "GeoModelKernel/GeoLogVol.h" -#include "AthenaBaseComps/AthMsgStreamMacros.h" -#include "G4PVPlacement.hh" -#include "G4ReflectionFactory.hh" -#include "globals.hh" -#include <iostream> - -LogicalVolume* SingleLVCopyBuilder::Build(const PVConstLink theGeoPhysVolume, OpticalVolumesMap* /*optical_volumes*/) const -{ - static SingleLogicalVolumeFactory LVFactory; - - const GeoLogVol * geoLog = theGeoPhysVolume->getLogVol(); - ATH_MSG_DEBUG ( "Start converting volume "<<geoLog->getName() ); - - G4LogicalVolume * theG4LogVolume = LVFactory.Build(geoLog); - - if (theGeoPhysVolume->getNChildVols()==0) return theG4LogVolume; - if (theG4LogVolume->GetNoDaughters() ) return theG4LogVolume; - // - // Loop over the children of the GeoVolume and place them - // - for(size_t ii = 0; ii<theGeoPhysVolume->getNChildVols(); ii++) { - std::string nameChild = theGeoPhysVolume->getNameOfChildVol(ii); - // - // Get the id from GeoModel - Query<int> Qint = theGeoPhysVolume->getIdOfChildVol(ii); - int id = 90999; - if(Qint.isValid() ) id = Qint; - // - // Get the child Phys volume ii - // - PVConstLink theGeoPhysChild = theGeoPhysVolume->getChildVol(ii); - // - // Build the child - // - G4LogicalVolume* theG4LogChild = Build(theGeoPhysChild); - - // Get its transform - const G4Transform3D theG4Position(theGeoPhysVolume->getXToChildVol(ii)); - - if (nameChild == "ANON") nameChild=theG4LogChild->GetName(); - // log<<MSG::VERBOSE<<"\t Positioning "<<theG4LogChild->GetName()<< - // " into "<<theG4LogVolume->GetName()<< " with name "<<nameChild - // << "and Id:" << id <<endreq; - G4ReflectionFactory::Instance()->Place(theG4Position, - nameChild, - theG4LogChild, - theG4LogVolume, - false, - id); - - } - - return theG4LogVolume; -} diff --git a/Simulation/G4Utilities/Geo2G4/src/SingleLVCopyBuilder.h b/Simulation/G4Utilities/Geo2G4/src/SingleLVCopyBuilder.h deleted file mode 100755 index 3f7de98beaa3d879018bc7e557f68f592fc9bd17..0000000000000000000000000000000000000000 --- a/Simulation/G4Utilities/Geo2G4/src/SingleLVCopyBuilder.h +++ /dev/null @@ -1,28 +0,0 @@ -/* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ - -#ifndef SingleLVCopyBuilder_H -#define SingleLVCopyBuilder_H - -#include "Geo2G4/VolumeBuilder.h" -#include "AthenaKernel/MsgStreamMember.h" -#include <string> - -class SingleLVCopyBuilder: public VolumeBuilder -{ - public: - SingleLVCopyBuilder(std::string n):VolumeBuilder(n),m_msg(n){} - - virtual LogicalVolume* Build(PVConstLink pv, OpticalVolumesMap* optical_volumes = 0) const; - /// Log a message using the Athena controlled logging system - MsgStream& msg( MSG::Level lvl ) const { return m_msg << lvl; } - /// Check whether the logging system is active at the provided verbosity level - bool msgLvl( MSG::Level lvl ) const { return m_msg.get().level() <= lvl; } - private: - /// Private message stream member - mutable Athena::MsgStreamMember m_msg; - -}; - -#endif diff --git a/Simulation/G4Utilities/Geo2G4/src/SingleLogicalVolumeFactory.cxx b/Simulation/G4Utilities/Geo2G4/src/SingleLogicalVolumeFactory.cxx old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/SingleLogicalVolumeFactory.h b/Simulation/G4Utilities/Geo2G4/src/SingleLogicalVolumeFactory.h old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/components/Geo2G4_entries.cxx b/Simulation/G4Utilities/Geo2G4/src/components/Geo2G4_entries.cxx old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/components/Geo2G4_load.cxx b/Simulation/G4Utilities/Geo2G4/src/components/Geo2G4_load.cxx old mode 100755 new mode 100644 diff --git a/Simulation/G4Utilities/Geo2G4/src/lcg_dict/selection.xml b/Simulation/G4Utilities/Geo2G4/src/lcg_dict/selection.xml new file mode 100644 index 0000000000000000000000000000000000000000..98e02df07a2c48bf657c52b7f5c1c1b3f31810af --- /dev/null +++ b/Simulation/G4Utilities/Geo2G4/src/lcg_dict/selection.xml @@ -0,0 +1,3 @@ +<lcgdict> + <class name="LArWheelSolidDDProxy" /> +</lcgdict>