From 0f68497aba442fa47021aebf8691eb0dbfd94df3 Mon Sep 17 00:00:00 2001
From: Markus Frank <Markus.Frank@cern.ch>
Date: Fri, 17 Mar 2023 16:46:00 +0100
Subject: [PATCH] Address issue https://github.com/AIDASoft/DD4hep/issues/1073

---
 DDCore/CMakeLists.txt                         |   5 +
 DDCore/include/DD4hep/FieldTypes.h            |   6 +-
 DDCore/src/FieldTypes.cpp                     |  49 ++--
 DDCore/src/Fields.cpp                         |   4 +-
 DDCore/src/graphics/DrawBField.cpp            | 269 ++++++++++++++++++
 DDCore/src/plugins/Compact2Objects.cpp        |   7 +-
 DDCore/src/plugins/DetectorFields.cpp         |   5 +-
 .../ClientTests/compact/QuadrupoleField.xml   |  54 +++-
 8 files changed, 371 insertions(+), 28 deletions(-)
 create mode 100644 DDCore/src/graphics/DrawBField.cpp

diff --git a/DDCore/CMakeLists.txt b/DDCore/CMakeLists.txt
index 021d25ca5..1d5621f60 100644
--- a/DDCore/CMakeLists.txt
+++ b/DDCore/CMakeLists.txt
@@ -107,6 +107,11 @@ endif()
 # Generate DDCore plugins---------------------------------------------------------
 dd4hep_add_plugin(DDCorePlugins SOURCES src/plugins/*.cpp USES DDCore)
 
+# Graphics DDCore plugins---------------------------------------------------------
+dd4hep_add_plugin(DDCoreGraphics
+  SOURCES src/graphics/*.cpp
+  USES DDCore ROOT::Core ROOT::Rint ROOT::Gui ROOT::Hist ROOT::HistPainter ROOT::ROOTHistDraw)
+
 # This plugins depends on the ROOT GDML readers. Hence, extra library
 IF(TARGET ROOT::Gdml)
   dd4hep_print("|++> Found Root::GDML: Creating DDGDMLPlugins")
diff --git a/DDCore/include/DD4hep/FieldTypes.h b/DDCore/include/DD4hep/FieldTypes.h
index d35aa5bf4..01b74f1e6 100644
--- a/DDCore/include/DD4hep/FieldTypes.h
+++ b/DDCore/include/DD4hep/FieldTypes.h
@@ -155,8 +155,12 @@ namespace dd4hep {
     Coefficents  skews       { };
     /// Boundary volume (optional)
     Solid        volume      { };
-    /// Position transformation of the field
+    /// Position transformation of the field. Only stored here for reference
     Transform3D  transform   { };
+    /// Inverse position transformation of the field
+    Transform3D  inverse     { };
+    /// The rotation part of the transformation. Need to rotate the field
+    Rotation3D   rotation    { };
     /// Constant Z field overlay
     double       B_z         { 0e0 };
 
diff --git a/DDCore/src/FieldTypes.cpp b/DDCore/src/FieldTypes.cpp
index 012eb3371..101b6dd81 100644
--- a/DDCore/src/FieldTypes.cpp
+++ b/DDCore/src/FieldTypes.cpp
@@ -49,7 +49,7 @@ void ConstantField::fieldComponents(const double* /* pos */, double* field) {
 SolenoidField::SolenoidField()
   : innerField(0), outerField(0), minZ(-INFINITY), maxZ(INFINITY), innerRadius(0), outerRadius(INFINITY)
 {
-  type = CartesianField::MAGNETIC;
+  field_type = CartesianField::MAGNETIC;
 }
 
 /// Compute  the field components at a given location and add to given field
@@ -67,7 +67,7 @@ void SolenoidField::fieldComponents(const double* pos, double* field) {
 
 /// Initializing constructor
 DipoleField::DipoleField() : zmax(INFINITY), zmin(-INFINITY), rmax(INFINITY) {
-  type = CartesianField::MAGNETIC;
+  field_type = CartesianField::MAGNETIC;
 }
 
 /// Compute  the field components at a given location and add to given field
@@ -96,35 +96,41 @@ namespace   {
   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;
+MultipoleField::MultipoleField()   {
+  field_type = CartesianField::MAGNETIC;
 }
 
 /// Compute  the field components at a given location and add to given field
 void MultipoleField::fieldComponents(const double* pos, double* field) {
-  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 ( (xx + yy + zz) < (3e0 - eps) )   {
+      flag |= FIELD_ROTATION_ONLY;
+    }
+    else  {
+      flag |= FIELD_POSITION_ONLY;
+      if ( (std::abs(dx) + std::abs(dy) + std::abs(dz)) < eps )
+	flag |= FIELD_IDENTITY;
+    }
+    this->inverse  = this->transform.Inverse();
+    this->transform.GetRotation(this->rotation);
+    this->transform.GetTranslation(this->translation);
   }
+  Transform3D::Point p, p0(pos[0],pos[1],pos[2]);
+  if      ( flag&FIELD_IDENTITY      ) p = p0;
+  else if ( flag&FIELD_POSITION_ONLY ) p = p0 - this->translation;
+  else      p = this->inverse * p0;
+
+  double x = p.X(), y = p.Y(), z = p.Z();
+  double coord[3] = {x, y, z};
+  ::printf("Pos: %+15.8e  %+15.8e  %+15.8e  Inverse: %+15.8e  %+15.8e  %+15.8e\n",
+	   pos[0]/dd4hep::cm,pos[1]/dd4hep::cm,pos[2]/dd4hep::cm, p.X()/dd4hep::cm, p.Y()/dd4hep::cm, p.Z()/dd4hep::cm);
 
   if ( 0 == volume.ptr() || volume->Contains(coord) )  {
     double bx = 0.0;
@@ -154,8 +160,9 @@ void MultipoleField::fieldComponents(const double* pos, double* field) {
     default:     // Error condition
       throw runtime_error("Invalid multipole field definition!");
     }
-    field[0] += bx;
-    field[1] += by;
-    field[2] += B_z;
+    Transform3D::Point f = this->rotation * Transform3D::Point(bx, by, B_z);
+    field[0] += f.X();
+    field[1] += f.Y();
+    field[2] += f.Z();
   }
 }
diff --git a/DDCore/src/Fields.cpp b/DDCore/src/Fields.cpp
index d66151e23..7bc8efeb1 100644
--- a/DDCore/src/Fields.cpp
+++ b/DDCore/src/Fields.cpp
@@ -33,6 +33,7 @@ namespace {
 
 /// Default constructor
 CartesianField::Object::Object() : TypedObject()  {
+  // The field_type MUST be overriden by the concrete sublass!
   field_type = UNKNOWN;
   InstanceCount::increment(this);
 }
@@ -126,8 +127,9 @@ void OverlayedField::add(CartesianField field) {
         o->field_type |= field.MAGNETIC;
         o->magnetic = (v.size() == 1) ? field : CartesianField();
       }
-      if (isMag || isEle)
+      if ( isMag || isEle )  {
         return;
+      }
       except("OverlayedField","add: Attempt to add an unknown field type.");
     }
     except("OverlayedField","add: Attempt to add an invalid object.");
diff --git a/DDCore/src/graphics/DrawBField.cpp b/DDCore/src/graphics/DrawBField.cpp
new file mode 100644
index 000000000..18b15f77f
--- /dev/null
+++ b/DDCore/src/graphics/DrawBField.cpp
@@ -0,0 +1,269 @@
+//==========================================================================
+//  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.
+//
+// Author     : M.Frank
+//
+//==========================================================================
+
+// Framework includes
+#include "DD4hep/Printout.h"
+#include "DD4hep/FieldTypes.h"
+#include "DD4hep/MatrixHelpers.h"
+#include "DD4hep/DetFactoryHelper.h"
+#include <cmath>
+
+using namespace dd4hep;
+
+#include <TH2.h>
+#include <TPad.h>
+#include <TRint.h>
+#include <TAxis.h>
+#include <TStyle.h>
+#include <TArrow.h>
+#include <TCanvas.h>
+#include <THistPainter.h>
+
+#include "Hoption.h"
+#include "Hparam.h"
+
+extern Hoption_t Hoption;
+extern Hparam_t  Hparam;
+
+
+namespace  {
+  void help_bfield(int argc, char** argv)   {
+    std::cout <<
+      "Usage: DD4hep_DrawBField -arg [-arg]                                 \n"
+      "  The plugin implicitly assumes that the geometry is already loaded. \n"
+      "     -area  <string>  Define area for which the B-field will be drawn\n"
+      "     -nx:   <number>  Number of bins in X direction                  \n"
+      "     -ny:   <number>  Number of bins in Y direction                  \n"
+      "     -dx:   <range>   Definition of the range in X                   \n"
+      "     -dy:   <range>   Definition of the range in Y                   \n"
+      "     -dz:   <range>   Definition of the range in Z                   \n"
+      "     Range definition as tuple: min,max                              \n"
+      "                 min:  lower boundary                                \n"
+      "                 max:  upper boundary                                \n"
+      "\tArguments given: " << arguments(argc,argv) << std::endl << std::flush;
+    ::exit(EINVAL);
+  }
+
+  struct range_t  {
+    double rmin=0e0, rmax=0e0;
+  };
+  struct field_t  {
+    Position  position;
+    Direction bfield;
+  };
+
+  range_t get_range(const std::string& s, int argc, char** argv)   {
+    int ival=0;
+    range_t range;
+    std::string val;
+    for(std::size_t i=0; i<s.length()+1; ++i)   {
+      char c = s.c_str()[i];
+      if ( c == ',' || i == s.length() )  {
+	switch(ival)   {
+	case 0:
+	  range.rmin = _toDouble(val);
+	  break;
+	case 1:
+	  range.rmax = _toDouble(val);
+	  break;
+	default:
+	  except("","+++ Too many value for range descriptor provided: %s", s.c_str());
+	  break;
+	}
+	val = "";
+	++ival;
+	continue;
+      }
+      val += c;
+    }
+    if ( ival != 2 )   {
+      help_bfield(argc, argv);
+    }
+    return range;
+  }
+
+  class MyHistPainter : public THistPainter {
+  public:
+    using THistPainter::THistPainter;
+    TH2* paintArrows(TH2* histo, const std::vector<field_t>& points)   {
+      fYaxis->SetTitle("Y coordinate [cm]");
+      fXaxis->SetTitle("X coordinate [cm]");
+      for(const auto& p : points)    {
+	double strength2 = p.bfield.mag2();
+	if ( strength2 > 1e-8 )   {
+	  double strength = sqrt(p.bfield.X()*p.bfield.X()+p.bfield.Y()*p.bfield.Y());
+	  int ix = fXaxis->FindBin(p.position.X());
+	  int iy = fYaxis->FindBin(p.position.Y());
+	  double xlo = fXaxis->GetBinLowEdge(ix);
+	  double xhi = fXaxis->GetBinLowEdge(ix+1);
+	  double ylo = fYaxis->GetBinLowEdge(iy);
+	  double yhi = fYaxis->GetBinLowEdge(iy+1);
+	  auto norm_field = p.bfield * ((xhi-xlo) / strength);
+	  double x1  = xlo + (xhi-xlo)/2.0;
+	  double x2  = x1 + norm_field.X();
+	  double y1  = ylo + (yhi-ylo)/2.0;
+	  double y2  = y1 + norm_field.Y();
+	  auto* arrow = new TArrow(x1, y1, x2, y2, 0.015, "|>");
+	  arrow->SetAngle(25);
+	  arrow->SetFillColor(kRed);
+	  arrow->SetLineColor(kRed);
+	  arrow->Draw();
+	  histo->Fill(p.position.X(), p.position.X(), strength);
+	  ::fprintf(stdout, "Arrow %+15.8e  %+15.8e  %+15.8e  %+15.8e    %+15.8e  %+15.8e\n",
+		    x1, y1, x2, y2, x2-x1, y2-y1);
+	  //delete arrow;
+	}
+      }
+      PaintPalette();
+      return histo;
+    }
+  };
+  
+  /// Basic entry point to interprete an XML document
+  /**
+   *  - The file URI to be opened 
+   *    is passed as first argument to the call.
+   *  - The processing hint (build type) is passed as optional 
+   *    second argument.
+   *
+   *  Factory: DD4hep_CompactLoader
+   *
+   *  \author  M.Frank
+   *  \version 1.0
+   *  \date    01/04/2014
+   */
+  static long draw_bfield(Detector& description, int argc, char** argv) {
+    std::vector<std::string> srange_x, srange_y, srange_z;
+    std::size_t nbin_x = 100, nbin_y = 100, nbin_z = 1;
+    std::string output;
+    double z_value = 0e0;
+    bool draw = false;
+
+    printout(WARNING,"DrawBField","This is a TEST. It does not function properly!");
+
+    for(int i = 0; i < argc && argv[i]; ++i)  {
+      if ( 0 == ::strncmp("-rx",argv[i],4) )
+	srange_x.push_back(argv[++i]);
+      else if ( 0 == ::strncmp("-ry",argv[i],3) )
+	srange_y.push_back(argv[++i]);
+      else if ( 0 == ::strncmp("-rz",argv[i],3) )
+	srange_z.push_back(argv[++i]);
+      else if ( 0 == ::strncmp("-nx",argv[i],3) )
+	nbin_x = _toULong(argv[++i]);
+      else if ( 0 == ::strncmp("-ny",argv[i],3) )
+	nbin_y = _toULong(argv[++i]);
+      else if ( 0 == ::strncmp("-nz",argv[i],3) )
+	nbin_z = _toULong(argv[++i]);
+      else if ( 0 == ::strncmp("-z",argv[i],2) )
+	z_value = _toDouble(argv[++i]);
+      else if ( 0 == ::strncmp("-draw",argv[i],4) )
+        draw = true;
+      else if ( 0 == ::strncmp("-output",argv[i],4) )
+        output = argv[++i];
+      else if ( 0 == ::strncmp("-help",argv[i],2) )
+	help_bfield(argc, argv);
+    }
+    if ( srange_x.empty() || srange_y.empty() )   {
+      help_bfield(argc, argv);
+    }
+    std::vector<range_t> range_x, range_y, range_z;
+    for(const auto& r : srange_x)  {
+      range_t rg = get_range(r, argc, argv);
+      range_x.push_back(rg);
+    }
+    for(const auto& r : srange_y)  {
+      range_t rg = get_range(r, argc, argv);
+      range_y.push_back(rg);
+    }
+    for(const auto& r : srange_z)  {
+      range_t rg = get_range(r, argc, argv);
+      range_z.push_back(rg);
+    }
+    
+    double ma = std::numeric_limits<double>::max();
+    range_t envelope_x  { ma, -ma };
+    range_t envelope_y  { ma, -ma };
+    range_t envelope_z  { ma, -ma };
+    for( const auto& range : range_x )   {
+      envelope_x.rmin = std::min(range.rmin, envelope_x.rmin);
+      envelope_x.rmax = std::max(range.rmax, envelope_x.rmax);
+    }
+    for( const auto& range : range_y )   {
+      envelope_y.rmin = std::min(range.rmin, envelope_y.rmin);
+      envelope_y.rmax = std::max(range.rmax, envelope_y.rmax);
+    }
+    for( const auto& range : range_z )   {
+      envelope_z.rmin = std::min(range.rmin, envelope_z.rmin);
+      envelope_z.rmax = std::max(range.rmax, envelope_z.rmax);
+    }
+
+    double dx = (envelope_x.rmax - envelope_x.rmin) / double(nbin_x);
+    double dy = (envelope_y.rmax - envelope_y.rmin) / double(nbin_y);
+    double dz = nbin_z == 1 ? 0e0 : (envelope_z.rmax - envelope_z.rmin) / double(nbin_z);
+    printf("Range(x) min:%4ld bins %+15.8e cm max:%+15.8e cm dx:%+15.8e cm\n",
+	   nbin_x, envelope_x.rmin/cm, envelope_x.rmax/cm, dx/cm);
+    printf("Range(y) min:%4ld bins %+15.8e cm max:%+15.8e cm dx:%+15.8e cm\n",
+	   nbin_y, envelope_y.rmin/cm, envelope_y.rmax/cm, dy/cm);
+    if ( nbin_z > 1 ) printf("Range(z) min:%4ld bins %+15.8e cm max:%+15.8e cm dx:%+15.8e cm\n",
+			     nbin_z, envelope_z.rmin/cm, envelope_z.rmax/cm, dz/cm);
+
+    FILE* out_file = ::fopen(output.empty() ? "/dev/null" : output.c_str(), "w");
+    ::fprintf(out_file,"#######################################################################################################\n");
+    ::fprintf(out_file,"      x[cm]            y[cm]            z[cm]          Bx[Tesla]        By[Tesla]        Bz[Tesla]     \n");
+    std::vector<field_t> field_values;
+    for( std::size_t i = 0; i < nbin_x; ++i )   {
+      float x = envelope_x.rmin + double(i)*dx + dx/2e0;
+      for( std::size_t j = 0; j < nbin_y; ++j )   {
+	float y = envelope_y.rmin + double(j)*dy + dy/2e0;
+	for( std::size_t k = 0; k < nbin_z; ++k )   {
+	  float z = nbin_z == 1 ? z_value : envelope_z.rmin + double(k)*dz + dz/2e0;
+	  field_t value;
+	  value.position = { x, y, z };
+	  value.bfield   = { 0e0, 0e0, 0e0 };
+	  value.bfield = description.field().magneticField(value.position);
+	  ::fprintf(out_file, " %+15.8e  %+15.8e  %+15.8e  %+15.8e  %+15.8e  %+15.8e\n",
+		 value.position.X()/cm, value.position.Y()/cm,  value.position.Z()/cm,
+		 value.bfield.X()/dd4hep::tesla, value.bfield.Y()/dd4hep::tesla, value.bfield.Z()/dd4hep::tesla);
+	  field_values.emplace_back(value);
+	}
+      }
+    }
+    ::fclose(out_file);
+    if ( draw )   {
+      if ( 0 == gApplication )  {
+	std::pair<int, char**> a(argc,argv);
+	gApplication = new TRint("DD4hepRint", &a.first, a.second);
+	printout(INFO,"DD4hepRint","++ Created ROOT interpreter instance for DD4hepUI.");
+      }
+      auto* histo = new TH2F("Bfield", "B-Field strength in Tesla",
+			    nbin_x, envelope_x.rmin, envelope_x.rmax,
+			    nbin_y, envelope_y.rmin, envelope_y.rmax);
+      MyHistPainter paint;
+      paint.SetHistogram(histo);
+      TCanvas* c1 = new TCanvas("B-Field");
+      c1->SetWindowSize(1000,1000);
+      //paint.Paint();
+      histo->SetStats(kFALSE);
+      //histo->Draw();
+      paint.paintArrows(histo, field_values);
+      histo->Draw("COLZ SAME");
+      gPad->SetGridx();
+      gPad->SetGridy();
+      if ( !gApplication->IsRunning() )  {
+	gApplication->Run();
+      }
+    }
+    return 1;
+  }
+}
+DECLARE_APPLY(DD4hep_DrawBField,draw_bfield)
diff --git a/DDCore/src/plugins/Compact2Objects.cpp b/DDCore/src/plugins/Compact2Objects.cpp
index bb4768fbf..6c979a911 100644
--- a/DDCore/src/plugins/Compact2Objects.cpp
+++ b/DDCore/src/plugins/Compact2Objects.cpp
@@ -137,7 +137,7 @@ static Ref_t create_ConstantField(Detector& /* description */, xml_h e) {
   xml_comp_t field(e), strength(e.child(_U(strength)));
   string t = e.attr<string>(_U(field));
   ConstantField* ptr = new ConstantField();
-  ptr->type = ::toupper(t[0]) == 'E' ? CartesianField::ELECTRIC : CartesianField::MAGNETIC;
+  ptr->field_type = ::toupper(t[0]) == 'E' ? CartesianField::ELECTRIC : CartesianField::MAGNETIC;
   ptr->direction.SetX(strength.x());
   ptr->direction.SetY(strength.y());
   ptr->direction.SetZ(strength.z());
@@ -188,6 +188,7 @@ static Ref_t create_SolenoidField(Detector& description, xml_h e) {
     ptr->minZ = c.attr<double>(_U(zmin));
   else
     ptr->minZ = -ptr->maxZ;
+  ptr->field_type = CartesianField::MAGNETIC;
   obj.assign(ptr, c.nameStr(), c.typeStr());
   return obj;
 }
@@ -219,6 +220,7 @@ static Ref_t create_DipoleField(Detector& /* description */, xml_h e) {
       val = _multiply<double>(coll.text(), mult);
     ptr->coefficents.emplace_back(val);
   }
+  ptr->field_type = CartesianField::MAGNETIC;
   obj.assign(ptr, c.nameStr(), c.typeStr());
   return obj;
 }
@@ -246,7 +248,7 @@ static Ref_t create_MultipoleField(Detector& description, xml_h e) {
     ptr->volume = xml::createShape(description, type, child);
   }
   ptr->B_z = bz;
-  ptr->transform = Transform3D(rot,pos).Inverse();
+  ptr->transform = Transform3D(rot,pos);
   for (xml_coll_t coll(c, _U(coefficient)); coll; ++coll, mult /= lunit) {
     xml_dim_t coeff = coll;
     if ( coll.hasAttr(_U(value)) )
@@ -257,6 +259,7 @@ static Ref_t create_MultipoleField(Detector& description, xml_h e) {
     val = coeff.skew(0.0) * mult;
     ptr->skews.emplace_back(val);
   }
+  ptr->field_type = CartesianField::MAGNETIC;
   obj.assign(ptr, c.nameStr(), c.typeStr());
   return obj;
 }
diff --git a/DDCore/src/plugins/DetectorFields.cpp b/DDCore/src/plugins/DetectorFields.cpp
index 0c7647129..87bc56319 100644
--- a/DDCore/src/plugins/DetectorFields.cpp
+++ b/DDCore/src/plugins/DetectorFields.cpp
@@ -12,9 +12,11 @@
 //==========================================================================
 
 // Framework includes
+#include "DD4hep/Printout.h"
 #include "DD4hep/FieldTypes.h"
 #include "DD4hep/MatrixHelpers.h"
 #include "DD4hep/DetFactoryHelper.h"
+#include <cmath>
 
 using namespace dd4hep;
 
@@ -96,8 +98,7 @@ static Handle<NamedObject> convert_multipole(Detector&, xml_h field, Handle<Name
   field.setAttr(_U(name), object->GetName());
   field.setAttr(_U(type), object->GetTitle());
   field.setAttr(_U(Z), ptr->B_z);
-  Transform3D inv = ptr->transform.Inverse();
-  detail::matrix::_decompose(inv, pos, rot);
+  detail::matrix::_decompose(ptr->transform, pos, rot);
   xml_elt_t x_pos = xml_elt_t(doc, _U(position));
   x_pos.setAttr(_U(x),pos.x());
   x_pos.setAttr(_U(y),pos.y());
diff --git a/examples/ClientTests/compact/QuadrupoleField.xml b/examples/ClientTests/compact/QuadrupoleField.xml
index 43b7bad68..e9dac494e 100644
--- a/examples/ClientTests/compact/QuadrupoleField.xml
+++ b/examples/ClientTests/compact/QuadrupoleField.xml
@@ -30,6 +30,7 @@
     <constant name="world_x"                value="20*cm"/>
     <constant name="world_y"                value="20*cm"/>
     <constant name="world_z"                value="20*cm"/>
+    <constant name="QC1_rmin"               value="2*cm"/>
   </define>
 
   <detectors>
@@ -40,8 +41,57 @@
   </detectors>
 
   <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_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_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>
 
-<fields>
+    <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>
+
+    -->
     <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"/>
@@ -73,5 +123,7 @@
       <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>
-- 
GitLab