FieldTypes.cpp

Go to the documentation of this file.

//==========================================================================
//  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
//
//==========================================================================
 
#include <DD4hep/FieldTypes.h>
#include <DD4hep/detail/Handle.inl>
#include <cmath>
 
using namespace dd4hep;
 
#ifndef INFINITY
#define INFINITY (numeric_limits<double>::max())
#endif
 
 
DD4HEP_INSTANTIATE_HANDLE(ConstantField);
DD4HEP_INSTANTIATE_HANDLE(SolenoidField);
DD4HEP_INSTANTIATE_HANDLE(DipoleField);
DD4HEP_INSTANTIATE_HANDLE(MultipoleField);
 
void ConstantField::fieldComponents(const double* pos, double* field) {
  if ( 0 == flag || nullptr == volume.ptr() )  {
    field[0] += direction.X();
    field[1] += direction.Y();
    field[2] += direction.Z();
  }
  else if ( flag&FIELD_LOCAL )  {
    Transform3D::Point p0(pos[0],pos[1],pos[2]);
    Transform3D::Point p = this->inverse_pos * p0;
    const Double_t coords[3] = {p.x(), p.y(), p.z()};
    if( volume->Contains(coords) )  {
      field[0] += direction.X();
      field[1] += direction.Y();
      field[2] += direction.Z();
    }
  }
}
 
SolenoidField::SolenoidField()
  : innerField(0), outerField(0), minZ(-INFINITY), maxZ(INFINITY), innerRadius(0), outerRadius(INFINITY)
{
  field_type = CartesianField::MAGNETIC;
}
 
void SolenoidField::fieldComponents(const double* pos, double* field) {
  double z = pos[2] ;
  //  std::cout << " field z=" << z << " maxZ=" << maxZ << " minZ = " << minZ << std::endl ;
  if( z > minZ && z < maxZ ){
    double radius = std::sqrt(pos[0] * pos[0] + pos[1] * pos[1]);
    if (radius < innerRadius)
      field[2] += innerField;
    else if (radius < outerRadius)
      field[2] += outerField;
  }
}
 
DipoleField::DipoleField() : zmax(INFINITY), zmin(-INFINITY), rmax(INFINITY) {
  field_type = CartesianField::MAGNETIC;
}
 
void DipoleField::fieldComponents(const double* pos, double* field) {
  double z = pos[2], r = std::sqrt(pos[0] * pos[0] + pos[1] * pos[1]);
  // Check if z coordinate is within dipole z bounds.
  if (z > zmin && z < zmax && r < rmax) {
    // Apply all coefficients to this z coordinate.
    double pp    = 1.0;
    double abs_z = fabs(z);
    double bx    = coefficents[0];
    for (size_t i = 1; i < coefficents.size(); ++i) {
      pp *= abs_z;
      bx += coefficents[i] * pp;
    }
    // Flip sign for negative z.
    if (z < 0)
      bx = -bx;
    // Add Bx to the input field.
    field[0] += bx;
  }
}
 
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;
  constexpr static unsigned char FIELD_HAS_AABB      = 1<<4;
}
 
MultipoleField::MultipoleField()   {
  field_type = CartesianField::MAGNETIC;
}
 
void MultipoleField::fieldComponents(const double* pos, double* field) {
  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;
    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);
    if ( volume.isValid() )   {
      // Obtain the bounding box from the shape. TGeoBBox is the base class
      // for most shapes, but not all.
      TGeoShape* shape = volume.ptr();
      auto* bbox = dynamic_cast<TGeoBBox*>(shape);
      if ( bbox )  {
        bbox->ComputeBBox();
        const double* orig = bbox->GetOrigin();
        double bdx = bbox->GetDX();
        double bdy = bbox->GetDY();
        double bdz = bbox->GetDZ();
        // Transform the local bbox center to world frame
        Transform3D::Point world_center =
          this->transform * Transform3D::Point(orig[0], orig[1], orig[2]);
        // World half-extents via rotation-matrix abs-sum formula (OBB -> AABB)
        double hx = std::abs(xx)*bdx + std::abs(xy)*bdy + std::abs(xz)*bdz;
        double hy = std::abs(yx)*bdx + std::abs(yy)*bdy + std::abs(yz)*bdz;
        double hz = std::abs(zx)*bdx + std::abs(zy)*bdy + std::abs(zz)*bdz;
        this->aabb_min[0] = world_center.X() - hx;
        this->aabb_max[0] = world_center.X() + hx;
        this->aabb_min[1] = world_center.Y() - hy;
        this->aabb_max[1] = world_center.Y() + hy;
        this->aabb_min[2] = world_center.Z() - hz;
        this->aabb_max[2] = world_center.Z() + hz;
        flag |= FIELD_HAS_AABB;
      }
    }
  }
  // AABB pre-filter: reject positions outside the world-frame bounding box
  if ( (flag&FIELD_HAS_AABB) &&
       (pos[0] < aabb_min[0] || pos[0] > aabb_max[0] ||
        pos[1] < aabb_min[1] || pos[1] > aabb_max[1] ||
        pos[2] < aabb_min[2] || pos[2] > aabb_max[2]) )
    return;
  // Transform the input position to the local frame of the field
  Transform3D::Point p, p0(pos[0],pos[1],pos[2]);
  if      ( flag&FIELD_IDENTITY      ) p = std::move(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};
  // Debug printout:
  //::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;
    double by = 0.0;
    double xy = x*y;
    double x2 = x*x;
    double y2 = y*y;
    switch(coefficents.size())  {
    case 4:      // Octupole momentum
      by += (1./6.) * ( coefficents[3] * (x2*x - 3.0*x*y2) + skews[3]*(y2*y - 3.0*x2*y) );
      bx += (1./6.) * ( coefficents[3] * (3.0*x2*y - y2*y) + skews[3]*(x2*x - 3.0*x*y2) );
      [[fallthrough]];
    case 3:      // Sextupole momentum:
      by +=  (1./2.) * ( coefficents[2] * (x2 - y2) - skews[2] * 2.0 * xy );
      bx +=  (1./2.) * ( coefficents[2] * 2.0 * xy + skews[2] * (x2 - y2) );
      [[fallthrough]];
    case 2:      // Quadrupole momentum:
      bx += coefficents[1] * y + skews[1]*x;
      by += coefficents[1] * x - skews[1]*y;
      [[fallthrough]];
    case 1:      // Dipole momentum:
      bx += skews[0];
      by += coefficents[0];
      [[fallthrough]];
    case 0:      // Nothing, but still valid
      break;
    default:     // Error condition
      throw std::runtime_error("Invalid multipole field definition!");
    }
    Transform3D::Point f = this->rotation * Transform3D::Point(bx, by, B_z);
    field[0] += f.X();
    field[1] += f.Y();
    field[2] += f.Z();
  }
}