VertexExtrapolate.cxx

Go to the documentation of this file.

/* <===<===<===<===<===<===<===<===<===~===>===>===>===>===>===>===>===>===>
 * File Name:    VertexExtrapolate.cxx
 * Author:       Hao-Kai SUN
 * Created:      2021-09-06 Mon 18:05:16 CST
 * <<=====================================>>
 * Last Updated: 2023-09-29 Thu 14:24:20 CST
 *           By: Hao-Kai SUN
 *     Update #: 213
 * ============================== CODES ==============================>>> */
#include "VertexFitRefine/VertexExtrapolate.h"
 
#include "G4Geo/BesG4Geo.h"
#include "G4Geo/MdcG4Geo.h"
#include "KalFitAlg/KalFitTrack.h"
 
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/ISvcLocator.h"
 
#include "MagneticFieldSvc/IBesMagFieldSvc.h"
 
static const char* matNames[9] = {
    "MDCGas",          "LogicalMdcInnerFilm1", "logicalMdcSegment2", "LogicalMdcInnerFilm0",
    "logicalWorld",    "logicalouterBe",       "logicaloilLayer",    "logicalinnerBe",
    "logicalgoldLayer" };
 
VertexExtrapolate* VertexExtrapolate::m_instance = nullptr;
 
VertexExtrapolate* VertexExtrapolate::instance() {
  if ( m_instance == nullptr ) m_instance = new VertexExtrapolate();
  return m_instance;
}
 
VertexExtrapolate::VertexExtrapolate() {
  m_BesKalmanExtMaterials.clear();
  m_BesKalmanExtWalls.clear();
 
  constructWallsFromGdml();
 
  m_helixVector = CLHEP::Hep3Vector( 5, 0 );
  m_errorMatrix = CLHEP::HepSymMatrix( 5, 0 );
 
  ISvcLocator* svcLocator = Gaudi::svcLocator();
 
  IBesMagFieldSvc* IMFSvc = nullptr;
 
  StatusCode sc = svcLocator->service( "MagneticFieldSvc", IMFSvc );
  if ( sc.isFailure() )
  { std::cerr << MSG::ERROR << "Unable to open Magnetic field service" << std::endl; }
  KalFitTrack::numf_ = 21; // ???
  KalFitElement::muls( 1 );
  KalFitElement::loss( 1 );
  KalFitTrack::setMagneticFieldSvc( IMFSvc );
  KalFitTrack::Bznom_ = -10.; // ???
}
 
void VertexExtrapolate::G4MtovKalFitM( G4Material* g4m, const std::string& name ) {
  double Z = 0.0;
  double A = 0.0;
  for ( unsigned i = 0; i != g4m->GetElementVector()->size(); ++i )
  {
    Z += ( g4m->GetElement( i )->GetZ() ) * ( g4m->GetFractionVector()[i] );
    A += ( g4m->GetElement( i )->GetA() ) * ( g4m->GetFractionVector()[i] );
  }
  double Ionization = g4m->GetIonisation()->GetMeanExcitationEnergy();
  double Density    = g4m->GetDensity() / ( g / cm3 );
  double Radlen     = g4m->GetRadlen();
  std::cout << " " << name << ": Z:  " << Z << " A: " << ( A / ( g / mole ) )
            << " Ionization: " << ( Ionization / eV ) << " Density: " << Density
            << " Radlen:     " << Radlen << std::endl;
  KalFitMaterial FitMdcMaterial( Z, A / g / mole, Ionization / eV, Density, Radlen / 10. );
  m_BesKalmanExtMaterials.push_back( FitMdcMaterial );
}
 
void VertexExtrapolate::AddWalls( int index, double radius, double thick, double length,
                                  double z0 ) {
  std::cout << " " << matNames[index] << ": radius: " << radius << ",  thick: " << thick
            << ",  length: " << length << std::endl;
 
  KalFitMaterial& tmp = m_BesKalmanExtMaterials[index];
  KalFitCylinder  innerwallCylinder( &tmp, radius, thick, length, z0 );
  m_BesKalmanExtWalls.push_back( innerwallCylinder );
}
 
void VertexExtrapolate::AddWalls( int index ) {
  G4Tubs* g4t    = getTubs( matNames[index] );
  double  radius = g4t->GetInnerRadius() / ( cm );
  double  thick  = g4t->GetOuterRadius() / (cm)-g4t->GetInnerRadius() / ( cm );
  double  length = 2.0 * g4t->GetZHalfLength() / ( cm );
  double  z0     = 0.0;
 
  AddWalls( index, radius, thick, length, z0 );
}
 
void VertexExtrapolate::testMW( int index ) {
  KalFitMaterial m = m_BesKalmanExtMaterials[index];
  KalFitCylinder w = m_BesKalmanExtWalls[index];
  std::cout << index << " ==========================" << std::endl;
  std::cout << "TEST==Radlen:  " << m.X0() << std::endl;
  std::cout << "TEST==Radlen:  " << w.material().X0() << "      radius: " << w.radius()
            << std::endl;
}
 
