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    : [email protected]
//-----------------------------------------------------------------------------
// 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 <float.h>
#include <iostream>
#include "TrkReco/TRungeFitter.h"
#include "TrkReco/TRunge.h"
#include "GaudiKernel/Service.h"                        
#include "GaudiKernel/ISvcLocator.h"                    
#include "GaudiKernel/SvcFactory.h" 
#define HEP_SHORT_NAMES
//#include "panther/panther.h"
//#include MDC_H
#include "MdcTables/MdcTables.h"
#include "CLHEP/Matrix/Vector.h"
#include "CLHEP/Matrix/SymMatrix.h"
#include "CLHEP/Matrix/Matrix.h"
#include "TrkReco/TMLink.h"
#include "TrkReco/TTrackBase.h"
#include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
#include "MdcCalibFunSvc/MdcCalibFunSvc.h"
#include "GaudiKernel/StatusCode.h"
#include "GaudiKernel/IInterface.h"
#include "GaudiKernel/Kernel.h"
#include "GaudiKernel/Service.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/SvcFactory.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/IMessageSvc.h"
#include "MdcTables/MdcTables.h"
#include "CLHEP/Matrix/Vector.h"
#include "CLHEP/Matrix/SymMatrix.h"
#include "CLHEP/Matrix/Matrix.h"
#include "TrkReco/TMLink.h"
#include "TrkReco/TTrackBase.h"
//#include "TrkReco/THistAddItem.h"
#include "GaudiKernel/ISvcLocator.h"
#include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
#include "MdcCalibFunSvc/MdcCalibFunSvc.h"
#include "GaudiKernel/Service.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/SvcFactory.h"
#include "AIDA/IHistogram1D.h"
 
#include "AIDA/IHistogram2D.h"
#include "G4Geo/MdcG4Geo.h"
#include "G4Geo/BesG4Geo.h"
#include "G4Material.hh"
#include "G4Tubs.hh"
 
#include <time.h>
 
using CLHEP::HepVector;
using CLHEP::HepMatrix;
using CLHEP::HepSymMatrix;
/*
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){
 
  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){
 
  StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF); 
  if(scmgn!=StatusCode::SUCCESS) { 
    std::cout<< "Unable to open Magnetic field service"<<std::endl;
  }
}
 
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*(m_pmgnIMF->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);
    }
}