From 43c3e3a7aef695ea36b1999fce1bf46c25b54ea2 Mon Sep 17 00:00:00 2001
From: Stewart Martin-Haugh <stewart.martin-haugh@cern.ch>
Date: Wed, 16 Sep 2020 08:56:34 +0000
Subject: [PATCH] Magnetic field MT migration: remove AtlasFieldSvc
(ATLASRECTS-5432)
---
.../TrigServices/src/HltEventLoopMgr.cxx | 10 -
.../InDetTrigPriVxFinder/CMakeLists.txt | 3 +-
.../InDetTrigPriVxFinder/TrigVxPrimary.h | 10 +-
.../InDetTrigPriVxFinder/TrigVxPrimaryAllTE.h | 10 +-
.../src/TrigVxPrimary.cxx | 26 +-
.../src/TrigVxPrimaryAllTE.cxx | 23 +-
.../MagFieldServices/AtlasFieldSvc.h | 261 ---
.../MagFieldServices/AtlasFieldSvcTLS.h | 48 -
.../python/MagFieldServicesConfig.py | 17 +-
.../python/MagFieldServicesConfigDb.py | 3 +-
.../python/MagFieldServicesSetup.py | 20 +-
.../MagFieldServices/python/SetupField.py | 6 +-
.../MagFieldServices/src/AtlasFieldSvc.cxx | 1563 -----------------
.../components/MagFieldServices_entries.cxx | 2 -
14 files changed, 30 insertions(+), 1972 deletions(-)
delete mode 100644 MagneticField/MagFieldServices/MagFieldServices/AtlasFieldSvc.h
delete mode 100644 MagneticField/MagFieldServices/MagFieldServices/AtlasFieldSvcTLS.h
delete mode 100644 MagneticField/MagFieldServices/src/AtlasFieldSvc.cxx
diff --git a/HLT/Trigger/TrigControl/TrigServices/src/HltEventLoopMgr.cxx b/HLT/Trigger/TrigControl/TrigServices/src/HltEventLoopMgr.cxx
index 1751b8d1533..c619cb671b5 100644
--- a/HLT/Trigger/TrigControl/TrigServices/src/HltEventLoopMgr.cxx
+++ b/HLT/Trigger/TrigControl/TrigServices/src/HltEventLoopMgr.cxx
@@ -834,16 +834,6 @@ StatusCode HltEventLoopMgr::updateMagField(const ptree& pt) const
auto tor_cur = pt.get<float>("Magnets.ToroidsCurrent.value");
auto sol_cur = pt.get<float>("Magnets.SolenoidCurrent.value");
- // Set currents on service (deprecated: ATLASRECTS-4687)
- IProperty* fieldSvc{nullptr};
- service("AtlasFieldSvc", fieldSvc, /*createIf=*/false).ignore();
- if ( fieldSvc ) {
- ATH_MSG_INFO("Setting field currents on AtlasFieldSvc");
- ATH_CHECK( Gaudi::Utils::setProperty(fieldSvc, "UseSoleCurrent", sol_cur) );
- ATH_CHECK( Gaudi::Utils::setProperty(fieldSvc, "UseToroCurrent", tor_cur) );
- }
- else ATH_MSG_WARNING("Cannot retrieve AtlasFieldSvc");
-
// Set current on conditions alg
const IAlgManager* algMgr = Gaudi::svcLocator()->as<IAlgManager>();
IAlgorithm* fieldAlg{nullptr};
diff --git a/InnerDetector/InDetTrigRecAlgs/InDetTrigPriVxFinder/CMakeLists.txt b/InnerDetector/InDetTrigRecAlgs/InDetTrigPriVxFinder/CMakeLists.txt
index 6d9444aab64..dbea24520c3 100644
--- a/InnerDetector/InDetTrigRecAlgs/InDetTrigPriVxFinder/CMakeLists.txt
+++ b/InnerDetector/InDetTrigRecAlgs/InDetTrigPriVxFinder/CMakeLists.txt
@@ -13,7 +13,6 @@ atlas_depends_on_subdirs( PUBLIC
Control/StoreGate
Event/EventPrimitives
InnerDetector/InDetRecTools/InDetRecToolInterfaces
- MagneticField/MagFieldInterfaces
InnerDetector/InDetConditions/InDetBeamSpotService
Tracking/TrkEvent/TrkEventPrimitives
Tracking/TrkEvent/TrkParameters
@@ -27,7 +26,7 @@ atlas_depends_on_subdirs( PUBLIC
atlas_add_component( InDetTrigPriVxFinder
src/*.cxx
src/components/*.cxx
- LINK_LIBRARIES GaudiKernel TrigInterfacesLib StoreGateLib SGtests EventPrimitives InDetRecToolInterfaces MagFieldInterfaces TrkEventPrimitives TrkParameters TrkParticleBase TrkTrack VxVertex TrigParticle xAODTracking InDetBeamSpotServiceLib )
+ LINK_LIBRARIES GaudiKernel TrigInterfacesLib StoreGateLib SGtests EventPrimitives InDetRecToolInterfaces TrkEventPrimitives TrkParameters TrkParticleBase TrkTrack VxVertex TrigParticle xAODTracking InDetBeamSpotServiceLib MagFieldConditions )
# Install files from the package:
atlas_install_headers( InDetTrigPriVxFinder )
diff --git a/InnerDetector/InDetTrigRecAlgs/InDetTrigPriVxFinder/InDetTrigPriVxFinder/TrigVxPrimary.h b/InnerDetector/InDetTrigRecAlgs/InDetTrigPriVxFinder/InDetTrigPriVxFinder/TrigVxPrimary.h
index c9e97d42b12..da19a2c0da1 100755
--- a/InnerDetector/InDetTrigRecAlgs/InDetTrigPriVxFinder/InDetTrigPriVxFinder/TrigVxPrimary.h
+++ b/InnerDetector/InDetTrigRecAlgs/InDetTrigPriVxFinder/InDetTrigPriVxFinder/TrigVxPrimary.h
@@ -27,6 +27,8 @@
#include <vector>
+#include "MagFieldConditions/AtlasFieldCacheCondObj.h"
+
/** Primary Vertex Finder. InDetTrigPriVxFinder uses the
InDetPrimaryVertexFinderTool in the package
InnerDetector/InDetRecTools/InDetPriVxFinderTool. It only gives the
@@ -34,12 +36,6 @@
VxContainer.
*/
-/* Forward declarations */
-
-namespace MagField {
- class IMagFieldSvc;
-}
-
class IBeamCondSvc;
namespace InDet
@@ -67,7 +63,7 @@ namespace InDet
bool m_createVtxTPLinks;
ToolHandle< IVertexFinder > m_VertexFinderTool;
- ServiceHandle< MagField::IMagFieldSvc> m_fieldSvc;
+ SG::ReadCondHandleKey<AtlasFieldCacheCondObj> m_fieldCondObjInputKey {this, "AtlasFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"};
ServiceHandle<IBeamCondSvc> m_BeamCondSvc;
};
diff --git a/InnerDetector/InDetTrigRecAlgs/InDetTrigPriVxFinder/InDetTrigPriVxFinder/TrigVxPrimaryAllTE.h b/InnerDetector/InDetTrigRecAlgs/InDetTrigPriVxFinder/InDetTrigPriVxFinder/TrigVxPrimaryAllTE.h
index 701d3618512..ec14dcae228 100755
--- a/InnerDetector/InDetTrigRecAlgs/InDetTrigPriVxFinder/InDetTrigPriVxFinder/TrigVxPrimaryAllTE.h
+++ b/InnerDetector/InDetTrigRecAlgs/InDetTrigPriVxFinder/InDetTrigPriVxFinder/TrigVxPrimaryAllTE.h
@@ -15,6 +15,8 @@
#include "TrigInterfaces/AllTEAlgo.h"
#include <vector>
+#include "MagFieldConditions/AtlasFieldCacheCondObj.h"
+
/** Primary Vertex Finder. InDetTrigPriVxFinder uses the
InDetPrimaryVertexFinderTool in the package
InnerDetector/InDetRecTools/InDetPriVxFinderTool. It only gives the
@@ -22,12 +24,6 @@
VxContainer.
*/
-/* Forward declarations */
-
-namespace MagField {
- class IMagFieldSvc;
-}
-
class IBeamCondSvc;
@@ -53,7 +49,7 @@ namespace InDet
bool m_runWithoutField;
ToolHandle< IVertexFinder > m_VertexFinderTool;
- ServiceHandle< MagField::IMagFieldSvc> m_fieldSvc;
+ SG::ReadCondHandleKey<AtlasFieldCacheCondObj> m_fieldCondObjInputKey {this, "AtlasFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"};
ServiceHandle<IBeamCondSvc> m_BeamCondSvc;
bool m_retrieve_tracks_from_SG;
diff --git a/InnerDetector/InDetTrigRecAlgs/InDetTrigPriVxFinder/src/TrigVxPrimary.cxx b/InnerDetector/InDetTrigRecAlgs/InDetTrigPriVxFinder/src/TrigVxPrimary.cxx
index 79fe7b91b97..664013376ca 100755
--- a/InnerDetector/InDetTrigRecAlgs/InDetTrigPriVxFinder/src/TrigVxPrimary.cxx
+++ b/InnerDetector/InDetTrigRecAlgs/InDetTrigPriVxFinder/src/TrigVxPrimary.cxx
@@ -28,7 +28,6 @@
#include "TrkEventPrimitives/ParamDefs.h"
-#include "MagFieldInterfaces/IMagFieldSvc.h"
#include "EventPrimitives/EventPrimitivesHelpers.h"
#include "InDetBeamSpotService/IBeamCondSvc.h"
@@ -43,7 +42,6 @@ namespace InDet
m_useTracksAsInput(false),
m_createVtxTPLinks(false),
m_VertexFinderTool("InDet::InDetPriVxFinderTool/InDetTrigPriVxFinderTool"),
- m_fieldSvc("AtlasFieldSvc", n),
m_BeamCondSvc("BeamCondSvc", n)
{
declareProperty("VertexFinderTool",m_VertexFinderTool);
@@ -73,26 +71,16 @@ namespace InDet
ATH_MSG_FATAL("Failed to retrieve tool " << m_VertexFinderTool);
return HLT::ErrorCode(HLT::Action::ABORT_JOB, HLT::Reason::BAD_JOB_SETUP);
}
- else{
- ATH_MSG_INFO("Retrieved tool " << m_VertexFinderTool);
- }
- if (m_fieldSvc.retrieve().isFailure()){
- ATH_MSG_FATAL("Failed to retrieve tool " << m_fieldSvc);
+ if (m_fieldCondObjInputKey.initialize().isFailure()){
+ ATH_MSG_FATAL("Failed to initialize " << m_fieldCondObjInputKey);
return HLT::ErrorCode(HLT::Action::ABORT_JOB, HLT::Reason::BAD_JOB_SETUP);
}
- else {
- ATH_MSG_INFO("Retrieved service " << m_fieldSvc);
- }
if (m_BeamCondSvc.retrieve().isFailure()){
ATH_MSG_FATAL("Failed to retrieve tool " << m_BeamCondSvc);
return HLT::ErrorCode(HLT::Action::ABORT_JOB, HLT::Reason::BAD_JOB_SETUP);
}
- else {
- ATH_MSG_INFO("Retrieved service " << m_fieldSvc);
-
- }
ATH_MSG_INFO("UseTracksAsInput: " << m_useTracksAsInput << " CreateVtxTrackParticleLinks: " << m_createVtxTPLinks);
return HLT::OK;
@@ -112,9 +100,15 @@ namespace InDet
bool runVtx(true);
- ATH_MSG_DEBUG(" In execHLTAlgorithm()");
+ EventContext ctx = Gaudi::Hive::currentContext();
+ MagField::AtlasFieldCache fieldCache;
+ SG::ReadCondHandle<AtlasFieldCacheCondObj> readHandle{m_fieldCondObjInputKey, ctx};
+ const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
+ fieldCondObj->getInitializedCache (fieldCache);
+
+
- if (!m_runWithoutField && !m_fieldSvc->solenoidOn()){
+ if (!m_runWithoutField && !fieldCache.solenoidOn()){
ATH_MSG_DEBUG("Solenoid Off and RunWithoutField=False. Algorithm not executed!");
runVtx = false;
}
diff --git a/InnerDetector/InDetTrigRecAlgs/InDetTrigPriVxFinder/src/TrigVxPrimaryAllTE.cxx b/InnerDetector/InDetTrigRecAlgs/InDetTrigPriVxFinder/src/TrigVxPrimaryAllTE.cxx
index f3815c2efb7..07a13e42256 100755
--- a/InnerDetector/InDetTrigRecAlgs/InDetTrigPriVxFinder/src/TrigVxPrimaryAllTE.cxx
+++ b/InnerDetector/InDetTrigRecAlgs/InDetTrigPriVxFinder/src/TrigVxPrimaryAllTE.cxx
@@ -14,7 +14,6 @@
#include "TrkEventPrimitives/ParamDefs.h"
#include "AthContainers/ConstDataVector.h"
-#include "MagFieldInterfaces/IMagFieldSvc.h"
#include "InDetBeamSpotService/IBeamCondSvc.h"
namespace InDet
@@ -24,7 +23,6 @@ namespace InDet
: HLT::AllTEAlgo(n, pSvcLoc),
m_runWithoutField(false),
m_VertexFinderTool("InDet::InDetPriVxFinderTool/InDetTrigPriVxFinderTool"),
- m_fieldSvc("AtlasFieldSvc", n),
m_BeamCondSvc("BeamCondSvc", n),
m_retrieve_tracks_from_SG(false),
m_track_collection_from_SG("Unconfigured_TrigVxPrimaryAllTE_For_SG_Access")
@@ -59,25 +57,16 @@ namespace InDet
msg() << MSG::FATAL << "Failed to retrieve tool " << m_VertexFinderTool << endmsg;
return HLT::ErrorCode(HLT::Action::ABORT_JOB, HLT::Reason::BAD_JOB_SETUP);
}
- else{
- msg() << MSG::INFO << "Retrieved tool " << m_VertexFinderTool << endmsg;
- }
- if (m_fieldSvc.retrieve().isFailure()){
- msg() << MSG::FATAL << "Failed to retrieve service " << m_fieldSvc << endmsg;
+ if (m_fieldCondObjInputKey.initialize().isFailure()){
+ ATH_MSG_FATAL("Failed to initialize " << m_fieldCondObjInputKey);
return HLT::ErrorCode(HLT::Action::ABORT_JOB, HLT::Reason::BAD_JOB_SETUP);
}
- else {
- msg() << MSG::INFO << "Retrieved tool " << m_fieldSvc << endmsg;
- }
if (m_BeamCondSvc.retrieve().isFailure()){
msg() << MSG::FATAL << "Failed to retrieve tool " << m_BeamCondSvc << endmsg;
return HLT::ErrorCode(HLT::Action::ABORT_JOB, HLT::Reason::BAD_JOB_SETUP);
}
- else {
- msg() << MSG::INFO << "Retrieved service " << m_fieldSvc << endmsg;
- }
return HLT::OK;
}
@@ -164,7 +153,13 @@ namespace InDet
runVtx = false;
}
- if (!m_runWithoutField && !m_fieldSvc->solenoidOn()){
+ EventContext ctx = Gaudi::Hive::currentContext();
+ MagField::AtlasFieldCache fieldCache;
+ SG::ReadCondHandle<AtlasFieldCacheCondObj> readHandle{m_fieldCondObjInputKey, ctx};
+ const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
+ fieldCondObj->getInitializedCache (fieldCache);
+
+ if (!m_runWithoutField && !fieldCache.solenoidOn()){
if(outputLevel <= MSG::DEBUG)
msg() << MSG::DEBUG << "Solenoid Off and RunWithoutField=False. Algorithm not executed!" << endmsg;
runVtx = false;
diff --git a/MagneticField/MagFieldServices/MagFieldServices/AtlasFieldSvc.h b/MagneticField/MagFieldServices/MagFieldServices/AtlasFieldSvc.h
deleted file mode 100644
index 15a5e27d213..00000000000
--- a/MagneticField/MagFieldServices/MagFieldServices/AtlasFieldSvc.h
+++ /dev/null
@@ -1,261 +0,0 @@
-/*
- Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
-*/
-
-///////////////////////////////////////////////////////////////////
-// AtlasFieldSvc.h, (c) ATLAS Detector software
-///////////////////////////////////////////////////////////////////
-
-#ifndef MAGFIELDSERVICES_ATLASFIELDSVC_H
-#define MAGFIELDSERVICES_ATLASFIELDSVC_H 1
-
-// FrameWork includes
-#include "AthenaBaseComps/AthService.h"
-#include "GaudiKernel/ToolHandle.h"
-#include "GaudiKernel/IIncidentListener.h"
-
-// MagField includes
-#include "MagFieldInterfaces/IMagFieldSvc.h"
-#include "MagFieldServices/AtlasFieldSvcTLS.h"
-#include "MagFieldElements/BFieldCache.h"
-#include "MagFieldElements/BFieldCacheZR.h"
-#include "MagFieldElements/BFieldCond.h"
-#include "MagFieldElements/BFieldZone.h"
-#include "MagFieldElements/BFieldMeshZR.h"
-
-// STL includes
-#include <vector>
-#include <iostream>
-
-#include "CxxUtils/checker_macros.h"
-
-// forward declarations
-class CondAttrListCollection;
-class BFieldZone;
-class TFile;
-class Incident;
-
-namespace MagField {
-
- /** @class AtlasFieldSvc
- @author Elmar.Ritsch -at- cern.ch
- */
-
- class ATLAS_NOT_THREAD_SAFE AtlasFieldSvc : public extends<AthService, IMagFieldSvc, IIncidentListener> {
- public:
-
- //** Constructor with parameters */
- AtlasFieldSvc( const std::string& name, ISvcLocator* pSvcLocator );
-
- /** Destructor */
- virtual ~AtlasFieldSvc();
-
- /** Athena algorithm's interface methods */
- virtual StatusCode initialize() override;
- virtual StatusCode finalize() override;
-
- /** Read **/
- virtual void handle(const Incident& runIncident) override;
-
- /** Call back for possible magnet current update **/
- StatusCode updateCurrent(IOVSVC_CALLBACK_ARGS);
-
- /** Call back for possible magnet filename update **/
- StatusCode updateMapFilenames(IOVSVC_CALLBACK_ARGS);
-
- /** get B field value at given position */
- /** xyz[3] is in mm, bxyz[3] is in kT */
- /** if deriv[9] is given, field derivatives are returned in kT/mm */
- virtual void getField(const double *xyz, double *bxyz, double *deriv = nullptr) const override final;
- virtual void getFieldZR(const double *xyz, double *bxyz, double *deriv = nullptr) const override final;
-
- private:
- /** Retrieve, initialize and return a thread-local storage object */
- inline struct AtlasFieldSvcTLS &getAtlasFieldSvcTLS() const;
-
- /* // Methods called to get field
- // pointer to actual function
- typedef void (AtlasFieldSvc::*FuncPtr)(const double *, double *, double *);
- FuncPtr getFieldActual;
- // standard calculation
- void getFieldStandard(const double *xyz, double *bxyz, double *deriv = 0);
- // manipulated field
- void getFieldManipulated(const double *xyz, double *bxyz, double *deriv = 0);
- */
-
- // Functions used by getField[ZR]
- // search for a "zone" to which the point (z,r,phi) belongs
- inline const BFieldZone* findZone( double z, double r, double phi ) const;
- // slow search is used during initialization to build the LUT
- BFieldZone* findZoneSlow( double z, double r, double phi );
- // fill the cache. return true if successful
- // return false if the position is outside the valid map volume
- inline bool fillFieldCache(double z, double r, double phi, AtlasFieldSvcTLS &tls) const;
- inline bool fillFieldCacheZR(double z, double r, AtlasFieldSvcTLS &tls) const;
-
- // set currents from when DCS is not read
- StatusCode importCurrents(AtlasFieldSvcTLS &tls);
- // initialize map
- StatusCode initializeMap(AtlasFieldSvcTLS &tls);
- // read the field map from an ASCII or ROOT file
- StatusCode readMap( const char* filename );
- StatusCode readMap( std::istream& input );
- StatusCode readMap( TFile* rootfile );
- // write the field map to a ROOT file
- void writeMap( TFile* rootfile ) const;
- // clear the field map
- void clearMap(AtlasFieldSvcTLS &tls);
-
- // utility functions used by readMap
- int read_packed_data( std::istream& input, std::vector<int>& data ) const;
- int read_packed_int( std::istream& input, int &n ) const;
- void buildLUT();
- void buildZR();
-
- // scale field according to the current
- void scaleField();
-
- /** approximate memory footprint in bytes */
- int memSize() const;
-
- /** Data Members **/
-
- // field map names
- std::string m_fullMapFilename; // all magnets on
- std::string m_soleMapFilename; // solenoid on / toroid off
- std::string m_toroMapFilename; // toroid on / solenoid off
- // current associated with the map
- double m_mapSoleCurrent; // solenoid current in A
- double m_mapToroCurrent; // toroid current in A
- // threshold below which currents are considered zero
- double m_soleMinCurrent; // minimum solenoid current to be considered ON
- double m_toroMinCurrent; // minimum toroid current to be considered ON
- // flag to use magnet current from DCS in COOL
- bool m_useDCS;
- // COOL folder name containing current information
- std::string m_coolCurrentsFolderName;
- // flag to read magnet map filenames from COOL
- bool m_useMapsFromCOOL;
- // COOL folder name containing field maps
- std::string m_coolMapsFolderName;
- // actual current if DCS is not in use
- double m_useSoleCurrent; // solenoid current in A
- double m_useToroCurrent; // toroid current in A
- // flag to skip current rescale and use map currents as they are
- bool m_lockMapCurrents;
-
- // handle for COOL field map filenames
- const DataHandle<CondAttrListCollection> m_mapHandle;
- // handle for COOL currents
- const DataHandle<CondAttrListCollection> m_currentHandle;
-
- // full 3d map (made of multiple zones)
- std::vector<BFieldZone> m_zone;
-
- // fast 2d map (made of one zone)
- BFieldMeshZR *m_meshZR;
-
- // data members used in zone-finding
- std::vector<double> m_edge[3]; // zone boundaries in z, r, phi
- std::vector<int> m_edgeLUT[3]; // look-up table for zone edges
- double m_invq[3]; // 1/stepsize in m_edgeLUT
- std::vector<const BFieldZone*> m_zoneLUT; // look-up table for zones
- // more data members to speed up zone-finding
- double m_zmin; // minimum z
- double m_zmax; // maximum z
- int m_nz; // number of z bins in zoneLUT
- double m_rmax; // maximum r
- int m_nr; // number of r bins in zoneLUT
- int m_nphi; // number of phi bins in zoneLUT
-
- /* // handle for field manipulator, if any
- bool m_doManipulation;
- ToolHandle<IMagFieldManipulator> m_manipulator; */
-
- };
-}
-
-// inline functions
-
-//
-// Initialize and return a thread-local storage object
-//
-struct MagField::AtlasFieldSvcTLS&
-MagField::AtlasFieldSvc::getAtlasFieldSvcTLS() const {
- static thread_local AtlasFieldSvcTLS tls = AtlasFieldSvcTLS();
- // return thread-local object
- return tls;
-}
-
-//
-// Search for the zone that contains a point (z, r, phi)
-// Fast version utilizing the LUT.
-//
-const BFieldZone*
-MagField::AtlasFieldSvc::findZone( double z, double r, double phi ) const
-{
- // make sure it's inside the largest zone
- // NB: the sign of the logic is chosen in order to return 0 on NaN inputs
- if ( z >= m_zmin && z <= m_zmax && r <= m_rmax ) {
- // find the edges of the zone
- // z
- const std::vector<double>& edgez(m_edge[0]);
- int iz = int(m_invq[0]*(z-m_zmin)); // index to LUT
- iz = m_edgeLUT[0][iz]; // tentative index from LUT
- if ( z > edgez[iz+1] ) iz++;
- // r
- const std::vector<double>& edger(m_edge[1]);
- int ir = int(m_invq[1]*r); // index to LUT - note minimum r is always 0
- ir = m_edgeLUT[1][ir]; // tentative index from LUT
- if ( r > edger[ir+1] ) ir++;
- // phi
- const std::vector<double>& edgephi(m_edge[2]);
- int iphi = int(m_invq[2]*(phi+M_PI)); // index to LUT - minimum phi is -pi
- iphi = m_edgeLUT[2][iphi]; // tentative index from LUT
- if ( phi > edgephi[iphi+1] ) iphi++;
- // use LUT to get the zone
- return m_zoneLUT[(iz*m_nr+ir)*m_nphi+iphi];
- } else {
- return nullptr;
- }
-}
-
-/** fill given magnetic field zone */
-bool
-MagField::AtlasFieldSvc::fillFieldCache(double z, double r, double phi, AtlasFieldSvcTLS &tls) const
-{
- // search for the zone
- const BFieldZone* zone = findZone( z, r, phi );
- if ( zone == nullptr ) {
- // outsize all zones
- return false;
- }
- // fill the cache
- zone->getCache( z, r, phi, tls.cache );
-
- // pointer to the conductors in the zone
- tls.cond = zone->condVector();
-
- // set a flag that the thread-local storage is initialized
- tls.isInitialized = true;
-
- return true;
-}
-
-/** fill Z-R cache for solenoid */
-bool
-MagField::AtlasFieldSvc::fillFieldCacheZR(double z, double r, AtlasFieldSvcTLS &tls) const
-{
- // is it inside the solenoid zone?
- if ( m_meshZR && m_meshZR->inside( z, r ) ) {
- // fill the cache
- m_meshZR->getCache( z, r, tls.cacheZR );
- } else {
- // outside solenoid
- return false;
- }
-
- return true;
-}
-
-#endif //> !MAGFIELDSERVICES_ATLASFIELDSVC_H
diff --git a/MagneticField/MagFieldServices/MagFieldServices/AtlasFieldSvcTLS.h b/MagneticField/MagFieldServices/MagFieldServices/AtlasFieldSvcTLS.h
deleted file mode 100644
index 1b9784384ed..00000000000
--- a/MagneticField/MagFieldServices/MagFieldServices/AtlasFieldSvcTLS.h
+++ /dev/null
@@ -1,48 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-/**
- * @author Elmar.Ritsch -at- cern.ch
- * @author Robert.Langenberg -at- cern.ch
- * @date March 2017
- * @brief Thread-local storage object used by MagField::AtlasFieldSvc
- */
-
-#ifndef MAGFIELDSERVICES_ATLASFIELDSVCTLS_H
-#define MAGFIELDSERVICES_ATLASFIELDSVCTLS_H 1
-
-// MagField includes
-#include "MagFieldElements/BFieldCond.h"
-#include "MagFieldElements/BFieldZone.h"
-#include "MagFieldElements/BFieldMeshZR.h"
-
-namespace MagField {
-
-/** @class AtlasFieldSvcTLS
- *
- * @brief Thread-local storage object used by MagField::AtlasFieldSvc
- *
- * @author Elmar.Ritsch -at- cern.ch
- * @author Robert.Langenberg -at- cern.ch
- */
-struct AtlasFieldSvcTLS {
-
- /// Constructor
- AtlasFieldSvcTLS() : isInitialized(false), cond(nullptr), cache(), cacheZR() { ; }
-
- /// Is the current AtlasFieldSvcTLS object properly initialized
- bool isInitialized;
-
- /// Pointer to the conductors in the current field zone (to compute Biot-Savart component)
- const std::vector<BFieldCond> *cond;
-
- /// Full 3d field
- BFieldCache cache;
- /// Fast 2d field
- BFieldCacheZR cacheZR;
-};
-
-} // namespace MagField
-
-#endif // MAGFIELDSERVICES_ATLASFIELDSVCTLS_H
diff --git a/MagneticField/MagFieldServices/python/MagFieldServicesConfig.py b/MagneticField/MagFieldServices/python/MagFieldServicesConfig.py
index 6183fb6b4db..1de8bd82eab 100644
--- a/MagneticField/MagFieldServices/python/MagFieldServicesConfig.py
+++ b/MagneticField/MagFieldServices/python/MagFieldServicesConfig.py
@@ -20,21 +20,6 @@ def MagneticFieldSvcCfg(flags, **kwargs):
result.merge(addFolders(flags, ['/EXT/DCS/MAGNETS/SENSORDATA'], detDb='DCS_OFL', className="CondAttrListCollection") )
- # AtlasFieldSvc - old one
- afsArgs = {
- "name": "AtlasFieldSvc",
- }
- if flags.Common.isOnline:
- afsArgs.update( UseDCS = False )
- afsArgs.update( UseSoleCurrent = 7730 )
- afsArgs.update( UseToroCurrent = 20400 )
- else:
- afsArgs.update( UseDCS = True )
- if 'UseDCS' in kwargs:
- afsArgs['UseDCS'] = kwargs['UseDCS']
- mag_field_svc = CompFactory.MagField.AtlasFieldSvc(**afsArgs)
- result.addService(mag_field_svc, primary=True)
-
# AtlasFieldMapCondAlg - for reading in map
afmArgs = {
"name": "AtlasFieldMapCondAlg",
@@ -91,7 +76,7 @@ if __name__=="__main__":
cfg=ComponentAccumulator()
acc = MagneticFieldSvcCfg(ConfigFlags)
- log.verbose(acc.getPrimary())
+ log.verbose(acc)
cfg.merge(acc)
diff --git a/MagneticField/MagFieldServices/python/MagFieldServicesConfigDb.py b/MagneticField/MagFieldServices/python/MagFieldServicesConfigDb.py
index 7ae3d142084..bf66b127b07 100644
--- a/MagneticField/MagFieldServices/python/MagFieldServicesConfigDb.py
+++ b/MagneticField/MagFieldServices/python/MagFieldServicesConfigDb.py
@@ -2,8 +2,7 @@
# database entries for https://twiki.cern.ch/twiki/bin/view/AtlasComputing/ConfiguredFactory#Factory_functions_vs_derived_cla
# Valerio Ippolito - Harvard University
-from AthenaCommon.CfgGetter import addService,addAlgorithm
+from AthenaCommon.CfgGetter import addAlgorithm
-addService('MagFieldServices.MagFieldServicesSetup.GetFieldSvc', 'AtlasFieldSvc')
addAlgorithm('MagFieldServices.MagFieldServicesSetup.GetFieldMapCondAlg', 'AtlasFieldMapCondAlg')
addAlgorithm('MagFieldServices.MagFieldServicesSetup.GetFieldCacheCondAlg', 'AtlasFieldCacheCondAlg')
diff --git a/MagneticField/MagFieldServices/python/MagFieldServicesSetup.py b/MagneticField/MagFieldServices/python/MagFieldServicesSetup.py
index 832602411d1..1f9275c7a64 100644
--- a/MagneticField/MagFieldServices/python/MagFieldServicesSetup.py
+++ b/MagneticField/MagFieldServices/python/MagFieldServicesSetup.py
@@ -1,6 +1,6 @@
# Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
-# JobOption fragment to set up the AtlasFieldSvc
+# JobOption fragment to set up the magnetic field services and algorithms
# Valerio Ippolito - Harvard University
# inspired by https://svnweb.cern.ch/trac/atlasoff/browser/MuonSpectrometer/MuonCnv/MuonCnvExample/trunk/python/MuonCalibConfig.py
@@ -9,22 +9,10 @@ from AthenaCommon.Logging import logging
logging.getLogger().info("Importing %s", __name__)
from AthenaCommon.AthenaCommonFlags import athenaCommonFlags
-from AthenaCommon.GlobalFlags import GlobalFlags
from AthenaCommon import CfgMgr
-MagField__AtlasFieldSvc=CfgMgr.MagField__AtlasFieldSvc
#--------------------------------------------------------------
-def AtlasFieldSvc(name="AtlasFieldSvc",**kwargs):
- if athenaCommonFlags.isOnline():
- kwargs.setdefault( "UseDCS", False )
- kwargs.setdefault( "UseSoleCurrent", 7730 )
- kwargs.setdefault( "UseToroCurrent", 20400 )
- else:
- kwargs.setdefault( "UseDCS", True )
-
- return CfgMgr.MagField__AtlasFieldSvc(name,**kwargs)
-
def AtlasFieldCacheCondAlg(name="AtlasFieldCacheCondAlg",**kwargs):
if athenaCommonFlags.isOnline():
kwargs.setdefault( "UseDCS", False )
@@ -50,12 +38,6 @@ def AtlasFieldMapCondAlg(name="AtlasFieldMapCondAlg",**kwargs):
def H8FieldSvc(name="H8FieldSvc",**kwargs):
return CfgMgr.MagField__H8FieldSvc(name,**kwargs)
-def GetFieldSvc(name="AtlasFieldSvc",**kwargs):
- if GlobalFlags.DetGeo == 'ctbh8':
- return H8FieldSvc(name, **kwargs)
- else:
- return AtlasFieldSvc(name, **kwargs)
-
def GetFieldCacheCondAlg(name="AtlasFieldCacheCondAlg",**kwargs):
return AtlasFieldCacheCondAlg(name, **kwargs)
diff --git a/MagneticField/MagFieldServices/python/SetupField.py b/MagneticField/MagFieldServices/python/SetupField.py
index 0806fe6e2ed..0b8bcf4d0dc 100755
--- a/MagneticField/MagFieldServices/python/SetupField.py
+++ b/MagneticField/MagFieldServices/python/SetupField.py
@@ -15,11 +15,7 @@ conddb.addFolderSplitMC('GLOBAL','/GLOBAL/BField/Maps <noover/>','/GLOBAL/BField
if not athenaCommonFlags.isOnline():
conddb.addFolder('DCS_OFL','/EXT/DCS/MAGNETS/SENSORDATA', className="CondAttrListCollection")
-# import the field service and conditions algs (service is 'to be removed')
-# (it is MagFieldServices.MagFieldServicesConfig which takes care of configuration)
-from AthenaCommon.CfgGetter import getService,getAlgorithm
-getService('AtlasFieldSvc')
-
+from AthenaCommon.CfgGetter import getAlgorithm
from AthenaCommon.AlgSequence import AthSequencer
condSequence = AthSequencer("AthCondSeq")
condSequence += getAlgorithm( "AtlasFieldMapCondAlg" )
diff --git a/MagneticField/MagFieldServices/src/AtlasFieldSvc.cxx b/MagneticField/MagFieldServices/src/AtlasFieldSvc.cxx
deleted file mode 100644
index 883c5901c23..00000000000
--- a/MagneticField/MagFieldServices/src/AtlasFieldSvc.cxx
+++ /dev/null
@@ -1,1563 +0,0 @@
-/*
- Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
-*/
-
-///////////////////////////////////////////////////////////////////
-// AtlasFieldSvc.cxx, (c) ATLAS Detector software
-///////////////////////////////////////////////////////////////////
-
-#include <iostream>
-
-// ISF_Services include
-#include "MagFieldServices/AtlasFieldSvc.h"
-
-// PathResolver
-#include "PathResolver/PathResolver.h"
-
-// StoreGate
-#include "StoreGate/StoreGateSvc.h"
-
-// Athena Pool
-#include "AthenaPoolUtilities/AthenaAttributeList.h"
-#include "AthenaPoolUtilities/CondAttrListCollection.h"
-
-// IncidentSvc
-#include "GaudiKernel/IIncidentSvc.h"
-
-// CLHEP
-#include "CLHEP/Units/SystemOfUnits.h"
-
-// ROOT
-#include "TFile.h"
-#include "TTree.h"
-
-/** Constructor **/
-MagField::AtlasFieldSvc::AtlasFieldSvc(const std::string& name,ISvcLocator* svc) :
- base_class(name,svc),
- m_fullMapFilename("MagneticFieldMaps/bfieldmap_7730_20400_14m.root"),
- m_soleMapFilename("MagneticFieldMaps/bfieldmap_7730_0_14m.root"),
- m_toroMapFilename("MagneticFieldMaps/bfieldmap_0_20400_14m.root"),
- m_mapSoleCurrent(7730.),
- m_mapToroCurrent(20400.),
- m_soleMinCurrent(1.0),
- m_toroMinCurrent(1.0),
- m_useDCS(false),
- m_coolCurrentsFolderName("/EXT/DCS/MAGNETS/SENSORDATA"),
- m_useMapsFromCOOL(true),
- m_coolMapsFolderName("/GLOBAL/BField/Maps"),
- m_useSoleCurrent(7730.),
- m_useToroCurrent(20400.),
- m_lockMapCurrents(false),
- m_mapHandle(),
- m_currentHandle(),
- m_zone(),
- m_meshZR(nullptr),
- m_edge(),
- m_edgeLUT(),
- m_invq(),
- m_zoneLUT(),
- m_zmin(0.),
- m_zmax(0.),
- m_nz(0),
- m_rmax(0.),
- m_nr(0),
- m_nphi(0)
- /* ,
- m_doManipulation(false),
- m_manipulator("undefined") */
-{
- declareProperty("FullMapFile", m_fullMapFilename, "File storing the full magnetic field map");
- declareProperty("SoleMapFile", m_soleMapFilename, "File storing the solenoid-only magnetic field map");
- declareProperty("ToroMapFile", m_toroMapFilename, "File storing the toroid-only magnetic field map");
- declareProperty("MapSoleCurrent", m_mapSoleCurrent, "Nominal solenoid current (A)");
- declareProperty("MapToroCurrent", m_mapToroCurrent, "Nominal toroid current (A)");
- declareProperty("SoleMinCurrent", m_soleMinCurrent, "Minimum solenoid current (A) for which solenoid is considered ON");
- declareProperty("ToroMinCurrent", m_toroMinCurrent, "Minimum toroid current (A) for which toroid is considered ON");
- declareProperty("UseDCS", m_useDCS, "Get magnet currents from DCS through COOL");
- declareProperty("COOLCurrentsFolderName", m_coolCurrentsFolderName, "Name of the COOL folder containing magnet currents");
- declareProperty("UseMapsFromCOOL", m_useMapsFromCOOL, "Get magnetic field map filenames from COOL");
- declareProperty("COOLMapsFolderName", m_coolMapsFolderName, "Name of the COOL folder containing field maps");
- declareProperty("UseSoleCurrent", m_useSoleCurrent, "Set actual solenoid current (A)");
- declareProperty("UseToroCurrent", m_useToroCurrent, "Set actual toroid current (A)");
- declareProperty("LockMapCurrents", m_lockMapCurrents, "Do not rescale currents (use the map values)");
- /* declareProperty("DoManipulation", m_doManipulation, "Apply field manipulation");
- declareProperty("ManipulatorTool", m_manipulator, "Tool handle for field manipulation"); */
-}
-
-MagField::AtlasFieldSvc::~AtlasFieldSvc()
-{
- delete m_meshZR;
-}
-
-/** framework methods */
-StatusCode MagField::AtlasFieldSvc::initialize( )
-{
- ATH_MSG_INFO( "initialize() ..." );
-
- // determine map location from COOL, if available
- if ( m_useMapsFromCOOL ) {
- // Register callback
- StoreGateSvc* detStore;
- if ( service( "DetectorStore", detStore ).isFailure() ) {
- ATH_MSG_FATAL( "Could not get detector store" );
- return StatusCode::FAILURE;
- }
- std::string folder( m_coolMapsFolderName );
- ATH_MSG_INFO("maps will be chosen reading COOL folder " << folder);
- if ( detStore->regFcn( &MagField::AtlasFieldSvc::updateMapFilenames, this,
- m_mapHandle, folder ).isFailure() ) {
- ATH_MSG_FATAL( "Could not book callback for " << folder );
- return StatusCode::FAILURE;
- }
- }
-
- // are we going to get the magnet currents from DCS?
- if ( m_useDCS ) {
- // Register callback
- StoreGateSvc* detStore;
- if ( service( "DetectorStore", detStore ).isFailure() ) {
- ATH_MSG_FATAL( "Could not get detector store" );
- return StatusCode::FAILURE;
- }
- std::string folder( m_coolCurrentsFolderName );
- ATH_MSG_INFO("magnet currents will be read from COOL folder " << folder);
- if ( detStore->regFcn( &MagField::AtlasFieldSvc::updateCurrent, this,
- m_currentHandle, folder ).isFailure() ) {
- ATH_MSG_FATAL( "Could not book callback for " << folder );
- return StatusCode::FAILURE;
- }
- ATH_MSG_INFO( "Booked callback for " << folder );
- // actual initalization has to wait for the fist callback
- } else {
- ATH_MSG_INFO( "Currents are set-up by jobOptions - delaying map initialization until BeginRun incident happens" );
-
- ServiceHandle<IIncidentSvc> incidentSvc("IncidentSvc", name());
- if (incidentSvc.retrieve().isFailure()) {
- ATH_MSG_FATAL( "Unable to retrieve the IncidentSvc" );
- return StatusCode::FAILURE;
- } else {
- incidentSvc->addListener( this, IncidentType::BeginRun );
- ATH_MSG_INFO( "Added listener to BeginRun incident" );
- }
- }
-
- // retrieve thread-local storage
- AtlasFieldSvcTLS &tls = getAtlasFieldSvcTLS();
-
- // clear the map for zero field
- clearMap(tls);
- setSolenoidCurrent(0.0);
- setToroidCurrent(0.0);
-
- /* // retrieve the manipulator tool
- if (m_doManipulation) {
- ATH_MSG_INFO( "field will be manipulated, retrieving tool" );
- if (m_manipulator.retrieve().isFailure()) {
- ATH_MSG_FATAL( "unable to retrieve manipulation tool" );
- } else {
- ATH_MSG_INFO( "manipulation tool retrieved" );
- getFieldActual = &MagField::AtlasFieldSvc::getFieldManipulated;
- }
- } else {
- ATH_MSG_INFO( "no manipulation set up" );
- getFieldActual = &MagField::AtlasFieldSvc::getFieldStandard;
- } */
-
- ATH_MSG_INFO( "initialize() successful" );
- return StatusCode::SUCCESS;
-}
-
-void MagField::AtlasFieldSvc::handle(const Incident& runIncident)
-{
- ATH_MSG_INFO( "handling incidents ..." );
- if ( !m_useDCS && runIncident.type() == IncidentType::BeginRun) {
- // get thread-local storage
- AtlasFieldSvcTLS &tls = getAtlasFieldSvcTLS();
-
- if ( importCurrents(tls).isFailure() ) {
- ATH_MSG_FATAL( "Failure in manual setting of currents" );
- } else {
- ATH_MSG_INFO( "BeginRun incident handled" );
- }
- }
- ATH_MSG_INFO( "incidents handled successfully" );
-}
-
-StatusCode MagField::AtlasFieldSvc::importCurrents(AtlasFieldSvcTLS &tls)
-{
- ATH_MSG_INFO( "importCurrents() ..." );
-
- // take the current values from JobOptions
- double solcur(m_useSoleCurrent);
- double torcur(m_useToroCurrent);
- if ( solcur < m_soleMinCurrent ) {
- solcur = 0.0;
- ATH_MSG_INFO( "Solenoid is off" );
- }
- if ( torcur < m_toroMinCurrent) {
- torcur = 0.0;
- ATH_MSG_INFO( "Toroids are off" );
- }
- setSolenoidCurrent(solcur);
- setToroidCurrent(torcur);
- // read the map file
- if ( initializeMap(tls).isFailure() ) {
- ATH_MSG_FATAL( "Failed to initialize field map" );
- return StatusCode::FAILURE;
- }
-
- ATH_MSG_INFO( "Currents imported and map initialized" );
- return StatusCode::SUCCESS;
-}
-
-/** callback for possible magnet current update **/
-StatusCode MagField::AtlasFieldSvc::updateCurrent(IOVSVC_CALLBACK_ARGS)
-{
- // get magnet currents from DCS
- double solcur(0.);
- double torcur(0.);
- bool gotsol(false);
- bool gottor(false);
-
- // due to inconsistencies between CONDBR2 and OFLP200/COMP200 (the former includes channel names
- // in the /EXT/DCS/MAGNETS/SENSORDATA folder, the latter don't), we try to read currents in
- // both ways
- bool hasChanNames(false);
-
- ATH_MSG_INFO( "Attempt 1 at reading currents from DCS (using channel name)" );
- for ( CondAttrListCollection::const_iterator itr = m_currentHandle->begin();
- itr != m_currentHandle->end(); ++itr ) {
-
- std::string name = m_currentHandle->chanName(itr->first);
- ATH_MSG_INFO( "Trying to read from DCS: [channel name, index, value] " << name << " , " << itr->first << " , " << itr->second["value"].data<float>() );
-
- if (name.compare("") != 0) {
- hasChanNames = true;
- }
-
- if ( name.compare("CentralSol_Current") == 0 ) {
- // channel 1 is solenoid current
- solcur = itr->second["value"].data<float>();
- gotsol = true;
- } else if ( name.compare("Toroids_Current") == 0 ) {
- // channel 3 is toroid current
- torcur = itr->second["value"].data<float>();
- gottor = true;
- }
- }
- if ( !hasChanNames ) {
- ATH_MSG_INFO( "Attempt 2 at reading currents from DCS (using channel index)" );
- // in no channel is named, try again using channel index instead
- for ( CondAttrListCollection::const_iterator itr = m_currentHandle->begin();
- itr != m_currentHandle->end(); ++itr ) {
-
- if ( itr->first == 1 ) {
- // channel 1 is solenoid current
- solcur = itr->second["value"].data<float>();
- gotsol = true;
- } else if ( itr->first == 3 ) {
- // channel 3 is toroid current
- torcur = itr->second["value"].data<float>();
- gottor = true;
- }
- }
- }
-
- if ( !gotsol || !gottor ) {
- if ( !gotsol ) ATH_MSG_ERROR( "Missing solenoid current in DCS information" );
- if ( !gottor ) ATH_MSG_ERROR( "Missing toroid current in DCS information" );
- return StatusCode::FAILURE;
- }
- ATH_MSG_INFO( "Currents read from DCS: solenoid " << solcur << " toroid " << torcur );
- // round to zero if close to zero
- if ( solcur < m_soleMinCurrent) {
- solcur = 0.0;
- ATH_MSG_INFO( "Solenoid is off" );
- }
- if ( torcur < m_toroMinCurrent) {
- torcur = 0.0;
- ATH_MSG_INFO( "Toroids are off" );
- }
- // did solenoid/toroid change status between on and off?
- bool solWasOn( solenoidOn() );
- bool torWasOn( toroidOn() );
- setSolenoidCurrent( solcur );
- setToroidCurrent( torcur );
- if ( solenoidOn() != solWasOn || toroidOn() != torWasOn ) {
- // get thread-local storage
- AtlasFieldSvcTLS &tls = getAtlasFieldSvcTLS();
-
- // map has changed. re-initialize the map
- if ( initializeMap(tls).isFailure() ) {
- ATH_MSG_ERROR( "Failed to re-initialize field map" );
- return StatusCode::FAILURE;
- }
- } else {
- // map is still valid. just scale the currents
- if (!m_lockMapCurrents)
- scaleField();
- else
- ATH_MSG_INFO( "Currents are NOT scaled - using map values sole=" << m_mapSoleCurrent << " toro=" << m_mapToroCurrent );
- }
-
- return StatusCode::SUCCESS;
-}
-
-/** callback for possible field map filenames update **/
-StatusCode MagField::AtlasFieldSvc::updateMapFilenames(IOVSVC_CALLBACK_ARGS)
-{
- ATH_MSG_INFO( "reading magnetic field map filenames from COOL" );
-
- ATH_MSG_DEBUG( "mapHandle isValid " << (int)m_mapHandle.isValid() );
- ATH_MSG_DEBUG( "mapHandle key " << m_mapHandle.key() );
- ATH_MSG_DEBUG( "mapHandle size " << m_mapHandle->size() );
-
-
- std::string fullMapFilename("");
- std::string soleMapFilename("");
- std::string toroMapFilename("");
-
- for (CondAttrListCollection::const_iterator itr = m_mapHandle->begin(); itr != m_mapHandle->end(); ++itr) {
- const coral::AttributeList &attr = itr->second;
- const std::string &mapType = attr["FieldType"].data<std::string>();
- const std::string &mapFile = attr["MapFileName"].data<std::string>();
- const float soleCur = attr["SolenoidCurrent"].data<float>();
- const float toroCur = attr["ToroidCurrent"].data<float>();
-
- ATH_MSG_INFO("found map of type " << mapType << " with soleCur=" << soleCur << " toroCur=" << toroCur << " (path " << mapFile << ")");
-
- // first 5 letters are reserved (like "file:")
- const std::string mapFile_decoded = mapFile.substr(5);
- if (mapType == "GlobalMap") {
- fullMapFilename = mapFile_decoded;
- m_mapSoleCurrent = soleCur;
- m_mapToroCurrent = toroCur;
- } else if (mapType == "SolenoidMap") {
- soleMapFilename = mapFile_decoded;
- } else if (mapType == "ToroidMap") {
- toroMapFilename = mapFile_decoded;
- }
- // note: the idea is that the folder contains exactly three maps
- // (if it contains more than 3 maps, then this logic doesn't work perfectly)
- // nominal currents are read from the global map
- }
-
- if (fullMapFilename.empty() || soleMapFilename.empty() || toroMapFilename.empty()) {
- ATH_MSG_ERROR("unexpected content in COOL field map folder");
- return StatusCode::FAILURE;
- }
-
- // check if maps really changed
- if (fullMapFilename != m_fullMapFilename || soleMapFilename != m_soleMapFilename || toroMapFilename != m_toroMapFilename) {
- ATH_MSG_INFO( "map set is new! reinitializing map");
- m_fullMapFilename = fullMapFilename;
- m_soleMapFilename = soleMapFilename;
- m_toroMapFilename = toroMapFilename;
-
- // retrieve the thread-local storage
- AtlasFieldSvcTLS &tls = getAtlasFieldSvcTLS();
-
- // trigger map reinitialization
- if ( initializeMap(tls).isFailure() ) {
- ATH_MSG_ERROR( "failed to re-initialize field map" );
- return StatusCode::FAILURE;
- }
- } else {
- ATH_MSG_INFO( "no need to update map set");
- }
-
- return StatusCode::SUCCESS;
-}
-
-//
-// read and initialize map
-//
-StatusCode MagField::AtlasFieldSvc::initializeMap(AtlasFieldSvcTLS &tls)
-{
- ATH_MSG_INFO( "Initializing the field map (solenoidCurrent=" << solenoidCurrent() << " toroidCurrent=" << toroidCurrent() << ")" );
- // empty the current map first
- clearMap(tls);
-
- // determine the map to load
- std::string mapFile("");
- if ( solenoidOn() && toroidOn() ) {
- mapFile = m_fullMapFilename;
- } else if ( solenoidOn() ) {
- mapFile = m_soleMapFilename;
- } else if ( toroidOn() ) {
- mapFile = m_toroMapFilename;
- } else {
- // all magnets OFF. no need to read map
- return StatusCode::SUCCESS;
- }
- // find the path to the map file
- std::string resolvedMapFile = PathResolver::find_file( mapFile.c_str(), "CALIBPATH" );
- if ( resolvedMapFile.empty() ) {
- ATH_MSG_ERROR( "Field map file " << mapFile << " not found" );
- return StatusCode::FAILURE;
- }
- // read the map file
- if ( readMap( resolvedMapFile.c_str() ).isFailure() ) {
- ATH_MSG_ERROR( "Something went wrong while trying to read the field map " << resolvedMapFile );
- return StatusCode::FAILURE;
- }
- ATH_MSG_INFO( "Initialized the field map from " << resolvedMapFile );
- // scale magnet current as needed
- if (!m_lockMapCurrents)
- scaleField();
- else
- ATH_MSG_INFO( "Currents are NOT scaled - using map values sole=" << m_mapSoleCurrent << " toro=" << m_mapToroCurrent );
-
- return StatusCode::SUCCESS;
-}
-
-void MagField::AtlasFieldSvc::scaleField()
-{
- BFieldZone *solezone(nullptr);
- //
- if ( solenoidOn() ) {
- solezone = findZoneSlow( 0.0, 0.0, 0.0 );
- if ( m_mapSoleCurrent > 0.0 &&
- std::abs( solenoidCurrent()/m_mapSoleCurrent - 1.0 ) > 0.001 ) {
- // scale the field in the solenoid zone
- double factor = solenoidCurrent()/m_mapSoleCurrent;
- solezone->scaleField( factor );
- // remake the fast map
- buildZR();
- ATH_MSG_INFO( "Scaled the solenoid field by a factor " << factor );
- }
- ATH_MSG_INFO( "Solenoid zone id " << solezone->id() );
- }
- //
- if ( toroidOn() ) {
- if ( m_mapToroCurrent > 0.0 &&
- std::abs( toroidCurrent()/m_mapToroCurrent - 1.0 ) > 0.001 ) {
- // scale the field in all zones except for the solenoid zone
- double factor = toroidCurrent()/m_mapToroCurrent;
- for ( unsigned i = 0; i < m_zone.size(); i++ ) {
- if ( &(m_zone[i]) != solezone ) {
- m_zone[i].scaleField( factor );
- }
- }
- ATH_MSG_INFO( "Scaled the toroid field by a factor " << factor );
- }
- // for ( unsigned i = 0; i < m_zone.size(); i++ ) {
- // ATH_MSG_INFO( "zone i,id " << i << ": " << m_zone[i].id() );
- // }
- }
-}
-
-/** framework methods */
-StatusCode MagField::AtlasFieldSvc::finalize()
-{
- //ATH_MSG_INFO( "finalize() ..." );
- //
- // finalization code would go here
- //
- ATH_MSG_INFO( "finalize() successful" );
- return StatusCode::SUCCESS;
-}
-
-/* void MagField::AtlasFieldSvc::getFieldStandard(const double *xyz, double *bxyz, double *deriv) */
-void MagField::AtlasFieldSvc::getField(const double *xyz, double *bxyz, double *deriv) const
-{
- const double &x(xyz[0]);
- const double &y(xyz[1]);
- const double &z(xyz[2]);
- double r = std::sqrt(x * x + y * y);
- double phi = std::atan2(y, x);
-
- // retrieve the thread-local storage
- AtlasFieldSvcTLS &tls = getAtlasFieldSvcTLS();
- BFieldCache &cache = tls.cache;
-
- // test if the TLS was initialized and the cache is valid
- if ( !tls.isInitialized || !cache.inside(z, r, phi) ) {
- // cache is invalid -> refresh cache
- if (!fillFieldCache(z, r, phi, tls)) {
- // caching failed -> outside the valid map volume
- // return default field (0.1 gauss)
- const double defaultB(0.1*CLHEP::gauss);
- bxyz[0] = bxyz[1] = bxyz[2] = defaultB;
- // return zero gradient if requested
- if ( deriv ) {
- for ( int i = 0; i < 9; i++ ) {
- deriv[i] = 0.;
- }
- }
- return;
- }
- }
-
- // do interpolation
- cache.getB(xyz, r, phi, bxyz, deriv);
-
- // add biot savart component
- if (tls.cond) {
- const size_t condSize = tls.cond->size();
- for (size_t i = 0; i < condSize; i++) {
- (*tls.cond)[i].addBiotSavart(1, xyz, bxyz, deriv); // added scale factor of 1 for compatibility with multi-threaded model
- }
- }
-}
-
-/*
-void MagField::AtlasFieldSvc::getFieldManipulated(const double *xyz, double *bxyz, double *deriv)
-{
- // this operation involves three steps:
- // - first we move the point at which the field is evaluated
- // ex1: solenoid translation by vector +a => xyz -= a
- // ex2: solenoid rotation R by angle +phi => rotate xyz by -phi (inverse rotation R_inv)
- // - then, we evaluate B in this new point
- // ex1: B(-a)
- // ex2: B(R_inv(xyz))
- // - then, we change the field
- // ex1: identity transformation
- // ex2: rotate the field properly
-
- // step 1
- double xyz_new[3];
- m_manipulator->modifyPosition(xyz, xyz_new);
-
- // step 2
- getFieldStandard(xyz_new, bxyz, deriv);
-
- // step 3
- m_manipulator->modifyField(bxyz, deriv);
-}
-
-void MagField::AtlasFieldSvc::getField(const double *xyz, double *bxyz, double *deriv) {
- (this->*this->getFieldActual)(xyz, bxyz, deriv);
-}
-*/
-
-void MagField::AtlasFieldSvc::getFieldZR(const double *xyz, double *bxyz, double *deriv) const
-{
- const double &x(xyz[0]);
- const double &y(xyz[1]);
- const double &z(xyz[2]);
- double r = sqrt(x * x + y * y);
-
- // get thread-local storage
- AtlasFieldSvcTLS &tls = getAtlasFieldSvcTLS();
- BFieldCacheZR &cacheZR = tls.cacheZR;
-
- // test if the TLS was initialized and the cache is valid
- if ( !tls.isInitialized || !cacheZR.inside(z, r) ) {
- // cache is invalid -> refresh cache
- if (!fillFieldCacheZR(z, r, tls)) {
- // caching failed -> outside the valid z-r map volume
- // call the full version of getField()
- getField(xyz, bxyz, deriv);
-
- return;
- }
- }
-
- // do interpolation
- cacheZR.getB(xyz, r, bxyz, deriv);
-
-}
-
-//
-// Clear the map.
-// Subsequent call should return zero magnetic field.
-//
-void MagField::AtlasFieldSvc::clearMap(AtlasFieldSvcTLS &tls)
-{
- tls.cache.invalidate();
- tls.cacheZR.invalidate();
-
- tls.cond = nullptr;
- // Next lines clear m_zone, m_edge[3], m_edgeLUT[3], and m_zoneLUT and deallocate their memory.
- std::vector<BFieldZone>().swap(m_zone);
- for ( int i = 0; i < 3; i++ ) {
- std::vector<double>().swap(m_edge[i]);
- std::vector<int>().swap(m_edgeLUT[i]);
- }
- std::vector<const BFieldZone*>().swap(m_zoneLUT);
- // Next lines ensure findZone() will fail
- m_zmin = 0.0;
- m_zmax = -1.0;
- m_rmax = -1.0;
- m_nz = m_nr = m_nphi = 0;
- delete m_meshZR;
- m_meshZR = nullptr;
-}
-
-//
-// Read the solenoid map from file.
-// The filename must end with ".root".
-//
-StatusCode MagField::AtlasFieldSvc::readMap( const char* filename )
-{
- if ( strstr(filename, ".root") == nullptr ) {
- ATH_MSG_ERROR("input file name '" << filename << "' does not end with .root");
- return StatusCode::FAILURE;
- }
- TFile* rootfile = new TFile(filename, "OLD");
- if ( ! rootfile ) {
- ATH_MSG_ERROR("failed to open " << filename);
- return StatusCode::FAILURE;
- }
- ATH_MSG_INFO("reading the map from " << filename);
- if ( readMap(rootfile).isFailure() ) {
- ATH_MSG_ERROR("something went wrong while trying to read the ROOT field map file");
- return StatusCode::FAILURE;
- }
-
- rootfile->Close();
- delete rootfile;
-
- return StatusCode::SUCCESS;
-}
-
-//
-// read an ASCII field map from istream
-// convert units m -> mm, and T -> kT
-//
-StatusCode MagField::AtlasFieldSvc::readMap( std::istream& input )
-{
- const std::string myname("readMap()");
- // first line contains version, date, time
- std::string word;
- int version;
- int date;
- int time;
- input >> word >> version;
- if ( word != "FORMAT-VERSION" ) {
- ATH_MSG_ERROR( myname << ": found '" << word << "' instead of 'FORMAT-VERION'");
- return StatusCode::FAILURE;
- }
- if ( version < 5 || version > 6) {
- ATH_MSG_ERROR( myname << ": version number is " << version << " instead of 5 or 6");
- return StatusCode::FAILURE;
- }
- input >> word >> date;
- if ( word != "DATE" ) {
- ATH_MSG_ERROR( myname << ": found '" << word << "' instead of 'DATE'" );
- return StatusCode::FAILURE;
- }
- input >> word >> time;
- if ( word != "TIME" ) {
- ATH_MSG_ERROR( myname << ": found '" << word << "' instead of 'TIME'" );
- return StatusCode::FAILURE;
- }
-
- // read and skip header cards
- int nheader;
- input >> word >> nheader;
- if ( word != "HEADERS" ) {
- ATH_MSG_ERROR( myname << ": found '" << word << "' instead of 'HEADERS'" );
- return StatusCode::FAILURE;
- }
- std::string restofline;
- getline( input, restofline );
- for ( int i = 0; i < nheader; i++ ) {
- std::string header;
- getline( input, header );
- }
-
- // read zone definitions
- int nzone;
- input >> word >> nzone;
- if ( word != "ZONES" ) {
- ATH_MSG_ERROR( myname << ": found '" << word << "' instead of 'ZONES'" );
- return StatusCode::FAILURE;
- }
- std::vector<int> jz(nzone);
-
- std::vector<int> nz(nzone);
- std::vector<int> jr(nzone);
-
- std::vector<int> nr(nzone);
- std::vector<int> jphi(nzone);
-
- std::vector<int> nphi(nzone);
- std::vector<int> jbs(nzone);
-
- std::vector<int> nbs(nzone);
- std::vector<int> jcoil(nzone);
-
- std::vector<int> ncoil(nzone);
- std::vector<int> jfield(nzone);
-
- std::vector<int> nfield(nzone);
- std::vector<int> jaux(nzone);
-
- std::vector<int> naux(nzone);
-
- for ( int i = 0; i < nzone; i++ )
- {
- int id;
- int nrep;
- int map; // unused
- double z1;
- double z2;
- double r1;
- double r2;
- double phi1;
- double phi2;
- int nzrphi0; // unused
- double tol; // unused
- int mzn;
- int mxsym;
- int mrefl;
- int mback; // unused
- double qz;
- double qr;
- double qphi; // unused
- double bscale;
- input >> id >> nrep;
- if ( version == 6 ) input >> map;
- input >> z1 >> z2 >> nz[i]
- >> r1 >> r2 >> nr[i]
- >> phi1 >> phi2 >> nphi[i]
- >> nzrphi0 >> tol
- >> jbs[i] >> nbs[i]
- >> jcoil[i] >> ncoil[i]
- >> jz[i] >> jr[i] >> jphi[i]
- >> jfield[i] >> nfield[i]
- >> mzn;
- if ( version == 6 ) input >> mxsym;
- input >> mrefl >> mback
- >> jaux[i] >> naux[i]
- >> qz >> qr >> qphi >> bscale;
- if ( id >= 0 ) { // remove dummy zone
- z1 *= CLHEP::meter;
- z2 *= CLHEP::meter;
- r1 *= CLHEP::meter;
- r2 *= CLHEP::meter;
- phi1 *= CLHEP::degree;
- phi2 *= CLHEP::degree;
- bscale *= CLHEP::tesla;
- BFieldZone zone( id, z1, z2, r1, r2, phi1, phi2, bscale );
- m_zone.push_back(zone);
- }
- }
-
- // read Biot-Savart data
- int nbiot;
- input >> word >> nbiot;
- if ( word != "BIOT" ) {
- ATH_MSG_ERROR( myname << ": found '" << word << "' instead of 'BIOT'" );
- return StatusCode::FAILURE;
- }
- std::vector<BFieldCond> bslist;
- for ( int i = 0; i < nbiot; i++ ) {
- char dummy; // unused
- char cfinite;
- double xyz1[3];
-
- double xyz2[3];
- double phirot; // unused
- double curr;
- input >> dummy >> cfinite
- >> xyz1[0] >> xyz1[1] >> xyz1[2]
- >> xyz2[0] >> xyz2[1] >> xyz2[2]
- >> phirot >> curr;
- bool finite = ( cfinite == 'T' );
- for ( int j = 0; j < 3; j++ ) {
- xyz1[j] *= CLHEP::meter;
- if ( finite ) xyz2[j] *= CLHEP::meter;
- }
- BFieldCond bs( finite, xyz1, xyz2, curr );
- bslist.push_back(bs);
- }
- // attach them to the zones
- for ( unsigned i = 0; i < m_zone.size(); i++ ) {
- // copy the range that belongs to this zone
- for ( int j = 0; j < nbs[i]; j++ ) {
- // Fortran -> C conversion requires "-1"
- m_zone[i].appendCond( bslist[jbs[i]+j-1] );
- }
- }
-
- // read and skip coil data
- int nc;
- input >> word >> nc;
- if ( word != "COIL" ) {
- ATH_MSG_ERROR( myname << ": found '" << word << "' instead of 'COIL'" );
- return StatusCode::FAILURE;
- }
- getline( input, restofline );
- for ( int i = 0; i < nc; i++ ) {
- std::string coildata;
- getline( input, coildata );
- }
-
- // read and skip auxiliary array = list of subzones
- int nauxarr;
- input >> word >> nauxarr;
- if ( word != "AUXARR" ) {
- ATH_MSG_ERROR( myname << ": found '" << word << "' instead of 'AUXARR'" );
- return StatusCode::FAILURE;
- }
- if ( version == 6 ) input >> word; // skip 'T'
- for ( int i = 0; i < nauxarr; i++ ) {
- int aux;
- input >> aux;
- }
-
- // read mesh definition
- int nmesh;
- input >> word >> nmesh;
- if ( word != "MESH" ) {
- ATH_MSG_ERROR( myname << ": found '" << word << "' instead of 'MESH'" );
- return StatusCode::FAILURE;
- }
- std::vector<double> meshlist;
- for ( int i = 0; i < nmesh; i++ ) {
- double mesh;
- input >> mesh;
- meshlist.push_back(mesh);
- }
- // attach them to the zones
- for ( unsigned i = 0; i < m_zone.size(); i++ ) {
- m_zone[i].reserve( nz[i], nr[i], nphi[i] );
- for ( int j = 0; j < nz[i]; j++ ) {
- m_zone[i].appendMesh( 0, meshlist[jz[i]+j-1]*CLHEP::meter );
- }
- for ( int j = 0; j < nr[i]; j++ ) {
- m_zone[i].appendMesh( 1, meshlist[jr[i]+j-1]*CLHEP::meter );
- }
- for ( int j = 0; j < nphi[i]; j++ ) {
- m_zone[i].appendMesh( 2, meshlist[jphi[i]+j-1] );
- }
- }
-
- // read field values
- int nf;
-
- int nzlist;
- std::string ftype;
-
- std::string bytype;
- input >> word >> nf >> nzlist >> ftype >> bytype;
- if ( word != "FIELD" ) {
- ATH_MSG_ERROR( myname << ": found '" << word << "' instead of 'FIELD'" );
- return StatusCode::FAILURE;
- }
- if ( ftype != "I2PACK" ) {
- ATH_MSG_ERROR( myname << ": found '" << ftype << "' instead of 'I2PACK'" );
- return StatusCode::FAILURE;
- }
- if ( bytype != "FBYTE" ) {
- ATH_MSG_ERROR( myname << ": found '" << bytype << "' instead of 'FBYTE'" );
- return StatusCode::FAILURE;
- }
- // read zone by zone
- for ( int i = 0; i < nzlist; i++ ) {
- int izone;
-
- int idzone;
-
- int nfzone;
- input >> izone >> idzone >> nfzone;
- izone--; // fortran -> C++
- if ( idzone != m_zone[izone].id() ) {
- ATH_MSG_ERROR( myname << ": zone id " << idzone << " != " << m_zone[izone].id() );
- return StatusCode(2);
- }
-
- std::vector<int> data[3];
-
- // for field data in 2 bytes
- for ( int j = 0; j < 3; j++ ) { // repeat z, r, phi
- int ierr = read_packed_data( input, data[j] );
- if ( ierr != 0 ) return StatusCode(ierr);
- for ( int k = 0; k < nfzone; k++ ) {
- // recover sign
- data[j][k] = ( data[j][k]%2==0 ) ? data[j][k]/2 : -(data[j][k]+1)/2;
- // second-order diff
- if ( k >= 2 ) data[j][k] += 2*data[j][k-1] - data[j][k-2];
- }
- }
- // store
- for ( int k = 0; k < nfzone; k++ ) {
- BFieldVector<short> B( data[0][k], data[1][k], data[2][k] );
- m_zone[izone].appendField( B );
- }
-
- // skip fbyte
- char c;
- while ( input.good() ) {
- input >> c;
- if ( input.eof() || c == '}' ) break;
- }
- }
-
- // build the LUTs and ZR zone
- buildLUT();
- buildZR();
-
- return StatusCode::SUCCESS;
-}
-
-//
-// wrire the map to a ROOT file
-//
-void MagField::AtlasFieldSvc::writeMap( TFile* rootfile ) const
-{
- if ( rootfile == nullptr ) return; // no file
- if ( !rootfile->cd() ) return; // could not make it current directory
- // define the tree
- TTree* tree = new TTree( "BFieldMap", "BFieldMap version 5" );
- TTree* tmax = new TTree( "BFieldMapSize", "Buffer size information" );
- int id;
- double zmin;
-
- double zmax;
-
- double rmin;
-
- double rmax;
-
- double phimin;
-
- double phimax;
- double bscale;
- int ncond;
- bool *finite;
- double *p1x;
-
- double *p1y;
-
- double *p1z;
-
- double *p2x;
-
- double *p2y;
-
- double *p2z;
- double *curr;
- int nmeshz;
-
- int nmeshr;
-
- int nmeshphi;
- double *meshz;
-
- double *meshr;
-
- double *meshphi;
- int nfield;
- short *fieldz;
-
- short *fieldr;
-
- short *fieldphi;
-
- // prepare arrays - need to know the maximum sizes
- unsigned maxcond(0);
-
- unsigned maxmeshz(0);
-
- unsigned maxmeshr(0);
-
- unsigned maxmeshphi(0);
-
- unsigned maxfield(0);
- for ( unsigned i = 0; i < m_zone.size(); i++ ) {
- maxcond = std::max( maxcond, m_zone[i].ncond() );
- maxmeshz = std::max( maxmeshz, m_zone[i].nmesh(0) );
- maxmeshr = std::max( maxmeshr, m_zone[i].nmesh(1) );
- maxmeshphi = std::max( maxmeshphi, m_zone[i].nmesh(2) );
- maxfield = std::max( maxfield, m_zone[i].nfield() );
- }
- // store the maximum sizes
- tmax->Branch( "maxcond", &maxcond, "maxcond/i");
- tmax->Branch( "maxmeshz", &maxmeshz, "maxmeshz/i");
- tmax->Branch( "maxmeshr", &maxmeshr, "maxmeshr/i");
- tmax->Branch( "maxmeshphi", &maxmeshphi, "maxmeshphi/i");
- tmax->Branch( "maxfield", &maxfield, "maxfield/i");
- tmax->Fill();
- // prepare buffers
- finite = new bool[maxcond];
- p1x = new double[maxcond];
- p1y = new double[maxcond];
- p1z = new double[maxcond];
- p2x = new double[maxcond];
- p2y = new double[maxcond];
- p2z = new double[maxcond];
- curr = new double[maxcond];
- meshz = new double[maxmeshz];
- meshr = new double[maxmeshr];
- meshphi = new double[maxmeshphi];
- fieldz = new short[maxfield];
- fieldr = new short[maxfield];
- fieldphi = new short[maxfield];
- // define the tree branches
- tree->Branch( "id", &id, "id/I" );
- tree->Branch( "zmin", &zmin, "zmin/D" );
- tree->Branch( "zmax", &zmax, "zmax/D" );
- tree->Branch( "rmin", &rmin, "rmin/D" );
- tree->Branch( "rmax", &rmax, "rmax/D" );
- tree->Branch( "phimin", &phimin, "phimin/D" );
- tree->Branch( "phimax", &phimax, "phimax/D" );
- tree->Branch( "bscale", &bscale, "bscale/D" );
- tree->Branch( "ncond", &ncond, "ncond/I" );
- tree->Branch( "finite", finite, "finite[ncond]/O" );
- tree->Branch( "p1x", p1x, "p1x[ncond]/D" );
- tree->Branch( "p1y", p1y, "p1y[ncond]/D" );
- tree->Branch( "p1z", p1z, "p1z[ncond]/D" );
- tree->Branch( "p2x", p2x, "p2x[ncond]/D" );
- tree->Branch( "p2y", p2y, "p2y[ncond]/D" );
- tree->Branch( "p2z", p2z, "p2z[ncond]/D" );
- tree->Branch( "curr", curr, "curr[ncond]/D" );
- tree->Branch( "nmeshz", &nmeshz, "nmeshz/I" );
- tree->Branch( "meshz", meshz, "meshz[nmeshz]/D" );
- tree->Branch( "nmeshr", &nmeshr, "nmeshr/I" );
- tree->Branch( "meshr", meshr, "meshr[nmeshr]/D" );
- tree->Branch( "nmeshphi", &nmeshphi, "nmeshphi/I" );
- tree->Branch( "meshphi", meshphi, "meshphi[nmeshphi]/D" );
- tree->Branch( "nfield", &nfield, "nfield/I" );
- tree->Branch( "fieldz", fieldz, "fieldz[nfield]/S" );
- tree->Branch( "fieldr", fieldr, "fieldr[nfield]/S" );
- tree->Branch( "fieldphi", fieldphi, "fieldphi[nfield]/S" );
- // loop over zones to write
- for ( unsigned i = 0; i < m_zone.size(); i++ ) {
- const BFieldZone z = m_zone[i];
- id = z.id();
- zmin = z.zmin(); zmax = z.zmax();
- rmin = z.rmin(); rmax = z.rmax();
- phimin = z.phimin(); phimax = z.phimax();
- bscale = z.bscale();
- ncond = z.ncond();
- for ( int j = 0; j < ncond; j++ ) {
- const BFieldCond& c = z.cond(j);
- finite[j] = c.finite();
- p1x[j] = c.p1(0);
- p1y[j] = c.p1(1);
- p1z[j] = c.p1(2);
- p2x[j] = c.p2(0);
- p2y[j] = c.p2(1);
- p2z[j] = c.p2(2);
- curr[j] = c.curr();
- }
- nmeshz = z.nmesh(0);
- for ( int j = 0; j < nmeshz; j++ ) {
- meshz[j] = z.mesh(0,j);
- }
- nmeshr = z.nmesh(1);
- for ( int j = 0; j < nmeshr; j++ ) {
- meshr[j] = z.mesh(1,j);
- }
- nmeshphi = z.nmesh(2);
- for ( int j = 0; j < nmeshphi; j++ ) {
- meshphi[j] = z.mesh(2,j);
- }
- nfield = z.nfield();
- for ( int j = 0; j < nfield; j++ ) {
- const BFieldVector<short> f = z.field(j);
- fieldz[j] = f.z();
- fieldr[j] = f.r();
- fieldphi[j] = f.phi();
- }
- tree->Fill();
- }
- rootfile->Write();
- // clean up
- delete[] finite;
- delete[] p1x;
- delete[] p1y;
- delete[] p1z;
- delete[] p2x;
- delete[] p2y;
- delete[] p2z;
- delete[] curr;
- delete[] meshz;
- delete[] meshr;
- delete[] meshphi;
- delete[] fieldz;
- delete[] fieldr;
- delete[] fieldphi;
-}
-
-//
-// read the map from a ROOT file.
-// returns 0 if successful.
-//
-StatusCode MagField::AtlasFieldSvc::readMap( TFile* rootfile )
-{
- if ( rootfile == nullptr ) {
- // no file
- ATH_MSG_ERROR("readMap(): unable to read field map, no TFile given");
- return StatusCode::FAILURE;
- }
- if ( !rootfile->cd() ) {
- // could not make it current directory
- ATH_MSG_ERROR("readMap(): unable to cd() into the ROOT field map TFile");
- return StatusCode::FAILURE;
- }
- // open the tree
- TTree* tree = (TTree*)rootfile->Get("BFieldMap");
- if ( tree == nullptr ) {
- // no tree
- ATH_MSG_ERROR("readMap(): TTree 'BFieldMap' does not exist in ROOT field map");
- return StatusCode::FAILURE;
- }
- int id;
- double zmin;
-
- double zmax;
-
- double rmin;
-
- double rmax;
-
- double phimin;
-
- double phimax;
- double bscale;
- int ncond;
- bool *finite;
- double *p1x;
-
- double *p1y;
-
- double *p1z;
-
- double *p2x;
-
- double *p2y;
-
- double *p2z;
- double *curr;
- int nmeshz;
-
- int nmeshr;
-
- int nmeshphi;
- double *meshz;
-
- double *meshr;
-
- double *meshphi;
- int nfield;
- short *fieldz;
-
- short *fieldr;
-
- short *fieldphi;
- // define the fixed-sized branches first
- tree->SetBranchAddress( "id", &id );
- tree->SetBranchAddress( "zmin", &zmin );
- tree->SetBranchAddress( "zmax", &zmax );
- tree->SetBranchAddress( "rmin", &rmin );
- tree->SetBranchAddress( "rmax", &rmax );
- tree->SetBranchAddress( "phimin", &phimin );
- tree->SetBranchAddress( "phimax", &phimax );
- tree->SetBranchAddress( "bscale", &bscale );
- tree->SetBranchAddress( "ncond", &ncond );
- tree->SetBranchAddress( "nmeshz", &nmeshz );
- tree->SetBranchAddress( "nmeshr", &nmeshr );
- tree->SetBranchAddress( "nmeshphi", &nmeshphi );
- tree->SetBranchAddress( "nfield", &nfield );
- // prepare arrays - need to know the maximum sizes
- // open the tree of buffer sizes (may not exist in old maps)
- unsigned maxcond(0);
-
- unsigned maxmeshz(0);
-
- unsigned maxmeshr(0);
-
- unsigned maxmeshphi(0);
-
- unsigned maxfield(0);
- TTree* tmax = (TTree*)rootfile->Get("BFieldMapSize");
- if ( tmax != nullptr ) {
- tmax->SetBranchAddress( "maxcond", &maxcond );
- tmax->SetBranchAddress( "maxmeshz", &maxmeshz );
- tmax->SetBranchAddress( "maxmeshr", &maxmeshr );
- tmax->SetBranchAddress( "maxmeshphi", &maxmeshphi );
- tmax->SetBranchAddress( "maxfield", &maxfield );
- tmax->GetEntry(0);
- } else { // "BFieldMapSize" tree does not exist
- for ( int i = 0; i < tree->GetEntries(); i++ ) {
- tree->GetEntry(i);
- maxcond = std::max( maxcond, (unsigned)ncond );
- maxmeshz = std::max( maxmeshz, (unsigned)nmeshz );
- maxmeshr = std::max( maxmeshr, (unsigned)nmeshr );
- maxmeshphi = std::max( maxmeshphi, (unsigned)nmeshphi );
- maxfield = std::max( maxfield, (unsigned)nfield );
- }
- }
- finite = new bool[maxcond];
- p1x = new double[maxcond];
- p1y = new double[maxcond];
- p1z = new double[maxcond];
- p2x = new double[maxcond];
- p2y = new double[maxcond];
- p2z = new double[maxcond];
- curr = new double[maxcond];
- meshz = new double[maxmeshz];
- meshr = new double[maxmeshr];
- meshphi = new double[maxmeshphi];
- fieldz = new short[maxfield];
- fieldr = new short[maxfield];
- fieldphi = new short[maxfield];
- // define the variable length branches
- tree->SetBranchAddress( "finite", finite );
- tree->SetBranchAddress( "p1x", p1x );
- tree->SetBranchAddress( "p1y", p1y );
- tree->SetBranchAddress( "p1z", p1z );
- tree->SetBranchAddress( "p2x", p2x );
- tree->SetBranchAddress( "p2y", p2y );
- tree->SetBranchAddress( "p2z", p2z );
- tree->SetBranchAddress( "curr", curr );
- tree->SetBranchAddress( "meshz", meshz );
- tree->SetBranchAddress( "meshr", meshr );
- tree->SetBranchAddress( "meshphi", meshphi );
- tree->SetBranchAddress( "fieldz", fieldz );
- tree->SetBranchAddress( "fieldr", fieldr );
- tree->SetBranchAddress( "fieldphi", fieldphi );
-
- // reserve the space for m_zone so that it won't move as the vector grows
- m_zone.reserve( tree->GetEntries() );
- // read all tree and store
- for ( int i = 0; i < tree->GetEntries(); i++ ) {
- tree->GetEntry(i);
- BFieldZone z( id, zmin, zmax, rmin, rmax, phimin, phimax, bscale );
- m_zone.push_back(z);
- m_zone.back().reserve( nmeshz, nmeshr, nmeshphi );
- for ( int j = 0; j < ncond; j++ ) {
- double p1[3];
-
- double p2[3];
- p1[0] = p1x[j];
- p1[1] = p1y[j];
- p1[2] = p1z[j];
- p2[0] = p2x[j];
- p2[1] = p2y[j];
- p2[2] = p2z[j];
- BFieldCond cond( finite[j], p1, p2, curr[j] );
- m_zone.back().appendCond(cond);
- }
- for ( int j = 0; j < nmeshz; j++ ) {
- m_zone.back().appendMesh( 0, meshz[j] );
- }
- for ( int j = 0; j < nmeshr; j++ ) {
- m_zone.back().appendMesh( 1, meshr[j] );
- }
- for ( int j = 0; j < nmeshphi; j++ ) {
- m_zone.back().appendMesh( 2, meshphi[j] );
- }
- for ( int j = 0; j < nfield; j++ ) {
- BFieldVector<short> field( fieldz[j], fieldr[j], fieldphi[j] );
- m_zone.back().appendField( field );
- }
- }
- // clean up
- tree->Delete();
- delete[] finite;
- delete[] p1x;
- delete[] p1y;
- delete[] p1z;
- delete[] p2x;
- delete[] p2y;
- delete[] p2z;
- delete[] curr;
- delete[] meshz;
- delete[] meshr;
- delete[] meshphi;
- delete[] fieldz;
- delete[] fieldr;
- delete[] fieldphi;
- // build the LUTs
- buildLUT();
- buildZR();
-
- return StatusCode::SUCCESS;
-}
-
-//
-// utility function used by readMap()
-//
-int MagField::AtlasFieldSvc::read_packed_data( std::istream& input, std::vector<int>& data ) const
-{
- const std::string myname("BFieldMap::read_packed_data()");
-
- data.resize(0);
- char mode = 'u';
- char c;
- while ( input.good() ) {
- input >> c;
- if ( input.eof() ) return 0;
- else if (c == '}') { // end of record
- return 0;
- }
- else if (c == 'z') { // series of zeros
- int n;
- int ierr = read_packed_int( input, n );
- if ( ierr != 0 ) return ierr;
- for ( int i = 0; i < n; i++ ) {
- data.push_back(0);
- }
- }
- else if (c >= 'u' && c <= 'y') { // mode change
- mode = c;
- }
- else if (c <= ' ' || c > 'z') {
- ATH_MSG_ERROR( myname << ": unexpected letter '" << c << "' in input" );
- return 3;
- }
- else { // normal letter in the range '!' - 't'
- switch (mode) {
- case 'u':
- {
- input.putback( c );
- int n;
- int ierr = read_packed_int( input, n );
- if ( ierr != 0 ) return ierr;
- data.push_back(n);
- }
- break;
- case 'v':
- {
- int n = c - '!';
- for ( int i = 0; i < 4; i++ ) {
- input >> c;
- data.push_back(c - '!' + 84*(n%3));
- n = n/3;
- }
- }
- break;
- case 'w':
- data.push_back(c - '!');
- break;
- case 'x':
- {
- int n = c - '!';
- data.push_back(n/9);
- data.push_back(n%9);
- }
- break;
- case 'y':
- {
- int n = c - '!';
- data.push_back(n/27);
- data.push_back((n/9)%3);
- data.push_back((n/3)%3);
- data.push_back(n%3);
- }
- break;
- }
- }
- }
- return 0;
-}
-
-//
-// utility function used by read_packed_data()
-//
-int MagField::AtlasFieldSvc::read_packed_int( std::istream &input, int &n ) const
-{
- const std::string myname("BFieldMap::read_packed_int()");
- n = 0;
- char c;
- input >> c;
- while ( c >= '!' && c <= 'J') {
- n = 42*n + c - '!';
- input >> c;
- }
- if ( c >= 'K' && c <= 't' ) {
- n = 42*n + c - 'K';
- } else {
- ATH_MSG_ERROR( myname << ": unexpected letter '" << c << "' in input" );
- return 4;
- }
- return 0;
-}
-
-//
-// Search for the zone that contains a point (z, r, phi)
-// This is a linear-search version, used only to construct the LUT.
-//
-BFieldZone* MagField::AtlasFieldSvc::findZoneSlow( double z, double r, double phi )
-{
- for ( int j = m_zone.size()-1; j >= 0; --j ) {
- if ( m_zone[j].inside( z, r, phi ) ) return &m_zone[j];
- }
- return nullptr;
-}
-
-//
-// Build the look-up table used by FindZone().
-// Called by readMap()
-// It also calls buildLUT() for all zones.
-//
-void MagField::AtlasFieldSvc::buildLUT()
-{
- // make lists of (z,r,phi) edges of all zones
- for ( int j = 0; j < 3; j++ ) { // z, r, phi
- for ( unsigned i = 0; i < m_zone.size(); i++ ) {
- double e[2];
- e[0] = m_zone[i].min(j);
- e[1] = m_zone[i].max(j);
- for ( int k = 0; k < 2; k++ ) {
- // for the phi edge, fit into [-pi,pi]
- if ( j==2 && e[k] > M_PI ) e[k] -= 2.0*M_PI;
- m_edge[j].push_back(e[k]);
- }
- }
- // add -pi and +pi to phi, just in case
- m_edge[2].push_back(-M_PI);
- m_edge[2].push_back(M_PI);
- // sort
- sort( m_edge[j].begin(), m_edge[j].end() );
- // remove duplicates
- // must do this by hand to allow small differences due to rounding in phi
- int i = 0;
- for ( unsigned k = 1; k < m_edge[j].size(); k++ ) {
- if ( fabs(m_edge[j][i] - m_edge[j][k]) < 1.0e-6 ) continue;
- m_edge[j][++i] = m_edge[j][k];
- }
- m_edge[j].resize( i+1 );
- // because of the small differences allowed in the comparison above,
- // m_edge[][] values do not exactly match the m_zone[] boundaries.
- // we have to fix this up in order to avoid invalid field values
- // very close to the zone boundaries.
- for ( unsigned i = 0; i < m_zone.size(); i++ ) {
- for ( unsigned k = 0; k < m_edge[j].size(); k++ ) {
- if ( fabs(m_zone[i].min(j) - m_edge[j][k]) < 1.0e-6 ) {
- m_zone[i].adjustMin(j,m_edge[j][k]);
- }
- if ( fabs(m_zone[i].max(j) - m_edge[j][k]) < 1.0e-6 ) {
- m_zone[i].adjustMax(j,m_edge[j][k]);
- }
- }
- }
- }
- // build LUT for edge finding
- for ( int j = 0; j < 3; j++ ) { // z, r, phi
- // find the size of the smallest interval
- double width = m_edge[j].back() - m_edge[j].front();
- double q(width);
- for ( unsigned i = 0; i < m_edge[j].size()-1; i++ ) {
- q = std::min( q, m_edge[j][i+1] - m_edge[j][i] );
- }
- // find the number of cells in the LUT
- int n = int(width/q) + 1;
- q = width/(n+0.5);
- m_invq[j] = 1.0/q;
- n++;
- // fill the LUT
- int m = 0;
- for ( int i = 0; i < n; i++ ) { // index of LUT
- if ( q*i+m_edge[j].front() > m_edge[j][m+1] ) m++;
- m_edgeLUT[j].push_back(m);
- }
- }
- // store min/max for speedup
- m_zmin=m_edge[0].front();
- m_zmax=m_edge[0].back();
- m_rmax=m_edge[1].back();
- // build LUT for zone finding
- m_nz = m_edge[0].size() - 1;
- m_nr = m_edge[1].size() - 1;
- m_nphi = m_edge[2].size() - 1;
- m_zoneLUT.reserve( m_nz*m_nr*m_nphi );
- for ( int iz = 0; iz < m_nz; iz++ ) {
- double z = 0.5*(m_edge[0][iz]+m_edge[0][iz+1]);
- for ( int ir = 0; ir < m_nr; ir++ ) {
- double r = 0.5*(m_edge[1][ir]+m_edge[1][ir+1]);
- for ( int iphi = 0; iphi < m_nphi; iphi++ ) {
- double phi = 0.5*(m_edge[2][iphi]+m_edge[2][iphi+1]);
- const BFieldZone* zone = findZoneSlow( z, r, phi );
- m_zoneLUT.push_back( zone );
- }
- }
- }
- // build LUT in each zone
- for ( unsigned i = 0; i < m_zone.size(); i++ ) {
- m_zone[i].buildLUT();
- }
-}
-
-//
-// Build the z-r 2d map for fast solenoid field
-//
-void MagField::AtlasFieldSvc::buildZR()
-{
- // delete if previously allocated
- delete m_meshZR;
-
- // find the solenoid zone
- // solenoid zone always covers 100 < R < 1000, but not necessarily R < 100
- // so we search for the zone that contains a point at R = 200, Z = 0
- const BFieldZone *solezone = findZone( 0.0, 200.0, 0.0 );
-
- // instantiate the new ZR map with the same external coverage as the solenoid zone
- // make sure R = 0 is covered
- m_meshZR = new BFieldMeshZR( solezone->zmin(), solezone->zmax(), 0.0, solezone->rmax() );
-
- // reserve the right amount of memroy
- unsigned nmeshz = solezone->nmesh(0);
- unsigned nmeshr = solezone->nmesh(1);
- if ( solezone->rmin() > 0.0 ) nmeshr++;
- m_meshZR->reserve( nmeshz, nmeshr );
-
- // copy the mesh structure in z/r
- // take care of R = 0 first
- if ( solezone->rmin() > 0.0 ) {
- m_meshZR->appendMesh( 1, 0.0 );
- }
- // copy the rest
- for ( int i = 0; i < 2; i++ ) { // z, r
- for ( unsigned j = 0; j < solezone->nmesh(i); j++ ) { // loop over mesh points
- m_meshZR->appendMesh( i, solezone->mesh( i, j ) );
- }
- }
-
- // loop through the mesh and compute the phi-averaged field
- for ( unsigned iz = 0; iz < m_meshZR->nmesh(0); iz++ ) { // loop over z
- double z = m_meshZR->mesh( 0, iz );
- for ( unsigned ir = 0; ir < m_meshZR->nmesh(1); ir++ ) { // loop over r
- double r = m_meshZR->mesh( 1, ir );
- const int nphi(200); // number of phi slices to average
- double Br = 0.0;
- double Bz = 0.0;
- for ( int iphi = 0; iphi < nphi; iphi++ ) {
- double phi = double(iphi)/double(nphi)*2.*M_PI;
- double xyz[3];
- xyz[0] = r*cos(phi);
- xyz[1] = r*sin(phi);
- xyz[2] = z;
- double B[3];
- solezone->getB( xyz, B, nullptr );
- Br += B[0]*cos(phi) + B[1]*sin(phi);
- Bz += B[2];
- }
- Br *= 1.0/double(nphi);
- Bz *= 1.0/double(nphi);
- m_meshZR->appendField( BFieldVectorZR( Bz, Br ) );
- }
- }
-
- // build the internal LUT
- m_meshZR->buildLUT();
-}
-
-//
-// Approximate memory footprint
-//
-int MagField::AtlasFieldSvc::memSize() const
-{
- int size = 0;
- size += sizeof(BFieldCache);
- size += sizeof(BFieldCacheZR);
- for ( unsigned i = 0; i < m_zone.size(); i++ ) {
- size += m_zone[i].memSize();
- }
- for ( int i = 0; i < 3; i++ ) {
- size += sizeof(double)*m_edge[i].capacity();
- size += sizeof(int)*m_edgeLUT[i].capacity();
- }
- size += sizeof(BFieldZone*)*m_zoneLUT.capacity();
- if ( m_meshZR ) {
- size += m_meshZR->memSize();
- }
- return size;
-}
-
diff --git a/MagneticField/MagFieldServices/src/components/MagFieldServices_entries.cxx b/MagneticField/MagFieldServices/src/components/MagFieldServices_entries.cxx
index 00ceb76cc36..db6ac48d80e 100644
--- a/MagneticField/MagFieldServices/src/components/MagFieldServices_entries.cxx
+++ b/MagneticField/MagFieldServices/src/components/MagFieldServices_entries.cxx
@@ -1,9 +1,7 @@
-#include "MagFieldServices/AtlasFieldSvc.h"
#include "MagFieldServices/AtlasFieldMapCondAlg.h"
#include "MagFieldServices/AtlasFieldCacheCondAlg.h"
#include "MagFieldServices/H8FieldSvc.h"
-DECLARE_COMPONENT( MagField::AtlasFieldSvc )
DECLARE_COMPONENT( MagField::AtlasFieldMapCondAlg )
DECLARE_COMPONENT( MagField::AtlasFieldCacheCondAlg )
DECLARE_COMPONENT( MagField::H8FieldSvc )
--
GitLab