diff --git a/DDCore/include/DD4hep/FieldTypes.h b/DDCore/include/DD4hep/FieldTypes.h index 963052fed79a4d1ad4253a18a43e25f0b3e88185..d35aa5bf43b884e13c06de6d7cb78e8b1f74fe60 100644 --- a/DDCore/include/DD4hep/FieldTypes.h +++ b/DDCore/include/DD4hep/FieldTypes.h @@ -149,12 +149,22 @@ namespace dd4hep { */ class MultipoleField : public CartesianField::Object { public: - Coefficents coefficents; - Coefficents skews; - Solid volume; - Transform3D transform; - double B_z; + /// Multi-pole coefficients + Coefficents coefficents { }; + /// Multi-pole skews + Coefficents skews { }; + /// Boundary volume (optional) + Solid volume { }; + /// Position transformation of the field + Transform3D transform { }; + /// Constant Z field overlay + double B_z { 0e0 }; + private: + /// The access to the field will be optimized. Remember properties. + unsigned char flag { 0 }; + /// Translation of the transformation + Transform3D::Point translation { }; public: /// Initializing constructor MultipoleField(); diff --git a/DDCore/include/DD4hep/Fields.h b/DDCore/include/DD4hep/Fields.h index bd1be4abccdcc4f8d961d4c0ca98246b449c8042..0a0d505c715bf578c8f2ba2f5c9c151f86dbba8b 100644 --- a/DDCore/include/DD4hep/Fields.h +++ b/DDCore/include/DD4hep/Fields.h @@ -40,7 +40,7 @@ namespace dd4hep { class CartesianField: public Handle<NamedObject> { public: enum FieldType { - UNKNOWN = 0, ELECTRIC = 0x1, MAGNETIC = 0x2 + UNKNOWN = 0, ELECTRIC = 0x1, MAGNETIC = 0x2, OVERLAY = 0x4, }; typedef std::map<std::string, std::map<std::string, std::string> > Properties; @@ -50,12 +50,24 @@ namespace dd4hep { * \version 1.0 * \ingroup DD4HEP_CORE */ - class Object: public NamedObject { + class TypedObject : public NamedObject { + public: + /// Field type + int field_type { UNKNOWN }; + /// Default constructor + using NamedObject::NamedObject; + }; + + /// Internal data class shared by all handles of a given type + /** + * \author M.Frank + * \version 1.0 + * \ingroup DD4HEP_CORE + */ + class Object: public TypedObject { public: /// Utility definition for concrete implementations typedef std::vector<double> Coefficents; - /// Field type - int type; /// Field extensions Properties properties; /// Default constructor @@ -85,7 +97,7 @@ namespace dd4hep { /// Access the field type int fieldType() const { - return data<Object>()->type; + return data<Object>()->field_type; } /// Access the field type (string) @@ -95,10 +107,10 @@ namespace dd4hep { bool changesEnergy() const; /// Returns the 3 field components (x, y, z). - void value(const Position& pos, Direction& field) const; // { value(pos,&field.x); } + void value(const Position& pos, Direction& field) const; /// Returns the 3 field components (x, y, z). - void value(const Position& pos, double* val) const; // { value(&pos.x,val); } + void value(const Position& pos, double* val) const; /// Returns the 3 field components (x, y, z). void value(const double* pos, double* val) const; @@ -126,7 +138,9 @@ namespace dd4hep { class OverlayedField: public Handle<NamedObject> { public: enum FieldType { - ELECTRIC = 0x1, MAGNETIC = 0x2 + ELECTRIC = CartesianField::ELECTRIC, + MAGNETIC = CartesianField::MAGNETIC, + OVERLAY = CartesianField::OVERLAY, }; typedef std::map<std::string, std::string> PropertyValues; typedef std::map<std::string, PropertyValues> Properties; @@ -137,15 +151,16 @@ namespace dd4hep { * \version 1.0 * \ingroup DD4HEP_CORE */ - class Object: public NamedObject { + class Object: public CartesianField::TypedObject { public: - int type; CartesianField electric; CartesianField magnetic; std::vector<CartesianField> electric_components; std::vector<CartesianField> magnetic_components; /// Field extensions Properties properties; + + public: /// Default constructor Object(); /// Default destructor @@ -162,8 +177,8 @@ namespace dd4hep { OverlayedField(const std::string& name); /// Access the field type - int type() const { - return data<Object>()->type; + int fieldType() const { + return data<Object>()->field_type; } /// Does the field change the energy of charged particles? @@ -173,50 +188,50 @@ namespace dd4hep { void add(CartesianField field); /// Returns the 3 electric field components (x, y, z) if many components are present - void combinedElectric(const Position& pos, double* field) const { - combinedElectric((const double*) &pos, field); - } + void combinedElectric(const Position& pos, double* field) const; /// Returns the 3 electric field components (x, y, z) if many components are present Direction combinedElectric(const Position& pos) const { - Direction field; - combinedElectric((const double*) &pos, (double*) &field); - return field; + double field[3] = { 0e0, 0e0, 0e0 }; + combinedElectric(pos, field); + return {field[0], field[1], field[2]}; } /// Returns the 3 electric field components (x, y, z) if many components are present - void combinedElectric(const double* pos, double* field) const; - - /// Returns the 3 magnetic field components (x, y, z) if many components are present - void combinedMagnetic(const Position& pos, double* field) const { - combinedMagnetic((const double*) &pos, field); + void combinedElectric(const double* pos, double* field) const { + combinedElectric(Position(pos[0], pos[1], pos[2]), field); } /// Returns the 3 magnetic field components (x, y, z) if many components are present - void combinedMagnetic(const double* pos, double* field) const; + void combinedMagnetic(const Position& pos, double* field) const; /// Returns the 3 magnetic field components (x, y, z) at a given position Direction combinedMagnetic(const Position& pos) const { - Direction field; - combinedMagnetic((const double*) &pos, (double*) &field); - return field; + double field[3] = { 0e0, 0e0, 0e0 }; + combinedMagnetic({pos.X(), pos.Y(), pos.Z()}, field); + return { field[0], field[1], field[2] }; } - /// Returns the 3 electric field components (x, y, z). - void electricField(const Position& pos, Direction& field) const { - electricField((const double*) &pos, (double*) &field); + /// Returns the 3 magnetic field components (x, y, z) if many components are present + void combinedMagnetic(const double* pos, double* field) const { + combinedMagnetic(Position(pos[0], pos[1], pos[2]), field); } /// Returns the 3 electric field components (x, y, z). - void electricField(const Position& pos, double* field) const { - electricField((double*) &pos, field); - } + void electricField(const Position& pos, double* field) const; /// Returns the 3 electric field components (x, y, z) at a given position Direction electricField(const Position& pos) const { - Direction field; - electricField((const double*) &pos, (double*) &field); - return field; + double field[3] = { 0e0, 0e0, 0e0 }; + electricField(pos, field); + return { field[0], field[1], field[2] }; + } + + /// Returns the 3 electric field components (x, y, z). + void electricField(const Position& pos, Direction& field) const { + double fld[3] = { 0e0, 0e0, 0e0 }; + electricField(Position(pos.X(), pos.Y(), pos.Z()), fld); + field = { fld[0], fld[1], fld[2] }; } /// Returns the 3 electric field components (x, y, z). @@ -226,37 +241,35 @@ namespace dd4hep { f.isValid() ? f.value(pos, field) : combinedElectric(pos, field); } - /// Returns the 3 magnetic field components (x, y, z). - void magneticField(const Position& pos, Direction& field) const { - magneticField((double*) &pos, (double*) &field); - } + /// Returns the 3 magnetic field components (x, y, z). + void magneticField(const Position& pos, double* field) const; - /// Returns the 3 electric field components (x, y, z) at a given position - Direction magneticField(const Position& pos) const { - Direction field; - magneticField((double*) &pos, (double*) &field); - return field; + /// Returns the 3 magnetic field components (x, y, z). + void magneticField(const double* pos, double* field) const { + magneticField(Position(pos[0], pos[1], pos[2]), field); } - /// Returns the 3 magnetic field components (x, y, z). - void magneticField(const Position& pos, double* field) const { - magneticField((double*) &pos, field); + /// Returns the 3 magnetic field components (x, y, z). + void magneticField(const double* pos, Direction& field) const { + double fld[3] = { 0e0, 0e0, 0e0 }; + magneticField(Position(pos[0], pos[1], pos[2]), fld); + field = { fld[0], fld[1], fld[2] }; } - /// Returns the 3 magnetic field components (x, y, z). - void magneticField(const double* pos, double* field) const { - field[0] = field[1] = field[2] = 0.0; - CartesianField f = data<Object>()->magnetic; - f.isValid() ? f.value(pos, field) : combinedMagnetic(pos, field); + /// Returns the 3 electric field components (x, y, z) at a given position + Direction magneticField(const Position& pos) const { + double field[3] = { 0e0, 0e0, 0e0 }; + magneticField(pos, field); + return { field[0], field[1], field[2] }; } /// Returns the 3 electric (val[0]-val[2]) and magnetic field components (val[3]-val[5]). - void electromagneticField(const Position& pos, double* val) const { - electromagneticField((double*) &pos, val); - } + void electromagneticField(const Position& pos, double* field) const; /// Returns the 3 electric (val[0]-val[2]) and magnetic field components (val[3]-val[5]). - void electromagneticField(const double* pos, double* val) const; + void electromagneticField(const double* pos, double* val) const { + electromagneticField(Position(pos[0], pos[1], pos[2]), val); + } /// Access to properties container Properties& properties() const; diff --git a/DDCore/src/FieldTypes.cpp b/DDCore/src/FieldTypes.cpp index d944f40860ea1f49a8ebaafb2a6a55a01a3d4d79..012eb3371dab5e567beb7d7545c458f18e7ea5f1 100644 --- a/DDCore/src/FieldTypes.cpp +++ b/DDCore/src/FieldTypes.cpp @@ -24,13 +24,13 @@ using namespace dd4hep; // fallthrough only exists from c++17 #if defined __has_cpp_attribute - #if __has_cpp_attribute(fallthrough) - #define ATTR_FALLTHROUGH [[fallthrough]] - #else - #define ATTR_FALLTHROUGH - #endif +#if __has_cpp_attribute(fallthrough) +#define ATTR_FALLTHROUGH [[fallthrough]] #else - #define ATTR_FALLTHROUGH +#define ATTR_FALLTHROUGH +#endif +#else +#define ATTR_FALLTHROUGH #endif DD4HEP_INSTANTIATE_HANDLE(ConstantField); @@ -91,6 +91,16 @@ void DipoleField::fieldComponents(const double* pos, double* field) { } } +namespace { + constexpr static unsigned char FIELD_INITIALIZED = 1<<0; + constexpr static unsigned char FIELD_IDENTITY = 1<<1; + constexpr static unsigned char FIELD_ROTATION_ONLY = 1<<2; + constexpr static unsigned char FIELD_POSITION_ONLY = 1<<3; + Transform3D::Point operator+(const Transform3D::Point& p0, const Transform3D::Point& p1) { + return Transform3D::Point(p0.X()+p1.X(), p0.Y()+p1.Y(),p0.Z()+p1.Z()); + } +} + /// Initializing constructor MultipoleField::MultipoleField() : coefficents(), skews(), volume(), transform(), B_z(0.0) { type = CartesianField::MAGNETIC; @@ -98,10 +108,24 @@ MultipoleField::MultipoleField() : coefficents(), skews(), volume(), transform() /// Compute the field components at a given location and add to given field void MultipoleField::fieldComponents(const double* pos, double* field) { - Transform3D::Point p = transform * Transform3D::Point(pos[0],pos[1],pos[2]); - //const Transform3D::Point::CoordinateType& c = p.Coordinates(); - double x=p.X(), y=p.Y(), z=p.Z(); - double coord[3] = {x,y,z}; + Transform3D::Point p0(pos[0],pos[1],pos[2]); + Transform3D::Point p = (flag&FIELD_IDENTITY) ? p0 + : (flag&FIELD_POSITION_ONLY) ? (p0 + translation) + : (transform * Transform3D::Point(pos[0],pos[1],pos[2])); + double x = p.X(), y = p.Y(), z = p.Z(); + double coord[3] = {x, y, z}; + + if ( 0 == flag ) { + constexpr static double eps = 1e-10; + double xx, xy, xz, dx, yx, yy, yz, dy, zx, zy, zz, dz; + transform.GetComponents(xx, xy, xz, dx, yx, yy, yz, dy, zx, zy, zz, dz); + flag |= FIELD_INITIALIZED; + flag = ((xx + yy + zz) < (3e0 - eps)) ? FIELD_ROTATION_ONLY : FIELD_POSITION_ONLY; + if ( (flag&FIELD_POSITION_ONLY) && (std::abs(dx) + std::abs(dy) + std::abs(dz)) < eps ) + flag |= FIELD_IDENTITY; + translation = Transform3D::Point(dx, dy, dz); + } + if ( 0 == volume.ptr() || volume->Contains(coord) ) { double bx = 0.0; double by = 0.0; diff --git a/DDCore/src/Fields.cpp b/DDCore/src/Fields.cpp index b642da53989ac90d445fadbabdc8821197b77eb3..d66151e23af8dc45d9950c12ac17ce0a6298a60d 100644 --- a/DDCore/src/Fields.cpp +++ b/DDCore/src/Fields.cpp @@ -12,6 +12,7 @@ //========================================================================== #include "DD4hep/Fields.h" +#include "DD4hep/Printout.h" #include "DD4hep/InstanceCount.h" #include "DD4hep/detail/Handle.inl" @@ -25,14 +26,14 @@ typedef OverlayedField::Object OverlayedFieldObject; DD4HEP_INSTANTIATE_HANDLE(OverlayedFieldObject); namespace { - void calculate_combined_field(vector<CartesianField>& v, const double* pos, double* field) { + void calculate_combined_field(vector<CartesianField>& v, const Position& pos, double* field) { for (const auto& i : v ) i.value(pos, field); } } /// Default constructor -CartesianField::Object::Object() - : NamedObject(), type(UNKNOWN) { +CartesianField::Object::Object() : TypedObject() { + field_type = UNKNOWN; InstanceCount::increment(this); } @@ -58,22 +59,27 @@ CartesianField::Properties& CartesianField::properties() const { /// Returns the 3 field components (x, y, z). void CartesianField::value(const Position& pos, Direction& field) const { - value(pos,(double*)&field); + double fld[3] = {0e0, 0e0, 0e0}; + double position[3] = {pos.X(), pos.Y(), pos.Z()}; + data<Object>()->fieldComponents(position, fld); + field = Direction(fld[0], fld[1], fld[2]); } /// Returns the 3 field components (x, y, z). -void CartesianField::value(const Position& pos, double* val) const { - value((double*)&pos,val); +void CartesianField::value(const Position& pos, double* field) const { + double position[3] = {pos.X(), pos.Y(), pos.Z()}; + data<Object>()->fieldComponents(position, field); } /// Returns the 3 field components (x, y, z). -void CartesianField::value(const double* pos, double* val) const { - data<Object>()->fieldComponents(pos, val); +void CartesianField::value(const double* pos, double* field) const { + data<Object>()->fieldComponents(pos, field); } /// Default constructor -OverlayedField::Object::Object() - : type(0), electric(), magnetic() { +OverlayedField::Object::Object() : TypedObject(), electric(), magnetic() +{ + field_type = CartesianField::OVERLAY; InstanceCount::increment(this); } @@ -83,9 +89,10 @@ OverlayedField::Object::~Object() { } /// Object constructor -OverlayedField::OverlayedField(const string& nam) - : Ref_t() { - assign(new Object(), nam, "overlay_field"); +OverlayedField::OverlayedField(const string& nam) : Ref_t() { + auto* obj = new Object(); + assign(obj, nam, "overlay_field"); + obj->field_type = CartesianField::OVERLAY; } /// Access to properties container @@ -95,7 +102,7 @@ OverlayedField::Properties& OverlayedField::properties() const { /// Does the field change the energy of charged particles? bool OverlayedField::changesEnergy() const { - int field = data<Object>()->type; + int field = data<Object>()->field_type; return CartesianField::ELECTRIC == (field & CartesianField::ELECTRIC); } @@ -103,45 +110,60 @@ bool OverlayedField::changesEnergy() const { void OverlayedField::add(CartesianField field) { if (field.isValid()) { Object* o = data<Object>(); - if (o) { - int typ = field.fieldType(); + if ( o ) { + int typ = field.fieldType(); bool isEle = field.ELECTRIC == (typ & field.ELECTRIC); bool isMag = field.MAGNETIC == (typ & field.MAGNETIC); if (isEle) { vector < CartesianField > &v = o->electric_components; v.emplace_back(field); - o->type |= field.ELECTRIC; - o->electric = v.size() == 1 ? field : CartesianField(); + o->field_type |= field.ELECTRIC; + o->electric = (v.size() == 1) ? field : CartesianField(); } if (isMag) { vector < CartesianField > &v = o->magnetic_components; v.emplace_back(field); - o->type |= field.MAGNETIC; - o->magnetic = v.size() == 1 ? field : CartesianField(); + o->field_type |= field.MAGNETIC; + o->magnetic = (v.size() == 1) ? field : CartesianField(); } if (isMag || isEle) return; - throw runtime_error("OverlayedField::add: Attempt to add an unknown field type."); + except("OverlayedField","add: Attempt to add an unknown field type."); } - throw runtime_error("OverlayedField::add: Attempt to add to an invalid object."); + except("OverlayedField","add: Attempt to add an invalid object."); } - throw runtime_error("OverlayedField::add: Attempt to add an invalid field."); + except("OverlayedField","add: Attempt to add an invalid field."); +} + +/// Returns the 3 magnetic field components (x, y, z). +void OverlayedField::magneticField(const Position& pos, double* field) const { + if ( isValid() ) { + field[0] = field[1] = field[2] = 0.0; + auto* obj = data<Object>(); + CartesianField f = obj->magnetic; + if ( f.isValid() ) + f.value(pos, field); + else + calculate_combined_field(obj->magnetic_components, pos, field); + return; + } + except("OverlayedField","add: Attempt to add an invalid field."); } /// Returns the 3 electric field components (x, y, z). -void OverlayedField::combinedElectric(const double* pos, double* field) const { +void OverlayedField::combinedElectric(const Position& pos, double* field) const { field[0] = field[1] = field[2] = 0.; calculate_combined_field(data<Object>()->electric_components, pos, field); } /// Returns the 3 magnetic field components (x, y, z). -void OverlayedField::combinedMagnetic(const double* pos, double* field) const { +void OverlayedField::combinedMagnetic(const Position& pos, double* field) const { field[0] = field[1] = field[2] = 0.; calculate_combined_field(data<Object>()->magnetic_components, pos, field); } /// Returns the 3 electric (val[0]-val[2]) and magnetic field components (val[3]-val[5]). -void OverlayedField::electromagneticField(const double* pos, double* field) const { +void OverlayedField::electromagneticField(const Position& pos, double* field) const { Object* o = data<Object>(); field[0] = field[1] = field[2] = 0.; calculate_combined_field(o->electric_components, pos, field); diff --git a/DDCore/src/RootDictionary.h b/DDCore/src/RootDictionary.h index 4088db8e18371324da274da7fd27936624ad3f46..432c312ca08ff75b1fc43de9d5710504dc1e2dcd 100644 --- a/DDCore/src/RootDictionary.h +++ b/DDCore/src/RootDictionary.h @@ -663,6 +663,7 @@ template class dd4hep::Handle<TNamed>; #pragma link C++ class std::map<dd4hep::VolumeID,dd4hep::VolumeManagerContext*>+; #pragma link C++ class dd4hep::CartesianField+; +#pragma link C++ class dd4hep::CartesianField::TypedObject+; #pragma link C++ class dd4hep::CartesianField::Object+; #pragma link C++ class dd4hep::Handle<dd4hep::CartesianField::Object>+; #pragma link C++ class dd4hep::OverlayedField+; diff --git a/examples/ClientTests/compact/QuadrupoleField.xml b/examples/ClientTests/compact/QuadrupoleField.xml new file mode 100644 index 0000000000000000000000000000000000000000..43b7bad68bfb9fd244664dc734197b64630f0038 --- /dev/null +++ b/examples/ClientTests/compact/QuadrupoleField.xml @@ -0,0 +1,77 @@ +<?xml version="1.0" encoding="UTF-8"?> +<lccdd> +<!-- #========================================================================== + # AIDA Detector description implementation + #========================================================================== + # Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN) + # All rights reserved. + # + # For the licensing terms see $DD4hepINSTALL/LICENSE. + # For the list of contributors see $DD4hepINSTALL/doc/CREDITS. + # + #========================================================================== +--> + + <info name="clic_sid_cdr" + title="CLIC Silicon Detector CDR" + author="Christian Grefe" + url="https://twiki.cern.ch/twiki/bin/view/CLIC/ClicSidCdr" + status="development" + version="$Id: compact.xml 1368 2014-11-03 20:15:27Z markus.frank@cern.ch $"> + <comment>The compact format for the CLIC Silicon Detector used for the conceptual design report</comment> + </info> + + <includes> + <gdmlFile ref="${DD4hepINSTALL}/DDDetectors/compact/elements.xml"/> + <gdmlFile ref="${DD4hepINSTALL}/DDDetectors/compact/materials.xml"/> + </includes> + + <define> + <constant name="world_x" value="20*cm"/> + <constant name="world_y" value="20*cm"/> + <constant name="world_z" value="20*cm"/> + </define> + + <detectors> + <detector id="1" name="Box" type="DD4hep_BoxSegment" vis="B2_vis"> + <material name="Air"/> + <box x="15*cm" y="15*cm" z="15*cm"/> + </detector> + </detectors> + + <fields> + +<fields> + <field name="QC1L1_field_0" type="MultipoleMagnet" Z="0.0*tesla"> + <position y="0*cm" x="5*cm" z="0*cm"/> + <rotation x="0" y="0" z="0.0"/> + <coefficient coefficient="0*tesla"/> + <coefficient coefficient="(-1)*(45.6)*(-0.273)/0.3*tesla/m"/> + <shape type="Tube" rmin="0.*cm" rmax="QC1_rmin" dz="2*cm" /> + </field> + + <field name="QC1L1_field_1" type="MultipoleMagnet" Z="0.0*tesla"> + <position y="5*cm" x="0*cm" z="0*cm"/> + <rotation x="0" y="pi/2." z="0.0"/> + <coefficient coefficient="0*tesla"/> + <coefficient coefficient="(-1)*(45.6)*(-0.273)/0.3*tesla/m"/> + <shape type="Tube" rmin="0.*cm" rmax="QC1_rmin" dz="2*cm" /> + </field> + + <field name="QC1L1_field_2" type="MultipoleMagnet" Z="0.0*tesla"> + <position y="0*cm" x="-5*cm" z="0*cm"/> + <rotation x="0" y="pi" z="0.0"/> + <coefficient coefficient="0*tesla"/> + <coefficient coefficient="(-1)*(45.6)*(-0.273)/0.3*tesla/m"/> + <shape type="Tube" rmin="0.*cm" rmax="QC1_rmin" dz="2*cm" /> + </field> + + <field name="QC1L1_field_3" type="MultipoleMagnet" Z="0.0*tesla"> + <position y="-5*cm" x="0*cm" z="0*cm"/> + <rotation x="0" y="pi*3./2." z="0.0"/> + <coefficient coefficient="0*tesla"/> + <coefficient coefficient="(-1)*(45.6)*(-0.273)/0.3*tesla/m"/> + <shape type="Tube" rmin="0.*cm" rmax="QC1_rmin" dz="2*cm" /> + </field> + </fields> +</lccdd>