void VertexExtrapolate::constructWallsFromGdml( void ) {
  /// mdcgas
  G4LogicalVolume* logicalMdc = 0;
  MdcG4Geo*        aMdcG4Geo  = new MdcG4Geo();
  logicalMdc                  = aMdcG4Geo->GetTopVolume();
  G4Material* mat             = logicalMdc->GetMaterial();
  KalFitTrack::mdcGasRadlen_  = mat->GetRadlen() / 10.;
  G4MtovKalFitM( mat, "MDCGas" );
 
  BesG4Geo* aBesG4Geo = new BesG4Geo();
  // G4LogicalVolume* logicalBes = aBesG4Geo->GetTopVolume();
 
  for ( int i = 1; i != 9; ++i )
  {
    G4LogicalVolume* logicalAirVolume = const_cast<G4LogicalVolume*>(
        GDMLProcessor::GetInstance()->GetLogicalVolume( matNames[i] ) );
    mat = logicalAirVolume->GetMaterial();
    G4MtovKalFitM( mat, matNames[i] );
  }
  delete aMdcG4Geo;
  delete aBesG4Geo;
  for ( int i = 1; i != 4; ++i ) AddWalls( i );
 
  /// air cylinder is special
  G4Tubs* innerwallTub = getTubs( "LogicalMdcInnerFilm0" );
  G4Tubs* outerBeTub   = getTubs( "logicalouterBe" );
  double  radius       = outerBeTub->GetOuterRadius() / ( cm );
  double  thick  = innerwallTub->GetInnerRadius() / (cm)-outerBeTub->GetOuterRadius() / ( cm );
  double  length = 2.0 * innerwallTub->GetZHalfLength() / ( cm );
  double  z0     = 0.0;
  AddWalls( 4, radius, thick, length, z0 );
 
  for ( int i = 5; i != 9; ++i ) AddWalls( i );
 
  /// the air in the innermost of the pipe
  G4Tubs* goldLayerTub = getTubs( "logicalgoldLayer" );
  radius               = 0.;
  thick                = goldLayerTub->GetInnerRadius() / ( cm );
  length               = 2.0 * goldLayerTub->GetZHalfLength() / ( cm );
  z0                   = 0.0;
  AddWalls( 4, radius, thick, length, z0 );
  // // debug
  // for (int i = 0; i < 9; ++i)
  //     testMW(i);
}
 
int VertexExtrapolate::getWallMdcNumber( const HepPoint3D& point ) const {
  int i = m_BesKalmanExtWalls.size() - 1;
  for ( ; i != -1; --i )
    if ( m_BesKalmanExtWalls[i].isInside2( point ) ) break;
 
  return i;
}
 
void VertexExtrapolate::extToAnyPoint( KalFitTrack& track, const HepPoint3D& point ) {
  const int index = getWallMdcNumber( point );
 
  // outside the inner MDC wall Film1
  if ( index == -1 ) return;
  // in the innermost pipe
  if ( index > 0 )
  {
    for ( int j = 0; j < index; j++ ) m_BesKalmanExtWalls[j].updateTrack( track, 1 );
  }
  int size = m_BesKalmanExtWalls.size() - 1;
  if ( index != size )
  {
    const KalFitMaterial& material = m_BesKalmanExtWalls[index].material();
    HepPoint3D            x;
 
    double path = m_BesKalmanExtWalls[index].intersect( track, x, point );
 
    if ( path > 0 )
    {
      // move pivot
      track.pivot_numf( x );
      // multiple scattering and energy loss
      int index_element( index );
      if ( index_element == 0 ) index_element = 1;
      track.ms( path, material, index_element );
      track.eloss( path, material, index_element );
    }
  }
}
 
void VertexExtrapolate::KalFitExt( const HepPoint3D& point, DstMdcKalTrack* kalTrack,
                                   const int pid ) {
  HepVector    tdsHelix  = kalTrack->getFHelix( pid );
  HepSymMatrix tdsMatrix = kalTrack->getFError( pid );
 
  HepPoint3D IP( 0., 0., 0. );
  // construct a KalFitTrack
  KalFitTrack fitTrack( IP, tdsHelix, tdsMatrix, pid, 0, 0 );
 
  // const double radius = 7.885; // centimeter first MDC wire?
  const double rp = point.perp();
  // const double radius = m_BesKalmanExtWalls[m_BesKalmanExtWalls.size() - 1]
  //                           .radius(); // outer radius
  const double     radius    = m_BesKalmanExtWalls[0].radius(); // outer radius
  const double     dphi      = fitTrack.intersect_cylinder( std::max( rp, radius ) );
  const HepPoint3D lastPivot = fitTrack.x( dphi );
  fitTrack.pivot( lastPivot );
  if ( rp <= radius ) extToAnyPoint( fitTrack, point );
 
  ///////// why pivot to IP? ///////
  // // set the pivot back to IP
  fitTrack.pivot( IP ); // for convention?
  setHelixVector( fitTrack.a() );
  setErrorMatrix( fitTrack.Ea() );
}
/* ===================================================================<<< */
/* ================== VertexExtrapolate.cxx ends here =================== */