Skip to content
Snippets Groups Projects
Commit 250d91c1 authored by Masahiro Morii's avatar Masahiro Morii Committed by Graeme Stewart
Browse files

'fix gcc4.7 warnings' (MagFieldUtils-00-01-06)

parent 95e4bb18
No related branches found
No related tags found
No related merge requests found
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
///////////////////////////////////////////////////////////////////
// G4FieldValidationAlg.h, (c) ATLAS Detector software
///////////////////////////////////////////////////////////////////
#ifndef MAGFIELDTESTBEDALG_H
#define MAGFIELDTESTBEDALG_H 1
// STL includes
#include <string>
#include <cmath>
// FrameWork includes
#include "GaudiKernel/ServiceHandle.h"
#include "AthenaBaseComps/AthAlgorithm.h"
#include "BFieldCore/AbstractMagneticField.h"
#include "BFieldAth/IMagFieldAthenaSvc.h"
// forward declarations
class ITHistSvc;
class TTree;
class G4Field;
namespace MagField {
class IMagFieldSvc;
}
namespace MagField {
/** @class MagFieldTestbedAlg
*/
class MagFieldTestbedAlg : public AthAlgorithm {
public:
//** Constructor with parameters */
MagFieldTestbedAlg( const std::string& name, ISvcLocator* pSvcLocator );
/** Destructor */
virtual ~MagFieldTestbedAlg();
/** Athena algorithm's interface method initialize() */
StatusCode initialize();
/** Athena algorithm's interface method execute() */
StatusCode execute();
/** Athena algorithm's interface method finalize() */
StatusCode finalize();
private:
bool checkWithReference(); //!< check current field with reference field
void getFieldValue(); //!< get Field value either by G4 or by MagFieldSvc
StatusCode fetchEnvironment(); //!< get environment either for g4 or for magFieldSvc
ServiceHandle<IMagFieldAthenaSvc> m_magFieldAthenaSvc; //!< service to get vanilla field svc
ServiceHandle<MagField::IMagFieldSvc> m_magFieldSvc; //!< service to get vanilla field svc
MagFieldAthena* p_MagField; //!< vanilla field service
const G4Field *p_g4field; //!< field service from G4
ServiceHandle<ITHistSvc> m_thistSvc; //!< the histogram service
std::string m_histStream; //!< THistSvc stream name
TTree *m_tree; //!< the ROOT tree containing the output
std::string m_treeName; //!< name of the Tree object
double m_halfX; //!< the half-lenght along x
double m_halfY; //!< the half-length along y
double m_halfZ; //!< the half-length along y
int m_stepsX; //!< the number of steps in X
int m_stepsY; //!< the number of steps in Y
int m_stepsZ; //!< the number of steps in Z
long m_numberOfReadings;//!< random readings instead of grid
// debug
int referenceCount;
std::string m_refFile; //!< reference field file name
std::string m_refTreeName; //!< TTree object in reference file
double m_absTolerance;//!< absolute tolerance in field against reference
double m_relTolerance;//!< relative tolerance in field against reference
double m_xyzt[4]; //!< stores the current xyzt position
double m_field[3]; //!< stores the field components
double m_deriv[9]; //!< stores derivatives TODO: currently only useMagFieldSvc = true
double m_explicitX; //!< if value to check is given explicitly
double m_explicitY; //!< if value to check is given explicitly
double m_explicitZ; //!< if value to check is given explicitly
bool m_complete; //!< validation completed already?
bool m_useG4Field; //!< G4
bool m_useOldMagFieldSvc; //!< MagFieldSvc
bool m_recordReadings; //!< record readings (for reproducability) or not (for speed tests)
bool m_useDerivatives; //!< only outer detector values will be tested
bool m_onlyCheckSolenoid; //!< only scan solenoid field
bool m_coordsAlongBeam; //!< coordinates on a beams in random directions
bool m_explicitCoords; //!< give values to check explicitly
};
} // end namespace MagField
#endif
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
///////////////////////////////////////////////////////////////////
// SolenoidTest.h, (c) ATLAS Detector software
///////////////////////////////////////////////////////////////////
#ifndef SOLENOIDTEST_H
#define SOLENOIDTEST_H 1
// STL includes
#include <string>
#include <cmath>
// FrameWork includes
#include "GaudiKernel/ServiceHandle.h"
#include "AthenaBaseComps/AthAlgorithm.h"
#include "BFieldCore/AbstractMagneticField.h"
#include "BFieldAth/IMagFieldAthenaSvc.h"
// forward declarations
class ITHistSvc;
class TTree;
namespace MagField {
class IMagFieldSvc;
}
namespace MagField {
class SolenoidTest : public AthAlgorithm {
public:
SolenoidTest( const std::string& name, ISvcLocator* pSvcLocator );
virtual ~SolenoidTest();
StatusCode initialize();
StatusCode execute();
StatusCode finalize();
private:
bool checkWithReference(); //!< check current field with reference field
void getFieldValue(); //!< get Field value either by G4 or by MagFieldSvc
StatusCode fetchEnvironment(); //!< get environment either for g4 or for magFieldSvc
ServiceHandle<IMagFieldAthenaSvc> m_magFieldAthenaSvc; //!< old field service
MagFieldAthena* m_magField; //!< old field service
ServiceHandle<MagField::IMagFieldSvc> m_magFieldSvc; //!< new field svc
ServiceHandle<ITHistSvc> m_thistSvc; //!< the histogram service
std::string m_histStream; //!< THistSvc stream name
TTree* m_tree; //!< the ROOT tree containing the output
std::string m_treeName; //!< name of the Tree object
double m_minR; //!< minimum radius to scan in mm
double m_maxR; //!< maximum radius to scan in mm
double m_minZ; //!< minimum z to scan in mm
double m_maxZ; //!< maximum z to scan in mm
int m_stepsR; //!< the number of steps in R
int m_stepsZ; //!< the number of steps in Z
int m_stepsPhi; //!< the number of steps in phi
double m_xyzt[4]; //!< stores the current xyzt position
double m_field[3]; //!< stores the field components
double m_fieldZR[3]; //!< stores the 2d field components
double m_fieldOld[3]; //!< stores the old field components
double m_deriv[9]; //!< stores derivatives
double m_derivZR[9]; //!< stores derivatives
double m_derivOld[9]; //!< stores derivatives
long m_event; //!< event counter
bool m_useFullField;
bool m_useFastField;
bool m_useOldField;
bool m_writeNtuple;
bool m_derivatives;
}; // end class SolenoidTest
} // end namespace MagField
#endif
package MagFieldUtils
author <elmar.ritsch@cern.ch>
manager Niels van Eldik <niels.van.eldik@cern.ch>
manager Robert Langenberg <robert.johannes.langenberg@cern.ch>
manager Masahiro Morii <masahiro_morii@harvard.edu>
manager Elmar Ritsch <elmar.ritsch@cern.ch>
#################################################################
# public use statements
#################################################################
public
use AtlasPolicy AtlasPolicy-*
use GaudiInterface GaudiInterface-* External
use AthenaBaseComps AthenaBaseComps-* Control
use BFieldCore BFieldCore-* MagneticField
use BFieldAth BFieldAth-* MagneticField
#################################################################
# private use statements
#################################################################
private
use Geant4 Geant4-* External
use AtlasROOT AtlasROOT-* External
use AtlasCLHEP AtlasCLHEP-* External
use MagFieldInterfaces MagFieldInterfaces-* MagneticField
public
library MagFieldUtils *.cxx components/*.cxx
apply_pattern component_library
apply_pattern declare_joboptions files="*.py"
apply_pattern declare_python_modules files="*.py"
private
# use this to activate debug info in this package:
#macro cppdebugflags '$(cppdebugflags_s)'
#macro_remove componentshr_linkopts "-Wl,-s"
# (0.) configure new MagFieldSvc
from AthenaCommon.AppMgr import ServiceMgr
from MagFieldServices.MagFieldServicesConf import MagField__AtlasFieldSvc
ServiceMgr += MagField__AtlasFieldSvc("AtlasFieldSvc");
#ServiceMgr.AtlasFieldSvc.FieldMapFile = '/afs/cern.ch/user/m/masahiro/public/BFieldMap/combinedMap.root'
from IOVDbSvc.CondDB import conddb
conddb.addOverride('/GLOBAL/BField/Map','BFieldMap-FullAsym-09-solTil3')
# (1.) setup the Athena algorithm
from AthenaCommon.AlgSequence import AlgSequence
topSequence = AlgSequence()
from MagFieldUtils.MagFieldUtilsConf import MagField__MagFieldTestbedAlg
topSequence += MagField__MagFieldTestbedAlg('MagFieldTestbedAlg')
topSequence.MagFieldTestbedAlg.useOldMagFieldSvc = True
topSequence.MagFieldTestbedAlg.numberOfReadings = 1000000
topSequence.MagFieldTestbedAlg.recordReadings = True
topSequence.MagFieldTestbedAlg.coordsAlongBeam = True
topSequence.MagFieldTestbedAlg.explicitCoords = True
topSequence.MagFieldTestbedAlg.explicitX = 4846.35
topSequence.MagFieldTestbedAlg.explicitY = 2037.96
topSequence.MagFieldTestbedAlg.explicitZ = -12625.6
topSequence.MagFieldTestbedAlg.useDerivatives = True
topSequence.MagFieldTestbedAlg.onlyCheckSolenoid = False
topSequence.MagFieldTestbedAlg.ReferenceFile = 'oldMagField.root'
# (2.) configure the THistStream
from GaudiSvc.GaudiSvcConf import THistSvc
if not hasattr(ServiceMgr, 'THistSvc'):
ServiceMgr += THistSvc("THistSvc")
ServiceMgr.THistSvc.Output = ["MagFieldTestbedAlg DATAFILE='comparedMagFields.root' OPT='RECREATE'"]
# (0.) configure the THistStream
from AthenaCommon.AppMgr import ServiceMgr
from GaudiSvc.GaudiSvcConf import THistSvc
if not hasattr(ServiceMgr, 'THistSvc'):
ServiceMgr += THistSvc("THistSvc")
ServiceMgr.THistSvc.Output = ["SolenoidTest DATAFILE='hurz.root' OPT='RECREATE'"]
#from MagFieldServices.MagFieldServicesConf import MagField__AtlasFieldSvc
#ServiceMgr += MagField__AtlasFieldSvc("AtlasFieldSvc");
# (1.) setup the SolenoidTest algorithm
from AthenaCommon.AlgSequence import AlgSequence
topSequence = AlgSequence()
from MagFieldUtils.MagFieldUtilsConf import MagField__SolenoidTest
testFull = 1000 # =1000 to measure the speed of the full 3d field
testFast = 1000 # =1000 to measure the speed of the fast 2d field
testOld = 1000 # =1000 to measure the speed of the old field
testHist = 0 # =100 to produce ROOT file to compare full vs. fast
if testFull:
topSequence += MagField__SolenoidTest('SolenoidTestFull')
topSequence.SolenoidTestFull.UseFullField = True
topSequence.SolenoidTestFull.UseFastField = False
topSequence.SolenoidTestFull.UseOldField = False
topSequence.SolenoidTestFull.WriteNtuple = False
topSequence.SolenoidTestFull.Derivatives = False
topSequence.SolenoidTestFull.StepsR = testFull
topSequence.SolenoidTestFull.StepsZ = testFull
topSequence.SolenoidTestFull.StepsPhi = testFull
if testFast:
topSequence += MagField__SolenoidTest('SolenoidTestFast')
topSequence.SolenoidTestFast.UseFullField = False
topSequence.SolenoidTestFast.UseFastField = True
topSequence.SolenoidTestFull.UseOldField = False
topSequence.SolenoidTestFast.WriteNtuple = False
topSequence.SolenoidTestFast.Derivatives = False
topSequence.SolenoidTestFast.StepsR = testFast
topSequence.SolenoidTestFast.StepsZ = testFast
topSequence.SolenoidTestFast.StepsPhi = testFast
if testOld:
topSequence += MagField__SolenoidTest('SolenoidTestOld')
topSequence.SolenoidTestOld.UseFullField = False
topSequence.SolenoidTestOld.UseFastField = False
topSequence.SolenoidTestOld.UseOldField = True
topSequence.SolenoidTestOld.WriteNtuple = False
topSequence.SolenoidTestOld.Derivatives = False
topSequence.SolenoidTestOld.StepsR = testOld
topSequence.SolenoidTestOld.StepsZ = testOld
topSequence.SolenoidTestOld.StepsPhi = testOld
if testHist:
topSequence += MagField__SolenoidTest('SolenoidTestHist')
topSequence.SolenoidTestHist.UseFullField = True
topSequence.SolenoidTestHist.UseFastField = True
topSequence.SolenoidTestHist.UseOldField = True
topSequence.SolenoidTestHist.WriteNtuple = True
topSequence.SolenoidTestHist.Derivatives = True
topSequence.SolenoidTestHist.StepsR = testHist
topSequence.SolenoidTestHist.StepsZ = testHist
topSequence.SolenoidTestHist.StepsPhi = testHist
theApp.EvtMax = 2
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
///////////////////////////////////////////////////////////////////
// MagFieldTestbedAlg.cxx, (c) ATLAS Detector software
///////////////////////////////////////////////////////////////////
/*
* Elmar Ritsch & Robert Langenberg
*
*/
// class include
#include "MagFieldUtils/MagFieldTestbedAlg.h"
// Framework
#include "GaudiKernel/ITHistSvc.h"
// CLHEP
#include "CLHEP/Vector/ThreeVector.h"
// MagneticField
#include "MagFieldInterfaces/IMagFieldSvc.h"
// G4includes
#include "G4TransportationManager.hh"
#include "G4FieldManager.hh"
#include "G4Field.hh"
// floating point precision
#include <limits>
//root
#include "TTree.h"
#include "TFile.h"
#include "TRandom3.h"
// performance test
#include <time.h>
#include <vector>
namespace {
const double solenoidRadius = 1075.;
const double solenoidZhalf = 2820.;
}
///////////////////////////////////////////////////////////////////
// Public methods:
///////////////////////////////////////////////////////////////////
// Constructor
////////////////
MagField::MagFieldTestbedAlg::MagFieldTestbedAlg(const std::string& name,
ISvcLocator* pSvcLocator) :
::AthAlgorithm(name, pSvcLocator), m_magFieldAthenaSvc(
"MagFieldAthenaSvc", name), m_magFieldSvc("AtlasFieldSvc",
name), m_thistSvc("THistSvc", name), m_histStream(
"MagFieldTestbedAlg"), m_tree(0), m_treeName("field"), m_halfX(
11999.), m_halfY(11999.), m_halfZ(22999.), m_stepsX(200), m_stepsY(
200), m_stepsZ(200), m_numberOfReadings(0), m_refFile(""), m_refTreeName(
"field"), m_absTolerance(1e-7), m_relTolerance(0.01), m_xyzt(), m_field(), m_explicitX(
0), m_explicitY(0), m_explicitZ(0), m_complete(false), m_useG4Field(
false), m_useOldMagFieldSvc(false), m_recordReadings(true), m_onlyCheckSolenoid(
false), m_coordsAlongBeam(false), m_explicitCoords(false) {
// histogram service
declareProperty("THistService", m_thistSvc, "The HistogramService");
declareProperty("HistogramStreamName", m_histStream = "MagFieldTestbedAlg",
"Name of the THistSvc output stream");
// TTree object name
declareProperty("ROOTTreeName", m_treeName = "field",
"Name of the TTree object in the output file.");
// boundaries for the magfield validation
declareProperty("HalfX", m_halfX = 11998., "half-length along x-direction");
declareProperty("HalfY", m_halfY = 11998., "half-length along y-direction");
declareProperty("HalfZ", m_halfZ = 22998., "half-length along z-direction");
// number of steps in each dimension (granularity)
declareProperty("StepsX", m_stepsX = 200,
"Number of steps along x-direction (granularity)");
declareProperty("StepsY", m_stepsY = 200,
"Number of steps along x-direction (granularity)");
declareProperty("StepsZ", m_stepsZ = 200,
"Number of steps along x-direction (granularity)");
// validate field against a reference file:
declareProperty("ReferenceFile", m_refFile = "",
"Filename of a reference file to compare the output with");
declareProperty("ReferenceTTreeName", m_refTreeName = "field",
"TTree object name in the reference file");
declareProperty("AbsTolerance", m_absTolerance = 1e-7,
"Numerical tolerance when comparing against reference.");
declareProperty("RelTolerance", m_relTolerance = 0.01,
"Numerical tolerance when comparing against reference.");
declareProperty("useG4Field", m_useG4Field);
declareProperty("useOldMagFieldSvc", m_useOldMagFieldSvc);
declareProperty("recordReadings", m_recordReadings = true);
declareProperty("numberOfReadings", m_numberOfReadings = 0);
declareProperty("useDerivatives", m_useDerivatives = false);
declareProperty("onlyCheckSolenoid", m_onlyCheckSolenoid = false);
declareProperty("coordsAlongBeam", m_coordsAlongBeam = false);
declareProperty("explicitCoords", m_explicitCoords = false);
declareProperty("explicitX", m_explicitX = 0.);
declareProperty("explicitY", m_explicitY = 0.);
declareProperty("explicitZ", m_explicitZ = 0.);
}
// Destructor
///////////////
MagField::MagFieldTestbedAlg::~MagFieldTestbedAlg() {
}
// Athena hook:
StatusCode MagField::MagFieldTestbedAlg::initialize() {
ATH_MSG_INFO("entering initialize()...");
if (m_useOldMagFieldSvc && m_magFieldAthenaSvc.retrieve().isFailure()) {
ATH_MSG_ERROR("Could not find MagFieldAthenaSvc");
return StatusCode::FAILURE;
} else if (!m_useG4Field && !m_useOldMagFieldSvc) {
//TODO: DUMMY reading to initialize newMagFieldSvc
m_xyzt[3] = 0.;
getFieldValue();
}
if (m_recordReadings || !(m_refFile.empty())) {
// retrieve the histogram service
if (m_thistSvc.retrieve().isSuccess()) {
// Create the prefix of histogram names for the THistSvc
std::string prefix = "/" + m_histStream + "/";
// the ROOT tree
m_tree = new TTree(m_treeName.c_str(), "Magnetic Field in Atlas");
m_tree->Branch("pos", &m_xyzt, "x/D:y/D:z/D");
m_tree->Branch("field", &m_field, "fx/D:fy/D:fz/D");
// G4Field doesn't allow to get derivatives
if (m_useDerivatives) {
m_tree->Branch("derivatives", &m_deriv,
"d1/D:d2/D:d3/D:d4/D:d5/D:d6/D:d7/D:d8/D:d9/D");
}
// register this ROOT TTree to the THistSvc
if (m_thistSvc->regTree(prefix + m_treeName, m_tree).isFailure()) {
ATH_MSG_ERROR("Unable to register TTree to THistSvc");
return StatusCode::FAILURE;
}
// failure in THistSvc retrieve
} else {
ATH_MSG_ERROR("Unable to retrieve HistogramSvc");
return StatusCode::FAILURE;
}
}
referenceCount = 0;
// success
ATH_MSG_INFO("end of initialize()");
return StatusCode::SUCCESS;
}
// Athena hook:
StatusCode MagField::MagFieldTestbedAlg::finalize() {
ATH_MSG_INFO("entering finalize()...");
ATH_MSG_INFO("end of finalize()");
return StatusCode::SUCCESS;
}
// Athena hook:
StatusCode MagField::MagFieldTestbedAlg::execute() {
// OUTPUT ONLY EXPLICIT COORDINATE THEN EXIT
if (m_explicitCoords) {
if (fetchEnvironment().isFailure()) {
ATH_MSG_ERROR("Unable to fetch the magnetic field!");
return StatusCode::FAILURE;
}
m_xyzt[0] = m_explicitX;
m_xyzt[1] = m_explicitY;
m_xyzt[2] = m_explicitZ;
getFieldValue();
ATH_MSG_INFO("Field Value at: " << m_explicitX << ", " << m_explicitY << ", " << m_explicitZ << "is: " << m_field[0] << ", " << m_field[1] << ", " << m_field[2] );
m_complete = true;
}
// if field validation was not done yet -> run it
if (!m_complete) {
if (fetchEnvironment().isFailure()) {
ATH_MSG_ERROR("Unable to fetch the magnetic field!");
return StatusCode::FAILURE;
}
// no reference file given?
if (m_refFile.empty()) {
//
// create normal validation output file
//
// set fixed time (will not be written to the ROOT tree)
std::vector < CLHEP::Hep3Vector > coordinatesToBeChecked;
double maxRadius = m_halfX;
double maxZ = m_halfZ;
if (m_onlyCheckSolenoid) {
maxRadius = solenoidRadius;
maxZ = solenoidZhalf;
}
if (m_numberOfReadings != 0 && !m_coordsAlongBeam) {
TRandom3 ran;
for (long i = 0; i < m_numberOfReadings; i++) {
CLHEP::Hep3Vector *temp;
temp = new CLHEP::Hep3Vector();
double maxRadius = m_halfX;
double maxZ = m_halfZ;
if (m_onlyCheckSolenoid) {
maxRadius = solenoidRadius;
maxZ = solenoidZhalf;
}
temp->setX(ran.Uniform(maxRadius * (-1), maxRadius));
temp->setY(ran.Uniform(maxRadius * (-1), maxRadius));
temp->setZ(ran.Uniform(maxZ * (-1), maxZ));
//only generate values within radius < 11.998; for onlyCheckSolenoid radius < 1075
double radius = sqrt(
pow(temp->getX(), 2) + pow(temp->getY(), 2));
if (radius > m_halfX) {
i--;
continue;
}
coordinatesToBeChecked.push_back(*temp);
}
}
if (m_coordsAlongBeam) {
TRandom3 rand;
CLHEP::Hep3Vector *temp;
temp = new CLHEP::Hep3Vector();
double r, phi, theta;
for (long i = 0; i < m_numberOfReadings;) {
r = 0.;
phi = rand.Uniform(3.14 * (-1), 3.14);
theta = rand.Uniform(0, 3.14);
while (r < maxRadius) {
temp->setX(r * sin(theta) * cos(phi));
temp->setY(r * sin(phi) * sin(theta));
temp->setZ(r * cos(phi));
r += 10.;
i++;
coordinatesToBeChecked.push_back(*temp);
}
}
}
// initialize performance timer
time_t seconds;
seconds = time(NULL);
if (m_numberOfReadings != 0) {
for (int i = 0; i < m_numberOfReadings; i++) {
m_xyzt[0] = coordinatesToBeChecked.at(i).getX();
m_xyzt[1] = coordinatesToBeChecked.at(i).getY();
m_xyzt[2] = coordinatesToBeChecked.at(i).getZ();
//TODO: delete random coordinates, this segfaults
// delete &randomCoordinates.at(i);
getFieldValue();
if (m_recordReadings) {
m_tree->Fill();
}
}
} else {
// loop over x
for (int stepx = 0; stepx < m_stepsX; stepx++) {
m_xyzt[0] = -maxRadius + (2 * maxRadius) / m_stepsX * stepx;
// loop over y
for (int stepy = 0; stepy < maxRadius; stepy++) {
m_xyzt[1] = -maxRadius
+ (2 * maxRadius) / maxRadius * stepy;
// loop over z
for (int stepz = 0; stepz < maxZ; stepz++) {
m_xyzt[2] = -maxZ + (2 * maxZ) / m_stepsZ * stepz;
// compute the field
getFieldValue();
// fill the ROOT Tree
if (m_recordReadings) {
m_tree->Fill();
}
}
}
}
}
seconds = time(NULL) - seconds;
if (m_numberOfReadings != 0) {
ATH_MSG_INFO(
"the reading of " << m_numberOfReadings
<< " random readings took " << seconds
<< "seconds.");
} else {
ATH_MSG_INFO(
"the reading of " << m_stepsX * m_stepsY * m_stepsZ
<< " values in a grid took " << seconds
<< "seconds.");
}
}
// reference file given
else {
//
// validate against reference file and create output file with positions
// from this file
//
ATH_MSG_INFO(
"will now run comparison against given reference file: "
<< m_refFile);
if (checkWithReference() == true)
ATH_MSG_INFO("comparison against reference file successful!");
else {
ATH_MSG_ERROR("comparison against reference file FAILED!!!!");
return StatusCode::FAILURE;
}
}
// never run the validation again
m_complete = true;
}
return StatusCode::SUCCESS;
}
StatusCode MagField::MagFieldTestbedAlg::fetchEnvironment() {
if (m_useG4Field) {
// use the transportation manager singleton to access the p_g4field
G4TransportationManager *transm =
G4TransportationManager::GetTransportationManager();
G4FieldManager *fieldm = (transm) ? transm->GetFieldManager() : 0;
p_g4field = (fieldm) ? fieldm->GetDetectorField() : 0;
if (!p_g4field) {
return StatusCode::FAILURE;
}
} else if (m_useOldMagFieldSvc) {
p_MagField = m_magFieldAthenaSvc->GetUpdatedMagFieldAthena();
}
return StatusCode::SUCCESS;
}
void MagField::MagFieldTestbedAlg::getFieldValue() {
if (m_useG4Field) {
p_g4field->GetFieldValue(m_xyzt, m_field);
} else if (m_useOldMagFieldSvc) {
// mm to cm unit conversion + double->float type conversion
float f_pos_in_cm[3];
for (int i = 0; i < 3; i++ ) {
f_pos_in_cm[i] = (float) (m_xyzt[i] * 0.1);
}
float f_bfield[12];
if (!m_useDerivatives) {
p_MagField->field_tesla_cm(f_pos_in_cm, f_bfield);
} else {
p_MagField->fieldGradient_tesla_cm(f_pos_in_cm, f_bfield);
}
for (int i = 0; i < 3; i++) {
m_field[i] = (double) f_bfield[i] * CLHEP::tesla;
}
if (m_useDerivatives) {
static const double tesla_cm( CLHEP::tesla / CLHEP::cm );
for (int i = 0; i < 9; i++) {
m_deriv[i] = (double) f_bfield[i+3] * tesla_cm;
}
}
} else {
// use new magnetic field service
m_magFieldSvc->getField(m_xyzt, m_field,
m_useDerivatives ? m_deriv : 0);
}
if (isnan(m_deriv[0]) || isnan(m_deriv[1]) || isnan(m_deriv[2])
|| isnan(m_deriv[3]) || isnan(m_deriv[4]) || isnan(m_deriv[5])
|| isnan(m_deriv[6]) || isnan(m_deriv[7]) || isnan(m_deriv[8])) {
ATH_MSG_INFO(
"at least one Derivative was NaN at: " << m_xyzt[0] << ", "
<< m_xyzt[1] << ", " << m_xyzt[2]);
}
}
// compare local field with reference file ROBERT: HISTOGRAM
bool MagField::MagFieldTestbedAlg::checkWithReference() {
// setup the reference TTree
TFile *refF = TFile::Open(m_refFile.c_str());
TTree *refT = (refF) ? (TTree*) refF->Get(m_refTreeName.c_str()) : 0;
//comparison with Masahiro's reference file
// std::string temp = "scan";
// TTree *refT = (refF) ? (TTree*) refF->Get(temp.c_str()) : 0;
if (!refT) {
ATH_MSG_ERROR(
"Unable to read the tree '" << m_refTreeName << "'"
<< "from file: " << m_refFile);
return false;
}
// enable all input branches
refT->SetBranchStatus("*", 1);
// setup the reference file branches
for (int i = 0; i < 4; i++)
m_xyzt[i] = 0.;
static double refField[3] = { 0., 0., 0. };
static double refDerivatives[9] = { 0., 0., 0., 0., 0., 0., 0., 0., 0. };
refT->SetBranchAddress("pos", &m_xyzt);
refT->SetBranchAddress("field", &refField);
if (m_useDerivatives) {
refT->SetBranchAddress("derivatives", &refDerivatives);
}
// number of reference TTree entries
Long64_t nentries = refT->GetEntries();
// setup the diff TTree
// TTree *diffT = new TTree("fieldDiff",
// "Magnetic Field Differences in AtlasG4");
static double fieldAbsDiff;
static double fieldRelDiff;
//
m_tree->Branch("ReferenceField", &refField, "x/D:y/D:z/D");
m_tree->Branch("fieldAbsDiff", &fieldAbsDiff, "diff/D");
m_tree->Branch("fieldRelDiff", &fieldRelDiff, "diff/D");
if (m_useDerivatives) {
m_tree->Branch("ReferenceDerivatives", &m_deriv,
"d1/D:d2/D:d3/D:d4/D:d5/D:d6/D:d7/D:d8/D:d9/D");
}
// register the diff ROOT TTree to the THistSvc
// std::string prefix = "/" + m_histStream + "/";
// if (m_thistSvc->regTree(prefix + "fieldDiff", diffT).isFailure()) {
// ATH_MSG_ERROR("Unable to register TTree to THistSvc");
// return false;
// }
// to keep track of the errors:
double biggestAbsDiff = 0.;
double biggestRelDiff = 0.;
// loop over all TTree entries
bool pass = true;
for (Long64_t i = 0; i < nentries; i++) {
// read the next entry in the reference ttree
// -> fill the m_xyzt and refField variables
refT->GetEntry(i);
// compute the field at the position given
// from the reference file
getFieldValue();
// compute the difference
double mFieldTotal = 0;
double refFieldTotal = 0;
for (int i = 0; i < 3; i++) {
mFieldTotal += m_field[i] * m_field[i];
refFieldTotal += refField[i] * refField[i];
}
mFieldTotal = sqrt(mFieldTotal);
refFieldTotal = sqrt(refFieldTotal);
fieldAbsDiff = fabs(mFieldTotal - refFieldTotal);
if (refFieldTotal != 0) {
fieldRelDiff = fieldAbsDiff / refFieldTotal;
}
if (m_useDerivatives) {
if (fabs((m_deriv[0] - refDerivatives[0]) / refDerivatives[0])
> m_relTolerance
|| fabs(
(m_deriv[1] - refDerivatives[1])
/ refDerivatives[1]) > m_relTolerance
|| fabs(
(m_deriv[2] - refDerivatives[2])
/ refDerivatives[2]) > m_relTolerance
|| fabs(
(m_deriv[3] - refDerivatives[3])
/ refDerivatives[3]) > m_relTolerance
|| fabs(
(m_deriv[4] - refDerivatives[4])
/ refDerivatives[4]) > m_relTolerance
|| fabs(
(m_deriv[5] - refDerivatives[5])
/ refDerivatives[5]) > m_relTolerance
|| fabs(
(m_deriv[6] - refDerivatives[6])
/ refDerivatives[6]) > m_relTolerance
|| fabs(
(m_deriv[7] - refDerivatives[7])
/ refDerivatives[7]) > m_relTolerance
|| fabs(
(m_deriv[8] - refDerivatives[8])
/ refDerivatives[8]) > m_relTolerance) {
ATH_MSG_INFO(
"referenceDerivatives 1-9: " << refDerivatives[0]
<< ", " << refDerivatives[1] << ", "
<< refDerivatives[2] << ", "
<< refDerivatives[3] << ", "
<< refDerivatives[4] << ", "
<< refDerivatives[5] << ", "
<< refDerivatives[6] << ", "
<< refDerivatives[7] << ", "
<< refDerivatives[8]);
ATH_MSG_INFO(
"derivatives 1-9: " << m_deriv[0] << ", "
<< m_deriv[1] << ", " << m_deriv[2] << ", "
<< m_deriv[3] << ", " << m_deriv[4] << ", "
<< m_deriv[5] << ", " << m_deriv[6] << ", "
<< m_deriv[7] << ", " << m_deriv[8]);
ATH_MSG_INFO(
"relative differences: "
<< ((m_deriv[0] - refDerivatives[0])
/ refDerivatives[0]) << ", "
<< ((m_deriv[1] - refDerivatives[1])
/ refDerivatives[1]) << ", "
<< ((m_deriv[2] - refDerivatives[2])
/ refDerivatives[2]) << ", "
<< ((m_deriv[3] - refDerivatives[3])
/ refDerivatives[3]) << ", "
<< ((m_deriv[4] - refDerivatives[4])
/ refDerivatives[4]) << ", "
<< ((m_deriv[5] - refDerivatives[5])
/ refDerivatives[5]) << ", "
<< ((m_deriv[6] - refDerivatives[6])
/ refDerivatives[6]) << ", "
<< ((m_deriv[7] - refDerivatives[7])
/ refDerivatives[7]) << ", "
<< ((m_deriv[8] - refDerivatives[8])
/ refDerivatives[8]));
ATH_MSG_INFO(
"coordinates: " << m_xyzt[0] << ", " << m_xyzt[1]
<< ", " << m_xyzt[2]);
referenceCount++;
}
}
// check if difference is smaller than one Gauss or if greater if < 1%
// only pass if below tolerances
if (fieldAbsDiff > m_absTolerance && fieldRelDiff > m_relTolerance) {
// if (pass) {
ATH_MSG_INFO(
"local reading" << m_field[0] << " " << m_field[1] << " "
<< m_field[2]);
ATH_MSG_INFO(
"reference reading" << refField[0] << " " << refField[1]
<< " " << refField[2]);
ATH_MSG_INFO(
"coordinates: " << m_xyzt[0] << ", " << m_xyzt[1] << ", "
<< m_xyzt[2]);
// }
pass = false;
}
if (fieldAbsDiff > biggestAbsDiff) {
biggestAbsDiff = fieldAbsDiff;
}
if (fieldRelDiff > biggestRelDiff) {
biggestRelDiff = fieldRelDiff;
}
// fill the magField ROOT Tree (using: m_xyzt, m_field)
m_tree->Fill();
// fill the diffT ROOT Tree (using: m_xyzt, refField)
// diffT->Fill();
}
// close the reference file
refF->Close();
ATH_MSG_INFO("number of values unequal: " << referenceCount);
ATH_MSG_INFO("number of readings from reference file: " << nentries);
ATH_MSG_INFO(
std::setprecision(20)
<< "1) biggest absolute difference in the mag field comparison: "
<< biggestAbsDiff << " (tolerance: " << m_absTolerance
<< ")." << endreq << "2) biggest relative difference: "
<< biggestRelDiff << " (tolerance: " << m_relTolerance
<< ")" << std::setprecision(-1));
// tests passed
return pass;
}
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
///////////////////////////////////////////////////////////////////
// SolenoidTest.cxx, (c) ATLAS Detector software
///////////////////////////////////////////////////////////////////
/*
* Masahiro Morii, Harvard University
*/
// class include
#include "MagFieldUtils/SolenoidTest.h"
// Framework
#include "GaudiKernel/ITHistSvc.h"
// CLHEP
#include "CLHEP/Vector/ThreeVector.h"
#include "CLHEP/Units/SystemOfUnits.h"
// MagneticField
#include "MagFieldInterfaces/IMagFieldSvc.h"
// floating point precision
#include <limits>
//root
#include "TTree.h"
#include "TFile.h"
#include "TRandom3.h"
// performance test
#include <time.h>
#include <vector>
///////////////////////////////////////////////////////////////////
// Public methods:
///////////////////////////////////////////////////////////////////
// Constructor
////////////////
MagField::SolenoidTest::SolenoidTest( const std::string& name, ISvcLocator* pSvcLocator ) :
::AthAlgorithm( name, pSvcLocator ),
m_magFieldAthenaSvc( "MagFieldAthenaSvc", name ),
m_magField(0),
m_magFieldSvc( "AtlasFieldSvc", name ),
m_thistSvc( "THistSvc", name ),
m_tree(0),
m_event(0)
{
// histogram service
declareProperty( "THistService", m_thistSvc, "The HistogramService" );
declareProperty( "HistogramStreamName", m_histStream = "SolenoidTest",
"Name of the THistSvc output stream" );
// TTree object name
declareProperty( "ROOTTreeName", m_treeName = "field",
"Name of the TTree object in the output file." );
// boundaries for the magfield validation
declareProperty( "MinimumR", m_minR = 0., "minimum R" );
declareProperty( "MaximumR", m_maxR = 1075., "maximum R" );
declareProperty( "MinimumZ", m_minZ = -2820., "minimum Z" );
declareProperty( "MaximumZ", m_maxZ = 2820., "maximum Z" );
// number of steps in each dimension (granularity)
declareProperty( "StepsR", m_stepsR = 200,
"Number of steps along radius (granularity)");
declareProperty( "StepsZ", m_stepsZ = 200,
"Number of steps along z (granularity)");
declareProperty( "StepsPhi", m_stepsPhi = 200,
"Number of steps along phi (granularity)");
declareProperty( "UseFullField", m_useFullField = true, "compute the full 3D field");
declareProperty( "UseFastField", m_useFastField = true, "compute the fast 2D field");
declareProperty( "UseOldField", m_useOldField = false, "compute the old field");
declareProperty( "WriteNtuple", m_writeNtuple = true, "produce the output ntuple");
declareProperty( "Derivatives", m_derivatives = false, "compute the derivatives");
}
// Destructor
///////////////
MagField::SolenoidTest::~SolenoidTest()
{;}
// Athena hook:
StatusCode MagField::SolenoidTest::initialize()
{
ATH_MSG_INFO("entering initialize()...");
if ( m_useOldField && m_magFieldAthenaSvc.retrieve().isFailure() ) {
ATH_MSG_ERROR( "Could not find MagFieldAthenaSvc" );
return StatusCode::FAILURE;
}
if ( m_writeNtuple ) {
// retrieve the histogram service
if ( m_thistSvc.retrieve().isSuccess() ) {
// Create the prefix of histogram names for the THistSvc
std::string prefix = "/" + m_histStream + "/";
// the ROOT tree
m_tree = new TTree( m_treeName.c_str(), "Magnetic field in the solenoid" );
m_tree->Branch( "pos", &m_xyzt, "x/D:y/D:z/D" );
if ( m_useFullField ) {
m_tree->Branch( "fieldFull", &m_field, "Bx/D:By/D:Bz/D" );
if ( m_derivatives ) {
m_tree->Branch( "deriv", &m_deriv,
"dBxdx/D:dBxdy/D:dBxdz/D:dBydx/D:dBydy/D:dBydz/D:dBzdx/D:dBzdy/D:dBzdz/D" );
}
}
if ( m_useFastField ) {
m_tree->Branch( "fieldFast", &m_fieldZR, "Bx2/D:By2/D:Bz2/D" );
if ( m_derivatives ) {
m_tree->Branch( "derivZR", &m_derivZR,
"dBxdx2/D:dBxdy2/D:dBxdz2/D:dBydx2/D:dBydy2/D:dBydz2/D:dBzdx2/D:dBzdy2/D:dBzdz2/D" );
}
}
if ( m_useOldField ) {
m_tree->Branch( "fieldOld", &m_fieldOld, "BxO/D:ByO/D:BzO/D" );
if ( m_derivatives ) {
m_tree->Branch( "derivOld", &m_derivOld,
"dBxdxO/D:dBxdyO/D:dBxdzO/D:dBydxO/D:dBydyO/D:dBydzO/D:dBzdxO/D:dBzdyO/D:dBzdzO/D" );
}
}
// register this ROOT TTree to the THistSvc
if ( m_thistSvc->regTree(prefix + m_treeName, m_tree).isFailure() ) {
ATH_MSG_ERROR("Unable to register TTree to THistSvc");
return StatusCode::FAILURE;
}
} else { // failure in THistSvc retrieve
ATH_MSG_ERROR("Unable to retrieve HistogramSvc");
return StatusCode::FAILURE;
}
}
// success
ATH_MSG_INFO("end of initialize()");
return StatusCode::SUCCESS;
}
// Athena hook:
StatusCode MagField::SolenoidTest::finalize() {
ATH_MSG_INFO("entering finalize()...");
ATH_MSG_INFO("end of finalize()");
return StatusCode::SUCCESS;
}
// Athena hook:
StatusCode MagField::SolenoidTest::execute() {
if ( m_useOldField ) {
m_magField = m_magFieldAthenaSvc->GetUpdatedMagFieldAthena();
}
if( m_event==0 ) { // event #0
// call field service once for initialization
m_xyzt[0] = m_xyzt[1] = m_xyzt[2] = 0;
if ( m_useFullField ) m_magFieldSvc->getField( m_xyzt, m_field, 0 );
if ( m_useFastField ) m_magFieldSvc->getFieldZR( m_xyzt, m_fieldZR, 0 );
// check the field status
ATH_MSG_INFO("New service solenoidOn = " << m_magFieldSvc->solenoidOn());
ATH_MSG_INFO("New service toroidOn = " << m_magFieldSvc->toroidOn());
if ( m_magField ) {
ATH_MSG_INFO("Old service getToroidBarrelOn = " << m_magField->getToroidBarrelOn());
ATH_MSG_INFO("Old service getToroidECTAOn = " << m_magField->getToroidECTAOn());
ATH_MSG_INFO("Old service getToroidECTCOn = " << m_magField->getToroidECTCOn());
ATH_MSG_INFO("Old service getSolenoidOn = " << m_magField->getSolenoidOn());
ATH_MSG_INFO("Old service getAllToroidOn = " << m_magField->getAllToroidOn());
ATH_MSG_INFO("Old service getFieldStatusOn = " << m_magField->getFieldStatusOn());
}
}
else if ( m_event==1 ) { // event #1
ATH_MSG_INFO( "solenoid field service test in progress." );
// initialize performance timer
time_t seconds;
seconds = time(NULL);
for ( int k = 0; k < m_stepsPhi; k++ ) { // loop over phi
double phi = 2.*M_PI*double(k)/double(m_stepsPhi);
for ( int j = 0; j < m_stepsZ; j++ ) { // loop over Z
m_xyzt[2] = (m_maxZ-m_minZ)*double(j)/double(m_stepsZ-1) + m_minZ;
for ( int i = 0; i < m_stepsR; i++ ) { // loop over R
double r = (m_maxR-m_minR)*double(i)/double(m_stepsR-1) + m_minR;
m_xyzt[0] = r*cos(phi);
m_xyzt[1] = r*sin(phi);
// compute the field
if ( m_useFullField ) {
if ( m_derivatives )
m_magFieldSvc->getField( m_xyzt, m_field, m_deriv );
else
m_magFieldSvc->getField( m_xyzt, m_field, 0 );
}
if ( m_useFastField ) {
if ( m_derivatives )
m_magFieldSvc->getFieldZR( m_xyzt, m_fieldZR, m_derivZR );
else
m_magFieldSvc->getFieldZR( m_xyzt, m_fieldZR, 0 );
}
// old field
if ( m_useOldField ) {
// unit conversion from mm to cm
float f_pos_in_cm[3] = { (float)(m_xyzt[0]*.1), (float)(m_xyzt[1]*.1), (float)(m_xyzt[2]*.1) };
float f_bfield[12];
if ( !m_derivatives ) {
m_magField->field_tesla_cm( f_pos_in_cm, f_bfield );
} else {
m_magField->fieldGradient_tesla_cm( f_pos_in_cm, f_bfield );
}
// unit conversion from kG to kT
for ( int i = 0; i < 3; i++ ) {
m_fieldOld[i] = (double)(f_bfield[i]*10.*CLHEP::kilogauss);
}
if ( m_derivatives ) {
// unit conversion from kG/mm to kT/mm
for ( int i = 0; i < 9; i++ ) {
m_derivOld[i] = (double)(f_bfield[i+3]*CLHEP::kilogauss);
}
}
}
// fill the ROOT Tree
if ( m_writeNtuple ) m_tree->Fill();
}
}
}
seconds = time(NULL) - seconds;
ATH_MSG_INFO( "accessing " << m_stepsR * m_stepsZ * m_stepsPhi
<< " values in a grid took " << seconds << " seconds." );
}
m_event++;
return StatusCode::SUCCESS;
}
#include "GaudiKernel/DeclareFactoryEntries.h"
#include "MagFieldUtils/MagFieldTestbedAlg.h"
#include "MagFieldUtils/SolenoidTest.h"
DECLARE_NAMESPACE_ALGORITHM_FACTORY( MagField , MagFieldTestbedAlg )
DECLARE_NAMESPACE_ALGORITHM_FACTORY( MagField , SolenoidTest )
DECLARE_FACTORY_ENTRIES( MagFieldUtils ) {
DECLARE_NAMESPACE_ALGORITHM( MagField , MagFieldTestbedAlg )
DECLARE_NAMESPACE_ALGORITHM( MagField , SolenoidTest )
}
#include "GaudiKernel/LoadFactoryEntries.h"
LOAD_FACTORY_ENTRIES( MagFieldUtils )
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment