TRungeFitter Class Reference

A class to fit a TTrackBase object to a 3D line. More...

#include <TRungeFitter.h>

Inheritance diagram for TRungeFitter:

Public Member Functions
	TRungeFitter (const std::string &name)
	Constructor.
	TRungeFitter (const std::string &name, bool m_sag, int m_prop, bool m_tof)
virtual	~TRungeFitter ()
	Destructor.
void	dump (const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
	dumps debug information.
virtual int	fit (TTrackBase &) const
virtual int	fit (TTrackBase &, float t0Offset, int layer) const
void	innerwall (TRunge &trk, int l_mass, double y[6]) const
void	sag (bool)
const RkFitMaterial	getMaterial (int i) const
void	propagation (int)
void	tof (bool)
void	setBesFromGdml (void)
Public Member Functions inherited from TMFitter
	TMFitter (const std::string &name)
	Constructor.
virtual	~TMFitter ()
	Destructor.
const std::string &	name (void) const
	returns name.
void	dump (const std::string &message=std::string(""), const std::string &prefix=std::string("")) const

Public Attributes
std::vector< RkFitCylinder >	_BesBeamPipeWalls
std::vector< RkFitMaterial >	_BesBeamPipeMaterials

Additional Inherited Members
Protected Member Functions inherited from TMFitter
void	fitDone (TTrackBase &) const
	sets the fitted flag. (Bad implementation)

Detailed Description

A class to fit a TTrackBase object to a 3D line.

Definition at line 34 of file TRungeFitter.h.

Constructor & Destructor Documentation

TRungeFitter() [1/2]

TRungeFitter::TRungeFitter

(

const std::string &

name

)

Constructor.

Definition at line 90 of file TRungeFitter.cxx.

  : 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() [2/2]

TRungeFitter::TRungeFitter	(	const std::string &	name,
		bool	m_sag,
		int	m_prop,
		bool	m_tof )

Definition at line 99 of file TRungeFitter.cxx.

  : 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::~TRungeFitter ( )

virtual

Destructor.

Definition at line 110 of file TRungeFitter.cxx.

                            {
}

Member Function Documentation

dump()

void TRungeFitter::dump	(	const std::string &	message = std::string(""),
		const std::string &	prefix = std::string("") ) const

dumps debug information.

fit() [1/2]

int TRungeFitter::fit ( TTrackBase & tb ) const

virtual

Implements TMFitter.

Definition at line 122 of file TRungeFitter.cxx.

                                         {
  return fit(tb,0,-1);
}

Referenced by TrkReco::execute(), and fit().

fit() [2/2]

int TRungeFitter::fit ( TTrackBase & tb, float t0Offset, int layer ) const

virtual

Definition at line 126 of file TRungeFitter.cxx.

                                                                    {
// 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;
}

getMaterial()

const RkFitMaterial TRungeFitter::getMaterial

(

int

i

)

const

innerwall()

void TRungeFitter::innerwall	(	TRunge &	trk,
		int	l_mass,
		double	y[6] ) const

Definition at line 673 of file TRungeFitter.cxx.

                                                                      {
    for(int i = 0; i < _BesBeamPipeWalls.size(); i++) {
        _BesBeamPipeWalls[i].updateTrack(track,y);
    }
}

Referenced by fit().

propagation()

void TRungeFitter::propagation

(

int

_in

)

Definition at line 116 of file TRungeFitter.cxx.

                                     {
  _propagation = _in;
}

sag()

void TRungeFitter::sag

(

bool

_in

)

Definition at line 113 of file TRungeFitter.cxx.

                              {
  _sag = _in;
}

setBesFromGdml()

void TRungeFitter::setBesFromGdml

(

void

)

mdcgas

inner wall aluminium

air

outer beryllium pipe

cooling oil

inner beryllium pipe

gold

now construct the cylinders

innerwall of inner drift chamber

outer air, be attention the calculation of the radius and thick of the air cylinder is special

outer Beryllium layer

oil layer

inner Beryllium layer

gold layer

Definition at line 461 of file TRungeFitter.cxx.

                                     {
    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

tof()

void TRungeFitter::tof

(

bool

_in

)

Definition at line 119 of file TRungeFitter.cxx.

                              {
  _tof = _in;
}

Referenced by fit().

Member Data Documentation

_BesBeamPipeMaterials

std::vector<RkFitMaterial> TRungeFitter::_BesBeamPipeMaterials

Definition at line 58 of file TRungeFitter.h.

Referenced by fit(), and setBesFromGdml().

_BesBeamPipeWalls

std::vector<RkFitCylinder> TRungeFitter::_BesBeamPipeWalls

Definition at line 57 of file TRungeFitter.h.

Referenced by innerwall(), and setBesFromGdml().

---

The documentation for this class was generated from the following files:

• TRungeFitter.h
• TRungeFitter.cxx