TRungeFitter.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TRungeFitter.cxx,v 1.39 2022/01/28 22:14:13 maqm Exp $
//-----------------------------------------------------------------------------
// Filename : TRungeFitter.cc
// Section  : Tracking
// Owner    : Kenji Inami
// Email    : inami@bmail.kek.jp
//-----------------------------------------------------------------------------
// Description : A class to fit a TTrackBase object to a Runge Kutta track
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
 
#ifdef TRKRECO_DEBUG_DETAIL
#  ifndef TRKRECO_DEBUG
#    define TRKRECO_DEBUG
#  endif
#endif
#include "TrkReco/TRungeFitter.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/Service.h"
#include "TrkReco/TRunge.h"
#include <float.h>
#include <iostream>
// #define HEP_SHORT_NAMES
// #include "panther/panther.h"
// #include MDC_H
#include "AIDA/IHistogram1D.h"
#include "AIDA/IHistogram2D.h"
#include "CLHEP/Matrix/Matrix.h"
#include "CLHEP/Matrix/SymMatrix.h"
#include "CLHEP/Matrix/Vector.h"
#include "G4Geo/BesG4Geo.h"
#include "G4Geo/MdcG4Geo.h"
#include "G4Material.hh"
#include "G4Tubs.hh"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/IInterface.h"
#include "GaudiKernel/IMessageSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/Kernel.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/Service.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/StatusCode.h"
#include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
#include "MdcTables/MdcTables.h"
#include "TrkReco/TMLink.h"
#include "TrkReco/TTrackBase.h"
 
#include <time.h>
 
using CLHEP::HepMatrix;
using CLHEP::HepSymMatrix;
using CLHEP::HepVector;
/*
extern "C"
void
calcdc_driftdist_(int *,
                  int *,
                  int *,
                  float[3],
                  float[3],
                  float *,
                  float *,
                  float *);
 
extern "C"
void
calcdc_tof2_(int *, float *, float *, float *);
*/
 
TRungeFitter::TRungeFitter( const std::string& name )
    : TMFitter( name ), _sag( true ), _propagation( 1 ), _tof( false ), m_pmgnIMF( nullptr ) {
  // StatusCode scmgn = Gaudi::svcLocator()->service( "MagneticFieldSvc", m_pmgnIMF );
  // if ( scmgn != StatusCode::SUCCESS )
  // { std::cout << __FILE__ << " Unable to open Magnetic field service" << std::endl; }
}
TRungeFitter::TRungeFitter( const std::string& name, bool m_sag, int m_prop, bool m_tof )
    : TMFitter( name )
    , _sag( m_sag )
    , _propagation( m_prop )
    , _tof( m_tof )
    , m_pmgnIMF( nullptr ) {}
 
TRungeFitter::~TRungeFitter() {}
 
void TRungeFitter::sag( bool _in ) { _sag = _in; }
void TRungeFitter::propagation( int _in ) { _propagation = _in; }
void TRungeFitter::tof( bool _in ) { _tof = _in; }
int  TRungeFitter::fit( TTrackBase& tb ) const { return fit( tb, 0, -1 ); }
 
int TRungeFitter::fit( TTrackBase& tb, float t0Offset, int layer ) const {
  // Argument layer is used for calibration interface in which doca is calculated without
  // measured hit. liucy
  IMdcCalibFunSvc* l_mdcCalFunSvc;
  Gaudi::svcLocator()->service( "MdcCalibFunSvc", l_mdcCalFunSvc );
  // std::cout<<"TRungeFitter::fit  start"<<std::endl;
  //...Type check...
  if ( tb.objectType() != Runge ) return TFitUnavailable;
  TRunge& t = (TRunge&)tb;
  //...Already fitted ?...
  if ( t.fitted() && layer == -1 ) return TFitAlreadyFitted;
  //...Count # of hits...
  const AList<TMLink>& cores        = t.cores();
  unsigned             nCores       = cores.length();
  unsigned             nStereoCores = NStereoHits( cores );
  TMLink*              last         = NULL;
  int                  lyr          = 0;
  int                  layerid      = 0;
  for ( int i = 0; i < nCores; i++ )
  {
    layerid = ( *cores[i]->hit() ).wire()->layerId();
    if ( layerid >= lyr )
    {
      lyr  = layerid;
      last = cores[i];
    }
  }
 
  //...Check # of hits...
  // if ((nStereoCores < 2) || (nCores - nStereoCores < 3))
  //   return TFitErrorFewHits;
 
  //...Move pivot to ORIGIN...
  const HepPoint3D pivot_bak = t.pivot();
  t.pivot( ORIGIN );
 
  //...Setup...
  Vector a( 5 ), da( 5 );
  a = t.a();
  Vector          ainitial( a );
  Vector          dDda( 5 );
  Vector          dchi2da( 5 );
  SymMatrix       d2chi2d2a( 5, 0 );
  const SymMatrix zero5( 5, 0 );
  double          chi2;
  double          chi2Old = 0;
  for ( int i = 0; i < t.links().length(); i++ ) { chi2Old = chi2Old + t.links()[i]->pull(); }
  chi2Old = DBL_MAX;
  int err = 0;
 
  double    factor = 0.1;
  Vector    beta( 5 );
  SymMatrix alpha( 5, 0 );
  SymMatrix alpha2( 5, 0 );
 
  double( *d )[5] = new double[nCores][5];
  //  double (*d2)[5]=new double[nCores][5];
  Vector ea( 5 );
 
  float       tof;
  HepVector3D p;
  unsigned    i;
 
  double distance;
  double eDistance;
 
  // ea... init
  const double ea_init = 0.000001;
  ea[0]                = ea_init; // dr
  ea[1]                = ea_init; // phi0
  ea[2]                = ea_init; // kappa
  ea[3]                = ea_init; // dz
  ea[4]                = ea_init; // tanl
 
  // std::cout<<"TRF ::"<<a[0]<<","<<a[1]<<","<<a[2]<<","<<a[3]<<","<<a[4]<<std::endl;
 
  // long int lclock0=clock()/1000;
  // long int lclock=lclock0;
 
  Vector def_a( t.a() );
  Vector ta( def_a );
 
  //...Fitting loop...
  unsigned nTrial = 0;
  while ( nTrial < 100 )
  {
 
    //...Set up...
    chi2 = 0;
    for ( unsigned j = 0; j < 5; j++ ) dchi2da[j] = 0;
    d2chi2d2a = zero5;
 
    def_a = t.a();
 
    // #### loop for shifted helix parameter set ####
    for ( unsigned j = 0; j < 5; j++ )
    {
      ta = def_a;
      ta[j] += ea[j];
      t.a( ta );
      //...Loop with hits...
      i = 0;
      // std::cout<<"TRF:: cores="<<cores.length()<<std::endl;
      while ( TMLink* l = cores[i++] )
      {
        //...Cal. closest points...
        t.approach( *l, tof, p, _BesBeamPipeMaterials[0], _sag );
        const HepPoint3D& onTrack = l->positionOnTrack();
        const HepPoint3D& onWire  = l->positionOnWire();
        d[i - 1][j]               = ( onTrack - onWire ).mag();
        // std::cout<<"TRF:: i="<<i<<std::endl;
      } // end of loop with hits
      // lclock=clock()/1000;
      // std::cout<<"TRF  clock="<<lclock-lclock0<<std::endl;
      // lclock0=lclock;
    }
    /*
    for(int j=0;j<5;j++){
      ta=def_a;
      ta[j]-=ea[j];
      t.a(ta);
      //...Loop with hits...
      float tof_dummy;
      Vector3 p_dummy;
      unsigned i=0;
      while(TMLink* l = cores[i++]){
        //...Cal. closest points...
        t.approach(*l,tof_dummy,p_dummy,_sag);
        const HepPoint3D& onTrack=l->positionOnTrack();
        const HepPoint3D& onWire=l->positionOnWire();
        d2[i-1][j]=(onTrack-onWire).mag();
      }//end of loop with hits
    }
    */
    t.a( def_a );
 
    // #### original parameter set  and  calc. chi2 ####
    //...Loop with hits...
    i = 0;
    while ( TMLink* l = cores[i++] )
    {
      const TMDCWireHit& h = *l->hit();
      //...Cal. closest points...
      if ( t.approach( *l, tof, p, _BesBeamPipeMaterials[0], _sag ) < 0 )
      {
        std::cout << "TrkReco:TRF::  bad wire" << std::endl;
        continue;
      }
      int               layerId   = h.wire()->layerId();
      const HepPoint3D& onTrack   = l->positionOnTrack();
      const HepPoint3D& onWire    = l->positionOnWire();
      unsigned          leftRight = WireHitRight;
      if ( onWire.cross( onTrack ).z() < 0. ) leftRight = WireHitLeft;
 
      //...Obtain drift distance and its error...
      if ( ( t0Offset == 0. ) && ( _propagation == 0 ) && ( !_tof ) )
      {
        //...No correction...
        distance  = l->drift( leftRight );
        eDistance = h.dDrift( leftRight );
      }
      else
      {
        //...T0 and propagation corrections...
        int wire      = h.wire()->id();
        int wirelocal = h.wire()->localId();
        int side      = leftRight;
        if ( side == 0 ) side = 0;
        // maqm float tp[3] = {p.x(),p.y(),p.z()};
        double tp[3] = { p.x(), p.y(), p.z() };
        // maqm float x[3] = {onWire.x(), onWire.y(), onWire.z()};
        double x[3]      = { onWire.x(), onWire.y(), onWire.z() };
        double tprop     = l_mdcCalFunSvc->getTprop( layerId, onWire.z() * 10. );
        double timeWalk  = l_mdcCalFunSvc->getTimeWalk( layerId, h.reccdc()->adc );
        double T0        = l_mdcCalFunSvc->getT0( layerId, wirelocal );
        double drifttime = h.reccdc()->tdc - tof - tprop - timeWalk - T0;
        l->setDriftTime( drifttime );
        float        dist;
        float        edist;
        int          prop = _propagation;
        const double x0   = t.helix().pivot().x();
        const double y0   = t.helix().pivot().y();
        Hep3Vector   pivot_helix( x0, y0, 0 );
        const double dr    = t.helix().a()[0];
        const double phi0  = t.helix().a()[1];
        const double kappa = t.helix().a()[2];
        const double dz    = t.helix().a()[3];
        const double tanl  = t.helix().a()[4];
 
        const double Bz    = -1000 * ( getPmgnIMF()->getReferField() );
        const double alpha = 333.564095 / Bz;
 
        const double cox = x0 + dr * cos( phi0 ) - alpha * cos( phi0 ) / kappa;
        const double coy = y0 + dr * sin( phi0 ) - alpha * sin( phi0 ) / kappa;
        HepPoint3D   dir( onTrack.y() - coy, cox - onTrack.x(), 0 );
        double       pos_phi = onWire.phi();
        double       dir_phi = dir.phi();
        while ( pos_phi > pi ) { pos_phi -= pi; }
        while ( pos_phi < 0 ) { pos_phi += pi; }
        while ( dir_phi > pi ) { dir_phi -= pi; }
        while ( dir_phi < 0 ) { dir_phi += pi; }
        double entrangle = dir_phi - pos_phi;
        while ( entrangle > pi / 2 ) entrangle -= pi;
        while ( entrangle < ( -1 ) * pi / 2 ) entrangle += pi;
        dist =
            l_mdcCalFunSvc->driftTimeToDist( drifttime, layerId, wirelocal, side, entrangle );
        dist = dist / 10; // mm->cmo
        if ( layer == -1 || layerId != layer )
        {
          edist = l_mdcCalFunSvc->getSigma( layerId, side, dist, entrangle, tanl,
                                            onWire.z() * 10, h.reccdc()->adc );
        }
        else if ( layerId == layer ) { edist = 99999999; }
        edist     = edist / 10;
        distance  = (double)dist;
        eDistance = (double)edist;
      }
      double eDistance2 = eDistance * eDistance;
 
      //...Residual...
      const double d0 = ( onTrack - onWire ).mag();
      // if(d0>2) std::cout<<"TRF:: strange dist.  d0="<<d0<<" x="<<distance
      //           <<" ex="<<eDistance<<std::endl;
      double dDistance = d0 - distance;
 
      //...dDda...
      for ( int j = 0; j < 5; j++ )
      {
        dDda[j] = ( d[i - 1][j] - d0 ) / ea[j];
        // if(dDda[j]==0) std::cout<<"TRF:: dDda=0 j="<<j<<" ea="<<ea[j]<<std::endl;
      }
      //      for(int j=0;j<5;j++) dDda[j]=(d[i-1][j]-d2[i-1][j])/ea[j]/2.;
      //...Chi2 related...
      dchi2da += ( dDistance / eDistance2 ) * dDda;
      d2chi2d2a += vT_times_v( dDda ) / eDistance2;
      double pChi2 = dDistance * dDistance / eDistance2;
      chi2 += pChi2;
      // if(!(pChi2<0 || pChi2>=0)){
      //    std::cout<<"TRF::  pChi2="<<pChi2<<" X="<<d0<<" fx="<<distance
      //        <<" ex="<<eDistance<<std::endl;
      // }
 
      //...Store results...
      if ( layer == -1 )
      {
        l->update( onTrack, onWire, leftRight, pChi2 );
        l->drift( distance, 0 );
        l->drift( distance, 1 );
        l->dDrift( eDistance, 0 );
        l->dDrift( eDistance, 1 );
      }
 
      else if ( layerId == layer ) { l->distance( ( onTrack - onWire ).mag() ); }
    } // end of loop with hits
 
    //...Check condition...
    double change = chi2Old - chi2;
    if ( 0 <= change && change < 0.01 ) break;
    if ( change < 0. )
    {
      factor *= 100;
      a += da; // recover
      if ( -0.01 < change )
      {
        d2chi2d2a = alpha;
        chi2      = chi2Old;
        break;
      }
    }
    else if ( change > 0. )
    {
      chi2Old = chi2;
      factor *= 0.1;
      alpha = d2chi2d2a;
      beta  = dchi2da;
    }
    else if ( nTrial == 0 )
    {
      chi2Old = chi2;
      factor *= 0.1;
      alpha = d2chi2d2a;
      beta  = dchi2da;
    }
    else
    {
      std::cout << "TrkReco:TRF::  bad chi2 = " << chi2 << std::endl;
      err = TFitFailed;
      //      break;    // protection for nan
      return err;
    }
    alpha2 = alpha;
    for ( unsigned j = 0; j < 5; j++ ) alpha2[j][j] *= ( 1 + factor );
    //...Cal. helix parameters for next loop...
    da = solve( alpha2, beta );
 
    // lclock=clock()/1000;
    // std::cout<<"TRF "<<nTrial<<": "
    //   <<"cl="<<lclock-lclock0<<": "
    //  <<a[0]<<","<<a[1]<<","<<a[2]<<","<<a[3]<<","<<a[4]<<" "
    //  <<chi2<<"/"<<nCores<<":"<<factor
    //  <<" :da "<<da[0]<<","<<da[1]<<","<<da[2]<<","<<da[3]<<","<<da[4]<<std::endl;
    // lclock0=lclock;
 
    a -= da;
    t.a( a );
    // ea = 0.0001*da;
    for ( i = 0; i < 5; i++ )
    {
      ea[i] = 0.0001 * abs( da[i] );
      if ( ea[i] > ea_init ) ea[i] = ea_init;
      if ( ea[i] < 1.0e-10 ) ea[i] = 1.0e-10;
    }
    ++nTrial;
  }
  // std::cout<<"TRungeFitter:: nTrial="<<nTrial<<std::endl;
 
  //...Cal. error matrix...
  SymMatrix Ea( 5, 0 );
  unsigned  dim;
  dim = 5;
  Ea  = d2chi2d2a.inverse( err );
  if ( nCores ) { t.approach( *last, tof, p, _BesBeamPipeMaterials[0], _sag ); }
  // consider the material effect beam pipe.
  double y_[6] = { 0, 0, 0, 0, 0, 0 };
  t.SetFirst( y_ ); // obtain the momentum of the first layer
  int lmass = 0;
  innerwall( t, lmass, y_ ); // consider the material layer by layer
  int nsign = a[2] / abs( a[2] );
  // Update the helix parameter (momentum part.)
  a[4] = y_[5] / sqrt( y_[4] * y_[4] + y_[3] * y_[3] );
  a[2] = nsign * 1 / sqrt( y_[4] * y_[4] + y_[3] * y_[3] );
  //...Store information...
  if ( !err && layer == -1 )
  {
    t.a( a );
    t.Ea( Ea );
    t._fitted = true;
  }
  else if ( !err && layer != -1 )
  {
    t.a( ainitial );
    //    t.Ea(Ea);
    t._fitted = true;
  }
  else { err = TFitFailed; }
 
  t._ndf  = nCores - dim;
  t._chi2 = chi2;
 
  //...Recover pivot...
  t.pivot( pivot_bak );
 
  delete[] d;
  //  delete [] d2;
 
  return err;
}
void TRungeFitter::setBesFromGdml( void ) {
  int    i( 0 );
  double Z( 0. ), A( 0. ), Ionization( 0. ), Density( 0. ), Radlen( 0. );
 
  G4LogicalVolume* logicalMdc = 0;
  MdcG4Geo*        aMdcG4Geo  = new MdcG4Geo();
  logicalMdc                  = aMdcG4Geo->GetTopVolume();
 
  /// mdcgas
  G4Material* mdcMaterial = logicalMdc->GetMaterial();
 
  for ( i = 0; i < mdcMaterial->GetElementVector()->size(); i++ )
  {
    Z += ( mdcMaterial->GetElement( i )->GetZ() ) * ( mdcMaterial->GetFractionVector()[i] );
    A += ( mdcMaterial->GetElement( i )->GetA() ) * ( mdcMaterial->GetFractionVector()[i] );
  }
  Ionization = mdcMaterial->GetIonisation()->GetMeanExcitationEnergy();
  Density    = mdcMaterial->GetDensity() / ( g / cm3 );
  Radlen     = mdcMaterial->GetRadlen();
  RkFitMaterial FitMdcMaterial( Z, A / g / mole, Ionization / eV, Density, Radlen / 10. );
  cout << "mdcgas: Z: " << Z << " A: " << ( A / ( g / mole ) )
       << " Ionization: " << ( Ionization / eV ) << " Density: " << Density
       << " Radlen: " << Radlen << endl;
  _BesBeamPipeMaterials.push_back( FitMdcMaterial );
  // RkFitTrack::mdcGasRadlen_ = Radlen/10.;
 
  /// inner wall aluminium
  G4LogicalVolume* innerwallVolume = const_cast<G4LogicalVolume*>(
      GDMLProcessor::GetInstance()->GetLogicalVolume( "logicalMdcSegment2" ) );
  G4Material* innerwallMaterial = innerwallVolume->GetMaterial();
  G4Tubs*     innerwallTub      = dynamic_cast<G4Tubs*>( innerwallVolume->GetSolid() );
 
  Z = 0.;
  A = 0.;
  for ( i = 0; i < innerwallMaterial->GetElementVector()->size(); i++ )
  {
    Z += ( innerwallMaterial->GetElement( i )->GetZ() ) *
         ( innerwallMaterial->GetFractionVector()[i] );
    A += ( innerwallMaterial->GetElement( i )->GetA() ) *
         ( innerwallMaterial->GetFractionVector()[i] );
  }
 
  Ionization = innerwallMaterial->GetIonisation()->GetMeanExcitationEnergy();
  Density    = innerwallMaterial->GetDensity() / ( g / cm3 );
  Radlen     = innerwallMaterial->GetRadlen();
  cout << "Mdc innerwall, Al: Z: " << Z << " A: " << ( A / ( g / mole ) )
       << " Ionization: " << ( Ionization / eV ) << " Density: " << Density
       << " Radlen: " << Radlen << endl;
  RkFitMaterial FitInnerwallMaterial( Z, A / g / mole, Ionization / eV, Density,
                                      Radlen / 10. );
  _BesBeamPipeMaterials.push_back( FitInnerwallMaterial );
 
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  G4LogicalVolume* logicalBes = 0;
  BesG4Geo*        aBesG4Geo  = new BesG4Geo();
  logicalBes                  = aBesG4Geo->GetTopVolume();
 
  /// air
  G4LogicalVolume* logicalAirVolume = const_cast<G4LogicalVolume*>(
      GDMLProcessor::GetInstance()->GetLogicalVolume( "logicalWorld" ) );
  G4Material* airMaterial = logicalAirVolume->GetMaterial();
  Z                       = 0.;
  A                       = 0.;
  for ( i = 0; i < airMaterial->GetElementVector()->size(); i++ )
  {
    Z += ( airMaterial->GetElement( i )->GetZ() ) * ( airMaterial->GetFractionVector()[i] );
    A += ( airMaterial->GetElement( i )->GetA() ) * ( airMaterial->GetFractionVector()[i] );
  }
  Ionization = airMaterial->GetIonisation()->GetMeanExcitationEnergy();
  Density    = airMaterial->GetDensity() / ( g / cm3 );
  Radlen     = airMaterial->GetRadlen();
  cout << "air: Z: " << Z << " A: " << ( A / ( g / mole ) )
       << " Ionization: " << ( Ionization / eV ) << " Density: " << Density
       << " Radlen: " << Radlen << endl;
  RkFitMaterial FitAirMaterial( Z, A / g / mole, Ionization / eV, Density, Radlen / 10. );
  _BesBeamPipeMaterials.push_back( FitAirMaterial );
 
  /// outer beryllium pipe
  G4LogicalVolume* logicalOuterBeVolume = const_cast<G4LogicalVolume*>(
      GDMLProcessor::GetInstance()->GetLogicalVolume( "logicalouterBe" ) );
  G4Material* outerBeMaterial = logicalOuterBeVolume->GetMaterial();
 
  G4Tubs* outerBeTub = dynamic_cast<G4Tubs*>( logicalOuterBeVolume->GetSolid() );
  Z                  = 0.;
  A                  = 0.;
  for ( i = 0; i < outerBeMaterial->GetElementVector()->size(); i++ )
  {
    Z += ( outerBeMaterial->GetElement( i )->GetZ() ) *
         ( outerBeMaterial->GetFractionVector()[i] );
    A += ( outerBeMaterial->GetElement( i )->GetA() ) *
         ( outerBeMaterial->GetFractionVector()[i] );
  }
  Ionization = outerBeMaterial->GetIonisation()->GetMeanExcitationEnergy();
  Density    = outerBeMaterial->GetDensity() / ( g / cm3 );
  Radlen     = outerBeMaterial->GetRadlen();
  cout << "outer beryllium: Z: " << Z << " A: " << ( A / ( g / mole ) )
       << " Ionization: " << ( Ionization / eV ) << " Density: " << Density
       << " Radlen: " << Radlen << endl;
  RkFitMaterial FitOuterBeMaterial( Z, A / g / mole, Ionization / eV, Density, Radlen / 10. );
  _BesBeamPipeMaterials.push_back( FitOuterBeMaterial );
 
  /// cooling oil
  G4LogicalVolume* logicalOilLayerVolume = const_cast<G4LogicalVolume*>(
      GDMLProcessor::GetInstance()->GetLogicalVolume( "logicaloilLayer" ) );
  G4Material* oilLayerMaterial = logicalOilLayerVolume->GetMaterial();
  G4Tubs*     oilLayerTub      = dynamic_cast<G4Tubs*>( logicalOilLayerVolume->GetSolid() );
 
  Z = 0.;
  A = 0.;
  for ( i = 0; i < oilLayerMaterial->GetElementVector()->size(); i++ )
  {
    Z += ( oilLayerMaterial->GetElement( i )->GetZ() ) *
         ( oilLayerMaterial->GetFractionVector()[i] );
    A += ( oilLayerMaterial->GetElement( i )->GetA() ) *
         ( oilLayerMaterial->GetFractionVector()[i] );
  }
  Ionization = oilLayerMaterial->GetIonisation()->GetMeanExcitationEnergy();
  Density    = oilLayerMaterial->GetDensity() / ( g / cm3 );
  Radlen     = oilLayerMaterial->GetRadlen();
  cout << "cooling oil: Z: " << Z << " A: " << ( A / ( g / mole ) )
       << " Ionization: " << ( Ionization / eV ) << " Density: " << Density
       << " Radlen: " << Radlen << endl;
  RkFitMaterial FitOilLayerMaterial( Z, A / g / mole, Ionization / eV, Density, Radlen / 10. );
  _BesBeamPipeMaterials.push_back( FitOilLayerMaterial );
 
  /// inner beryllium pipe
  G4LogicalVolume* logicalInnerBeVolume = const_cast<G4LogicalVolume*>(
      GDMLProcessor::GetInstance()->GetLogicalVolume( "logicalinnerBe" ) );
 
  G4Material* innerBeMaterial = logicalInnerBeVolume->GetMaterial();
  G4Tubs*     innerBeTub      = dynamic_cast<G4Tubs*>( logicalInnerBeVolume->GetSolid() );
  Z                           = 0.;
  A                           = 0.;
  for ( i = 0; i < innerBeMaterial->GetElementVector()->size(); i++ )
  {
    Z += ( innerBeMaterial->GetElement( i )->GetZ() ) *
         ( innerBeMaterial->GetFractionVector()[i] );
    A += ( innerBeMaterial->GetElement( i )->GetA() ) *
         ( innerBeMaterial->GetFractionVector()[i] );
  }
  Ionization = innerBeMaterial->GetIonisation()->GetMeanExcitationEnergy();
  Density    = innerBeMaterial->GetDensity() / ( g / cm3 );
  Radlen     = innerBeMaterial->GetRadlen();
  cout << "inner beryllium: Z: " << Z << " A: " << ( A / ( g / mole ) )
       << " Ionization: " << ( Ionization / eV ) << " Density: " << Density
       << " Radlen: " << Radlen << endl;
  RkFitMaterial FitInnerBeMaterial( Z, A / g / mole, Ionization / eV, Density, Radlen / 10. );
  _BesBeamPipeMaterials.push_back( FitInnerBeMaterial );
 
  /// gold
  G4LogicalVolume* logicalGoldLayerVolume = const_cast<G4LogicalVolume*>(
      GDMLProcessor::GetInstance()->GetLogicalVolume( "logicalgoldLayer" ) );
  G4Material* goldLayerMaterial = logicalGoldLayerVolume->GetMaterial();
  G4Tubs*     goldLayerTub      = dynamic_cast<G4Tubs*>( logicalGoldLayerVolume->GetSolid() );
 
  Z = 0.;
  A = 0.;
  for ( i = 0; i < goldLayerMaterial->GetElementVector()->size(); i++ )
  {
    Z += ( goldLayerMaterial->GetElement( i )->GetZ() ) *
         ( goldLayerMaterial->GetFractionVector()[i] );
    A += ( goldLayerMaterial->GetElement( i )->GetA() ) *
         ( goldLayerMaterial->GetFractionVector()[i] );
  }
  Ionization = goldLayerMaterial->GetIonisation()->GetMeanExcitationEnergy();
  Density    = goldLayerMaterial->GetDensity() / ( g / cm3 );
  Radlen     = goldLayerMaterial->GetRadlen();
  cout << "gold layer: Z: " << Z << " A: " << ( A / ( g / mole ) )
       << " Ionization: " << ( Ionization / eV ) << " Density: " << Density
       << " Radlen: " << Radlen << endl;
  RkFitMaterial FitGoldLayerMaterial( Z, A / g / mole, Ionization / eV, Density,
                                      Radlen / 10. );
  _BesBeamPipeMaterials.push_back( FitGoldLayerMaterial );
 
  /// now construct the cylinders
  double radius, thick, length, z0;
 
  /// innerwall of inner drift chamber
  radius = innerwallTub->GetInnerRadius() / ( cm );
  thick  = innerwallTub->GetOuterRadius() / (cm)-innerwallTub->GetInnerRadius() / ( cm );
  length = 2.0 * innerwallTub->GetZHalfLength() / ( cm );
  z0     = 0.0;
  cout << "innerwall: "
       << " radius: " << radius << " thick:" << thick << " length: " << length << endl;
  RkFitCylinder innerwallCylinder( &_BesBeamPipeMaterials[1], radius, thick, length, z0 );
  _BesBeamPipeWalls.push_back( innerwallCylinder );
 
  /// outer air, be attention the calculation of the radius and thick of the air cylinder is
  /// special
  radius = outerBeTub->GetOuterRadius() / ( cm );
  thick  = innerwallTub->GetInnerRadius() / (cm)-outerBeTub->GetOuterRadius() / ( cm );
  length = 2.0 * innerwallTub->GetZHalfLength() / ( cm );
  z0     = 0.0;
  cout << "outer air: "
       << " radius: " << radius << " thick:" << thick << " length: " << length << endl;
  RkFitCylinder outerAirCylinder( &_BesBeamPipeMaterials[2], radius, thick, length, z0 );
  _BesBeamPipeWalls.push_back( outerAirCylinder );
 
  /// outer Beryllium layer
  radius = outerBeTub->GetInnerRadius() / ( cm );
  thick  = outerBeTub->GetOuterRadius() / (cm)-outerBeTub->GetInnerRadius() / ( cm );
  length = 2.0 * outerBeTub->GetZHalfLength() / ( cm );
  z0     = 0.0;
  cout << "outer Be: "
       << " radius: " << radius << " thick:" << thick << " length: " << length << endl;
  RkFitCylinder outerBeCylinder( &_BesBeamPipeMaterials[3], radius, thick, length, z0 );
  _BesBeamPipeWalls.push_back( outerBeCylinder );
 
  /// oil layer
  radius = oilLayerTub->GetInnerRadius() / ( cm );
  thick  = oilLayerTub->GetOuterRadius() / (cm)-oilLayerTub->GetInnerRadius() / ( cm );
  length = 2.0 * oilLayerTub->GetZHalfLength() / ( cm );
  z0     = 0.0;
  cout << "oil layer: "
       << " radius: " << radius << " thick:" << thick << " length: " << length << endl;
  RkFitCylinder oilLayerCylinder( &_BesBeamPipeMaterials[4], radius, thick, length, z0 );
  _BesBeamPipeWalls.push_back( oilLayerCylinder );
 
  /// inner Beryllium layer
  radius = innerBeTub->GetInnerRadius() / ( cm );
  thick  = innerBeTub->GetOuterRadius() / (cm)-innerBeTub->GetInnerRadius() / ( cm );
  length = 2.0 * innerBeTub->GetZHalfLength() / ( cm );
  z0     = 0.0;
  cout << "inner Be: "
       << " radius: " << radius << " thick:" << thick << " length: " << length << endl;
  RkFitCylinder innerBeCylinder( &_BesBeamPipeMaterials[5], radius, thick, length, z0 );
  _BesBeamPipeWalls.push_back( innerBeCylinder );
 
  /// gold layer
  radius = goldLayerTub->GetInnerRadius() / ( cm );
  thick  = goldLayerTub->GetOuterRadius() / (cm)-goldLayerTub->GetInnerRadius() / ( cm );
  length = 2.0 * goldLayerTub->GetZHalfLength() / ( cm );
  z0     = 0.0;
  cout << "gold layer: "
       << " radius: " << radius << " thick:" << thick << " length: " << length << endl;
  RkFitCylinder goldLayerCylinder( &_BesBeamPipeMaterials[6], radius, thick, length, z0 );
  _BesBeamPipeWalls.push_back( goldLayerCylinder );
  delete aMdcG4Geo;
  delete aBesG4Geo;
} // end of setBesFromGdml
 
void TRungeFitter::innerwall( TRunge& track, int l_mass, double y[6] ) const {
  for ( int i = 0; i < _BesBeamPipeWalls.size(); i++ )
  { _BesBeamPipeWalls[i].updateTrack( track, y ); }
}
 
IBesMagFieldSvc* TRungeFitter::getPmgnIMF() const {
  if ( nullptr == m_pmgnIMF )
  {
    StatusCode sc = Gaudi::svcLocator()->service( "MagneticFieldSvc", m_pmgnIMF );
    if ( sc.isFailure() )
    { std::cout << "Unable to open Magnetic field service" << std::endl; }
  }
  return m_pmgnIMF;
}