TRungeFitter Class Reference

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

#include <TRungeFitter.h>

Inheritance diagram for TRungeFitter:

TMFitter List of all members.

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)
const std::stringname (void) const
 returns name.

Public Attributes

std::vector< RkFitCylinder_BesBeamPipeWalls
std::vector< RkFitMaterial_BesBeamPipeMaterials

Protected Member Functions

void fitDone (TTrackBase &) const
 sets the fitted flag. (Bad implementation)

Private Attributes

IMagneticFieldSvcm_pmgnIMF
bool _sag
int _propagation
bool _tof

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::TRungeFitter ( const std::string name  ) 

Constructor.

Definition at line 90 of file TRungeFitter.cxx.

References m_pmgnIMF.

00091   : TMFitter(name),
00092     _sag(true),_propagation(1),_tof(false){
00093 
00094   StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF); 
00095   if(scmgn!=StatusCode::SUCCESS) { 
00096     std::cout<< __FILE__<<" Unable to open Magnetic field service"<<std::endl;
00097   }
00098 }

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

Definition at line 99 of file TRungeFitter.cxx.

References m_pmgnIMF.

00101   : TMFitter(name),
00102     _sag(m_sag),_propagation(m_prop),_tof(m_tof){
00103 
00104   StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF); 
00105   if(scmgn!=StatusCode::SUCCESS) { 
00106     std::cout<< "Unable to open Magnetic field service"<<std::endl;
00107   }
00108 }

TRungeFitter::~TRungeFitter (  )  [virtual]

Destructor.

Definition at line 110 of file TRungeFitter.cxx.

00110                             {
00111 }


Member Function Documentation

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

dumps debug information.

Reimplemented from TMFitter.

int TRungeFitter::fit ( TTrackBase ,
float  t0Offset,
int  layer 
) const [virtual]

Definition at line 126 of file TRungeFitter.cxx.

References _BesBeamPipeMaterials, _propagation, _sag, _tof, abs, MdcRec_wirhit::adc, alpha, cos(), DBL_MAX, TMDCWireHit::dDrift(), check_raw_filter::dist, IMdcCalibFunSvc::driftTimeToDist(), showlog::err, IMagneticFieldSvc::getReferField(), IMdcCalibFunSvc::getSigma(), IMdcCalibFunSvc::getT0(), IMdcCalibFunSvc::getTimeWalk(), IMdcCalibFunSvc::getTprop(), genRecEmupikp::i, TMDCWire::id(), innerwall(), ganga-rec::j, TMDCWire::layerId(), TMDCWire::localId(), m_pmgnIMF, NStereoHits(), TTrackBase::objectType(), ORIGIN, phi0, pi, TMDCWireHit::reccdc(), Runge, sin(), t(), MdcRec_wirhit::tdc, TFitAlreadyFitted, TFitFailed, TFitUnavailable, tof(), TMDCWireHit::wire(), WireHitLeft, WireHitRight, and x.

00126                                                                     {
00127 // Argument layer is used for calibration interface in which doca is calculated without measured hit. liucy
00128   IMdcCalibFunSvc* l_mdcCalFunSvc;
00129   Gaudi::svcLocator() -> service("MdcCalibFunSvc", l_mdcCalFunSvc);
00130   //std::cout<<"TRungeFitter::fit  start"<<std::endl;
00131   //...Type check...
00132   if(tb.objectType() != Runge) return TFitUnavailable;
00133   TRunge& t = (TRunge&) tb;
00134   //...Already fitted ?...
00135   if(t.fitted()&&layer==-1) return TFitAlreadyFitted;
00136   //...Count # of hits...
00137   const AList<TMLink>& cores = t.cores();
00138   unsigned nCores = cores.length();
00139   unsigned nStereoCores = NStereoHits(cores);
00140   TMLink* last=NULL;
00141   int lyr=0;
00142   int layerid=0;
00143   for(int i=0;i<nCores;i++){
00144     layerid=(*cores[i]->hit()).wire()->layerId(); 
00145     if(layerid>=lyr){
00146         lyr=layerid;
00147         last=cores[i];
00148     }
00149   }
00150   
00151   //...Check # of hits...
00152  // if ((nStereoCores < 2) || (nCores - nStereoCores < 3))
00153  //   return TFitErrorFewHits;
00154 
00155   //...Move pivot to ORIGIN...
00156   const HepPoint3D pivot_bak = t.pivot();
00157   t.pivot(ORIGIN);
00158 
00159   //...Setup...
00160   Vector a(5),da(5);
00161   a = t.a();
00162   Vector ainitial(a);
00163   Vector dDda(5);
00164   Vector dchi2da(5);
00165   SymMatrix d2chi2d2a(5,0);
00166   const SymMatrix zero5(5,0);
00167   double chi2;
00168   double chi2Old = 0;
00169   for(int i=0;i<t.links().length();i++){
00170     chi2Old=chi2Old+t.links()[i]->pull();
00171   }
00172   chi2Old=DBL_MAX;
00173   int err = 0;
00174   
00175   double factor = 0.1;
00176   Vector beta(5);
00177   SymMatrix alpha(5,0);
00178   SymMatrix alpha2(5,0);
00179 
00180   double (*d)[5]=new double[nCores][5];
00181   //  double (*d2)[5]=new double[nCores][5];
00182   Vector ea(5);
00183 
00184   float tof;
00185   HepVector3D p;
00186   unsigned i;
00187 
00188   double distance;
00189   double eDistance;
00190 
00191   // ea... init
00192   const double ea_init=0.000001;
00193   ea[0]=ea_init;                //dr
00194   ea[1]=ea_init;                //phi0
00195   ea[2]=ea_init;                //kappa
00196   ea[3]=ea_init;                //dz
00197   ea[4]=ea_init;                //tanl
00198 
00199   //std::cout<<"TRF ::"<<a[0]<<","<<a[1]<<","<<a[2]<<","<<a[3]<<","<<a[4]<<std::endl;
00200 
00201   //long int lclock0=clock()/1000;
00202   //long int lclock=lclock0;
00203 
00204   Vector def_a(t.a());
00205   Vector ta(def_a);
00206 
00207   //...Fitting loop...
00208   unsigned nTrial = 0;
00209   while(nTrial < 100){
00210     
00211     //...Set up...
00212     chi2 = 0;
00213     for (unsigned j=0;j<5;j++) dchi2da[j]=0;
00214     d2chi2d2a=zero5;
00215 
00216     def_a=t.a();
00217 
00218     //#### loop for shifted helix parameter set ####
00219     for(unsigned j=0;j<5;j++){
00220       ta=def_a;
00221       ta[j]+=ea[j];
00222       t.a(ta);
00223       //...Loop with hits...
00224       i=0;
00225       //std::cout<<"TRF:: cores="<<cores.length()<<std::endl;
00226       while(TMLink * l = cores[i++]){
00227         //...Cal. closest points...
00228         t.approach(* l,tof,p,_BesBeamPipeMaterials[0],_sag);
00229         const HepPoint3D& onTrack=l->positionOnTrack();
00230         const HepPoint3D& onWire=l->positionOnWire();
00231         d[i-1][j]=(onTrack-onWire).mag();
00232         //std::cout<<"TRF:: i="<<i<<std::endl;
00233       }//end of loop with hits
00234       //lclock=clock()/1000;
00235       //std::cout<<"TRF  clock="<<lclock-lclock0<<std::endl;
00236       //lclock0=lclock;
00237     }
00238     /*
00239     for(int j=0;j<5;j++){
00240       ta=def_a;
00241       ta[j]-=ea[j];
00242       t.a(ta);
00243       //...Loop with hits...
00244       float tof_dummy;
00245       Vector3 p_dummy;
00246       unsigned i=0;
00247       while(TMLink* l = cores[i++]){
00248         //...Cal. closest points...
00249         t.approach(*l,tof_dummy,p_dummy,_sag);
00250         const HepPoint3D& onTrack=l->positionOnTrack();
00251         const HepPoint3D& onWire=l->positionOnWire();
00252         d2[i-1][j]=(onTrack-onWire).mag();
00253       }//end of loop with hits
00254     }
00255     */
00256     t.a(def_a);
00257       
00258     //#### original parameter set  and  calc. chi2 ####
00259     //...Loop with hits...
00260     i=0;
00261     while(TMLink * l = cores[i++]){
00262       const TMDCWireHit& h = * l->hit();
00263       //...Cal. closest points...
00264       if(t.approach(* l,tof,p,_BesBeamPipeMaterials[0],_sag)<0){
00265         std::cout<<"TrkReco:TRF::  bad wire"<<std::endl;
00266         continue;
00267       }
00268       int layerId=h.wire()->layerId();
00269       const HepPoint3D& onTrack=l->positionOnTrack();
00270       const HepPoint3D& onWire=l->positionOnWire();
00271       unsigned leftRight = WireHitRight;
00272       if (onWire.cross(onTrack).z() < 0.) leftRight = WireHitLeft;
00273 
00274       //...Obtain drift distance and its error...
00275       if ((t0Offset == 0.) && (_propagation==0) && (! _tof)) {
00276         //...No correction...
00277         distance = l->drift(leftRight);
00278         eDistance = h.dDrift(leftRight);
00279       }else{
00280         //...T0 and propagation corrections...
00281         int wire = h.wire()->id();
00282         int wirelocal=h.wire()->localId();
00283         int side = leftRight;
00284         if (side==0) side = 0;
00285         float tp[3] = {p.x(),p.y(),p.z()};
00286         float x[3] = {onWire.x(), onWire.y(), onWire.z()};
00287         double tprop = l_mdcCalFunSvc->getTprop(layerId,onWire.z()*10.);
00288         double timeWalk = l_mdcCalFunSvc->getTimeWalk(layerId, h.reccdc()->adc);
00289         double T0 = l_mdcCalFunSvc->getT0(layerId,wirelocal);
00290         double drifttime = h.reccdc()->tdc - tof - tprop - timeWalk -T0;
00291         l->setDriftTime(drifttime);
00292         float dist;
00293         float edist;
00294         int prop = _propagation;
00295         const double x0 = t.helix().pivot().x();
00296         const double y0 = t.helix().pivot().y();
00297         Hep3Vector pivot_helix(x0,y0,0);
00298         const double dr    =  t.helix().a()[0];
00299         const double phi0  =  t.helix().a()[1];
00300         const double kappa =  t.helix().a()[2];
00301         const double dz    =  t.helix().a()[3];
00302         const double tanl  =  t.helix().a()[4];
00303 
00304         const double Bz = -1000*(m_pmgnIMF->getReferField());
00305         const double alpha = 333.564095 / Bz;
00306 
00307         const double cox = x0 + dr*cos(phi0)- alpha*cos(phi0)/kappa;
00308         const double coy = y0 + dr*sin(phi0) - alpha*sin(phi0)/kappa;
00309         HepPoint3D dir(onTrack.y()-coy,cox-onTrack.x(),0);
00310         double pos_phi=onWire.phi();
00311         double dir_phi=dir.phi();
00312         while(pos_phi>pi){pos_phi-=pi;}
00313         while(pos_phi<0){pos_phi+=pi;}
00314         while(dir_phi>pi){dir_phi-=pi;}
00315         while(dir_phi<0){dir_phi+=pi;}
00316         double entrangle=dir_phi-pos_phi;
00317         while(entrangle>pi/2)entrangle-=pi;
00318         while(entrangle<(-1)*pi/2)entrangle+=pi;
00319         dist = l_mdcCalFunSvc->driftTimeToDist(drifttime,layerId, wirelocal, side, entrangle);
00320         dist= dist/10;//mm->cmo
00321         if(layer==-1||layerId!=layer){
00322         edist = l_mdcCalFunSvc->getSigma(layerId, side, dist, entrangle,tanl,onWire.z()*10,h.reccdc()->adc);   }
00323         else if(layerId==layer){edist=99999999;}
00324         edist = edist/10;
00325         distance = (double) dist;
00326         eDistance = (double) edist;
00327         
00328       }
00329       double eDistance2=eDistance*eDistance;
00330 
00331       //...Residual...
00332       const double d0 = (onTrack-onWire).mag();
00333       //if(d0>2) std::cout<<"TRF:: strange dist.  d0="<<d0<<" x="<<distance
00334       //                   <<" ex="<<eDistance<<std::endl;
00335       double dDistance = d0 - distance;
00336 
00337       //...dDda...
00338       for(int j=0;j<5;j++){
00339         dDda[j]=(d[i-1][j]-d0)/ea[j];
00340         //if(dDda[j]==0) std::cout<<"TRF:: dDda=0 j="<<j<<" ea="<<ea[j]<<std::endl;
00341       }
00342       //      for(int j=0;j<5;j++) dDda[j]=(d[i-1][j]-d2[i-1][j])/ea[j]/2.;
00343       //...Chi2 related...
00344       dchi2da += (dDistance / eDistance2) * dDda;
00345       d2chi2d2a += vT_times_v(dDda) / eDistance2;
00346       double pChi2 = dDistance * dDistance / eDistance2;
00347       chi2 += pChi2;
00348       //if(!(pChi2<0 || pChi2>=0)){
00349       //        std::cout<<"TRF::  pChi2="<<pChi2<<" X="<<d0<<" fx="<<distance
00350       //            <<" ex="<<eDistance<<std::endl;
00351       //}
00352 
00353       //...Store results...
00354       if(layer==-1){
00355       l->update(onTrack, onWire, leftRight, pChi2);
00356       l->drift(distance,0);
00357       l->drift(distance,1);
00358       l->dDrift(eDistance,0);
00359       l->dDrift(eDistance,1);
00360       }
00361 
00362       else if(layerId==layer){
00363         l->distance((onTrack-onWire).mag());
00364       }
00365     }//end of loop with hits
00366 
00367     //...Check condition...
00368     double change = chi2Old - chi2;
00369     if (0 <= change && change < 0.01) break;
00370     if (change < 0.) {
00371       factor *= 100;
00372       a += da;  //recover
00373       if(-0.01 < change){
00374         d2chi2d2a=alpha;
00375         chi2=chi2Old;
00376         break;
00377       }
00378     }else if(change > 0.){
00379       chi2Old = chi2;
00380       factor *= 0.1;
00381       alpha=d2chi2d2a;
00382       beta=dchi2da;
00383     }else if(nTrial==0){
00384     chi2Old = chi2;
00385           factor *= 0.1;
00386                 alpha=d2chi2d2a;
00387                       beta=dchi2da;
00388     }else{
00389       std::cout<<"TrkReco:TRF::  bad chi2 = "<<chi2<<std::endl;
00390       err=TFitFailed;
00391 //      break;  // protection for nan
00392     return err;
00393     }
00394     alpha2=alpha;
00395     for(unsigned j=0;j<5;j++) alpha2[j][j]*=(1+factor);
00396     //...Cal. helix parameters for next loop...
00397     da = solve(alpha2,beta);
00398 
00399     //lclock=clock()/1000;
00400     //std::cout<<"TRF "<<nTrial<<": "
00401     //  <<"cl="<<lclock-lclock0<<": "
00402     //  <<a[0]<<","<<a[1]<<","<<a[2]<<","<<a[3]<<","<<a[4]<<" "
00403     //  <<chi2<<"/"<<nCores<<":"<<factor
00404     //  <<" :da "<<da[0]<<","<<da[1]<<","<<da[2]<<","<<da[3]<<","<<da[4]<<std::endl;
00405     //lclock0=lclock;
00406 
00407     a -= da;
00408     t.a(a);
00409     //ea = 0.0001*da;
00410     for(i=0;i<5;i++){
00411       ea[i]=0.0001*abs(da[i]);
00412       if(ea[i]>ea_init) ea[i]=ea_init;
00413       if(ea[i]<1.0e-10) ea[i]=1.0e-10;
00414     }
00415     ++nTrial;
00416   }
00417   //std::cout<<"TRungeFitter:: nTrial="<<nTrial<<std::endl;
00418   
00419   //...Cal. error matrix...
00420   SymMatrix Ea(5,0);
00421   unsigned dim;
00422   dim=5;
00423   Ea = d2chi2d2a.inverse(err);
00424   if(nCores){
00425   t.approach(*last,tof,p,_BesBeamPipeMaterials[0],_sag);}
00426   // consider the material effect beam pipe.
00427   double y_[6]={0,0,0,0,0,0};
00428   t.SetFirst(y_);//obtain the momentum of the first layer
00429   int lmass=0;
00430   innerwall(t,lmass,y_);// consider the material layer by layer
00431   int nsign=a[2]/abs(a[2]);
00432   // Update the helix parameter (momentum part.)
00433   a[4]=y_[5]/sqrt(y_[4]*y_[4]+y_[3]*y_[3]);
00434   a[2]=nsign*1/sqrt(y_[4]*y_[4]+y_[3]*y_[3]);
00435   //...Store information...
00436   if(! err&&layer==-1){
00437     t.a(a);
00438     t.Ea(Ea);
00439     t._fitted = true;
00440   }else if(! err&&layer!=-1){
00441     t.a(ainitial);
00442 //    t.Ea(Ea);
00443     t._fitted = true;
00444   }else{
00445     err = TFitFailed;
00446   }
00447 
00448   t._ndf = nCores - dim;
00449   t._chi2 = chi2;
00450 
00451   //...Recover pivot...
00452   t.pivot(pivot_bak);
00453 
00454   delete [] d;
00455   //  delete [] d2;
00456   
00457   return err;
00458 }

int TRungeFitter::fit ( TTrackBase  )  const [virtual]

Implements TMFitter.

Definition at line 122 of file TRungeFitter.cxx.

Referenced by TrkReco::execute().

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

void TMFitter::fitDone ( TTrackBase  )  const [protected, inherited]

sets the fitted flag. (Bad implementation)

Definition at line 24 of file TMFitter.cxx.

References t().

Referenced by TRobustLineFitter::fit(), TLineFitter::fit(), and TCircleFitter::fit().

00024                                       {
00025     t._fitted = true;
00026 }

const RkFitMaterial TRungeFitter::getMaterial ( int  i  )  const

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

Definition at line 671 of file TRungeFitter.cxx.

References _BesBeamPipeWalls, and genRecEmupikp::i.

Referenced by fit().

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

const std::string & TMFitter::name ( void   )  const [inline, inherited]

returns name.

Definition at line 73 of file TMFitter.h.

References TMFitter::_name.

00073                          {
00074     return _name;
00075 }

void TRungeFitter::propagation ( int   ) 

Definition at line 116 of file TRungeFitter.cxx.

References _propagation.

00116                                      {
00117   _propagation = _in;
00118 }

void TRungeFitter::sag ( bool   ) 

Definition at line 113 of file TRungeFitter.cxx.

References _sag.

00113                               {
00114   _sag = _in;
00115 }

void TRungeFitter::setBesFromGdml ( void   ) 

Definition at line 459 of file TRungeFitter.cxx.

References _BesBeamPipeMaterials, _BesBeamPipeWalls, EvtCyclic3::A, SubDetectorG4Geo::GetTopVolume(), and genRecEmupikp::i.

Referenced by TrkReco::initPara().

00459                                      {
00460     int i(0);
00461     double Z(0.),A(0.),Ionization(0.),Density(0.),Radlen(0.);
00462 
00463     G4LogicalVolume *logicalMdc = 0;
00464     MdcG4Geo* aMdcG4Geo = new MdcG4Geo();
00465     logicalMdc = aMdcG4Geo->GetTopVolume();
00466 
00468     G4Material* mdcMaterial = logicalMdc->GetMaterial();
00469 
00470     for(i=0; i<mdcMaterial->GetElementVector()->size(); i++){
00471         Z += (mdcMaterial->GetElement(i)->GetZ())*
00472             (mdcMaterial->GetFractionVector()[i]);
00473         A += (mdcMaterial->GetElement(i)->GetA())*
00474             (mdcMaterial->GetFractionVector()[i]);
00475     }
00476     Ionization = mdcMaterial->GetIonisation()->GetMeanExcitationEnergy();
00477     Density = mdcMaterial->GetDensity()/(g/cm3);
00478     Radlen = mdcMaterial->GetRadlen();
00479     RkFitMaterial FitMdcMaterial(Z,A/g/mole,Ionization/eV,Density,Radlen/10.);
00480     cout<<"mdcgas: Z: "<<Z<<" A: "<<(A/(g/mole))<<" Ionization: "<<(Ionization/eV)<<" Density: "<<Density<<" Radlen: "<<Radlen<<endl;
00481     _BesBeamPipeMaterials.push_back(FitMdcMaterial);
00482      //RkFitTrack::mdcGasRadlen_ = Radlen/10.;
00483 
00485     G4LogicalVolume* innerwallVolume = const_cast<G4LogicalVolume*>(GDMLProcessor::GetInstance()->GetLogicalVolume("logicalMdcSegment2"));
00486     G4Material* innerwallMaterial = innerwallVolume->GetMaterial();
00487     G4Tubs* innerwallTub = dynamic_cast<G4Tubs*>(innerwallVolume->GetSolid());
00488 
00489     Z = 0.;
00490     A = 0.;
00491     for(i=0; i<innerwallMaterial->GetElementVector()->size(); i++){
00492         Z += (innerwallMaterial->GetElement(i)->GetZ())*
00493             (innerwallMaterial->GetFractionVector()[i]);
00494         A += (innerwallMaterial->GetElement(i)->GetA())*
00495             (innerwallMaterial->GetFractionVector()[i]);
00496     }
00497 
00498     Ionization = innerwallMaterial->GetIonisation()->GetMeanExcitationEnergy();
00499     Density = innerwallMaterial->GetDensity()/(g/cm3);
00500     Radlen = innerwallMaterial->GetRadlen();
00501     cout<<"Mdc innerwall, Al: Z: "<<Z<<" A: "<<(A/(g/mole))<<" Ionization: "<<(Ionization/eV)<<" Density: "<<Density<<" Radlen: "<<Radlen<<endl;
00502     RkFitMaterial FitInnerwallMaterial(Z,A/g/mole,Ionization/eV,Density,Radlen/10.);
00503     _BesBeamPipeMaterials.push_back(FitInnerwallMaterial);
00504 
00506     G4LogicalVolume *logicalBes = 0;
00507     BesG4Geo* aBesG4Geo = new BesG4Geo();
00508     logicalBes = aBesG4Geo->GetTopVolume();
00509 
00511     G4LogicalVolume* logicalAirVolume = const_cast<G4LogicalVolume*>(GDMLProcessor::GetInstance()->GetLogicalVolume("logicalWorld"));
00512     G4Material* airMaterial = logicalAirVolume->GetMaterial();
00513     Z = 0.;
00514     A = 0.;
00515     for(i=0; i<airMaterial->GetElementVector()->size(); i++){
00516         Z += (airMaterial->GetElement(i)->GetZ())*
00517             (airMaterial->GetFractionVector()[i]);
00518         A += (airMaterial->GetElement(i)->GetA())*
00519             (airMaterial->GetFractionVector()[i]);
00520     }
00521     Ionization = airMaterial->GetIonisation()->GetMeanExcitationEnergy();
00522     Density = airMaterial->GetDensity()/(g/cm3);
00523     Radlen = airMaterial->GetRadlen();
00524     cout<<"air: Z: "<<Z<<" A: "<<(A/(g/mole))<<" Ionization: "<<(Ionization/eV)<<" Density: "<<Density<<" Radlen: "<<Radlen<<endl;
00525     RkFitMaterial FitAirMaterial(Z,A/g/mole,Ionization/eV,Density,Radlen/10.);
00526     _BesBeamPipeMaterials.push_back(FitAirMaterial);
00527 
00529     G4LogicalVolume* logicalOuterBeVolume = const_cast<G4LogicalVolume*>(GDMLProcessor::GetInstance()->GetLogicalVolume("logicalouterBe"));
00530     G4Material* outerBeMaterial = logicalOuterBeVolume->GetMaterial();
00531 
00532     G4Tubs* outerBeTub = dynamic_cast<G4Tubs*>(logicalOuterBeVolume->GetSolid());
00533     Z = 0.;
00534     A = 0.;
00535     for(i=0; i<outerBeMaterial->GetElementVector()->size(); i++){
00536         Z += (outerBeMaterial->GetElement(i)->GetZ())*
00537             (outerBeMaterial->GetFractionVector()[i]);
00538         A += (outerBeMaterial->GetElement(i)->GetA())*
00539             (outerBeMaterial->GetFractionVector()[i]);
00540     }
00541     Ionization =  outerBeMaterial->GetIonisation()->GetMeanExcitationEnergy();
00542     Density = outerBeMaterial->GetDensity()/(g/cm3);
00543     Radlen = outerBeMaterial->GetRadlen();
00544     cout<<"outer beryllium: Z: "<<Z<<" A: "<<(A/(g/mole))<<" Ionization: "<<(Ionization/eV)<<" Density: "<<Density<<" Radlen: "<<Radlen<<endl;
00545     RkFitMaterial FitOuterBeMaterial(Z,A/g/mole,Ionization/eV,Density,Radlen/10.);
00546     _BesBeamPipeMaterials.push_back(FitOuterBeMaterial);
00547 
00549     G4LogicalVolume* logicalOilLayerVolume = const_cast<G4LogicalVolume*>(GDMLProcessor::GetInstance()->GetLogicalVolume("logicaloilLayer"));
00550     G4Material* oilLayerMaterial = logicalOilLayerVolume->GetMaterial();
00551     G4Tubs* oilLayerTub = dynamic_cast<G4Tubs*>(logicalOilLayerVolume->GetSolid());
00552 
00553     Z = 0.;
00554     A = 0.;
00555     for(i=0; i<oilLayerMaterial->GetElementVector()->size(); i++){
00556         Z += (oilLayerMaterial->GetElement(i)->GetZ())*
00557             (oilLayerMaterial->GetFractionVector()[i]);
00558         A += (oilLayerMaterial->GetElement(i)->GetA())*
00559             (oilLayerMaterial->GetFractionVector()[i]);
00560     }
00561     Ionization = oilLayerMaterial->GetIonisation()->GetMeanExcitationEnergy();
00562     Density = oilLayerMaterial->GetDensity()/(g/cm3);
00563     Radlen = oilLayerMaterial->GetRadlen();
00564     cout<<"cooling oil: Z: "<<Z<<" A: "<<(A/(g/mole))<<" Ionization: "<<(Ionization/eV)<<" Density: "<<Density<<" Radlen: "<<Radlen<<endl;
00565     RkFitMaterial FitOilLayerMaterial(Z,A/g/mole,Ionization/eV,Density,Radlen/10.);
00566     _BesBeamPipeMaterials.push_back(FitOilLayerMaterial);
00567 
00568 
00570     G4LogicalVolume* logicalInnerBeVolume = const_cast<G4LogicalVolume*>(GDMLProcessor::GetInstance()->GetLogicalVolume("logicalinnerBe"));
00571 
00572     G4Material* innerBeMaterial = logicalInnerBeVolume->GetMaterial();
00573     G4Tubs* innerBeTub = dynamic_cast<G4Tubs*>(logicalInnerBeVolume->GetSolid());
00574     Z = 0.;
00575     A = 0.;
00576     for(i=0; i<innerBeMaterial->GetElementVector()->size(); i++){
00577         Z += (innerBeMaterial->GetElement(i)->GetZ())*
00578             (innerBeMaterial->GetFractionVector()[i]);
00579         A += (innerBeMaterial->GetElement(i)->GetA())*
00580             (innerBeMaterial->GetFractionVector()[i]);
00581     }
00582     Ionization = innerBeMaterial->GetIonisation()->GetMeanExcitationEnergy();
00583     Density = innerBeMaterial->GetDensity()/(g/cm3);
00584     Radlen = innerBeMaterial->GetRadlen();
00585     cout<<"inner beryllium: Z: "<<Z<<" A: "<<(A/(g/mole))<<" Ionization: "<<(Ionization/eV)<<" Density: "<<Density<<" Radlen: "<<Radlen<<endl;
00586     RkFitMaterial FitInnerBeMaterial(Z,A/g/mole,Ionization/eV,Density,Radlen/10.);
00587     _BesBeamPipeMaterials.push_back(FitInnerBeMaterial);
00588 
00589 
00591     G4LogicalVolume* logicalGoldLayerVolume = const_cast<G4LogicalVolume*>(GDMLProcessor::GetInstance()->GetLogicalVolume("logicalgoldLayer"));
00592     G4Material* goldLayerMaterial = logicalGoldLayerVolume->GetMaterial();
00593     G4Tubs* goldLayerTub = dynamic_cast<G4Tubs*>(logicalGoldLayerVolume->GetSolid());
00594 
00595     Z = 0.;
00596     A = 0.;
00597     for(i=0; i<goldLayerMaterial->GetElementVector()->size(); i++){
00598         Z += (goldLayerMaterial->GetElement(i)->GetZ())*
00599             (goldLayerMaterial->GetFractionVector()[i]);
00600         A += (goldLayerMaterial->GetElement(i)->GetA())*
00601             (goldLayerMaterial->GetFractionVector()[i]);
00602     }
00603     Ionization = goldLayerMaterial->GetIonisation()->GetMeanExcitationEnergy();
00604     Density = goldLayerMaterial->GetDensity()/(g/cm3);
00605     Radlen = goldLayerMaterial->GetRadlen();
00606     cout<<"gold layer: Z: "<<Z<<" A: "<<(A/(g/mole))<<" Ionization: "<<(Ionization/eV)<<" Density: "<<Density<<" Radlen: "<<Radlen<<endl;
00607     RkFitMaterial FitGoldLayerMaterial(Z,A/g/mole,Ionization/eV,Density,Radlen/10.);
00608     _BesBeamPipeMaterials.push_back(FitGoldLayerMaterial);
00609 
00610 
00612     double radius, thick, length , z0;
00613 
00615     radius = innerwallTub->GetInnerRadius()/(cm);
00616     thick  = innerwallTub->GetOuterRadius()/(cm) - innerwallTub->GetInnerRadius()/(cm);
00617     length = 2.0*innerwallTub->GetZHalfLength()/(cm);
00618     z0     = 0.0;
00619     cout<<"innerwall: "<<" radius: "<<radius<<" thick:"<<thick<<" length: "<<length<<endl;
00620     RkFitCylinder innerwallCylinder(&_BesBeamPipeMaterials[1], radius, thick, length , z0);
00621     _BesBeamPipeWalls.push_back(innerwallCylinder);
00622 
00624     radius = outerBeTub->GetOuterRadius()/(cm);
00625     thick  = innerwallTub->GetInnerRadius()/(cm) - outerBeTub->GetOuterRadius()/(cm);
00626     length = 2.0*innerwallTub->GetZHalfLength()/(cm);
00627     z0     = 0.0;
00628     cout<<"outer air: "<<" radius: "<<radius<<" thick:"<<thick<<" length: "<<length<<endl;
00629     RkFitCylinder outerAirCylinder(&_BesBeamPipeMaterials[2], radius, thick, length , z0);
00630     _BesBeamPipeWalls.push_back(outerAirCylinder);
00631 
00633     radius = outerBeTub->GetInnerRadius()/(cm);
00634     thick  = outerBeTub->GetOuterRadius()/(cm) - outerBeTub->GetInnerRadius()/(cm);
00635     length = 2.0*outerBeTub->GetZHalfLength()/(cm);
00636     z0     = 0.0;
00637     cout<<"outer Be: "<<" radius: "<<radius<<" thick:"<<thick<<" length: "<<length<<endl;
00638     RkFitCylinder outerBeCylinder(&_BesBeamPipeMaterials[3], radius, thick, length , z0);
00639     _BesBeamPipeWalls.push_back(outerBeCylinder);
00640 
00642     radius = oilLayerTub->GetInnerRadius()/(cm);
00643     thick  = oilLayerTub->GetOuterRadius()/(cm) - oilLayerTub->GetInnerRadius()/(cm);
00644     length = 2.0*oilLayerTub->GetZHalfLength()/(cm);
00645     z0     = 0.0;
00646     cout<<"oil layer: "<<" radius: "<<radius<<" thick:"<<thick<<" length: "<<length<<endl;
00647     RkFitCylinder oilLayerCylinder(&_BesBeamPipeMaterials[4], radius, thick, length , z0);
00648     _BesBeamPipeWalls.push_back(oilLayerCylinder);
00649 
00651     radius = innerBeTub->GetInnerRadius()/(cm);
00652     thick  = innerBeTub->GetOuterRadius()/(cm) - innerBeTub->GetInnerRadius()/(cm);
00653     length = 2.0*innerBeTub->GetZHalfLength()/(cm);
00654     z0     = 0.0;
00655     cout<<"inner Be: "<<" radius: "<<radius<<" thick:"<<thick<<" length: "<<length<<endl;
00656     RkFitCylinder innerBeCylinder(&_BesBeamPipeMaterials[5], radius, thick, length , z0);
00657     _BesBeamPipeWalls.push_back(innerBeCylinder);
00658 
00660     radius = goldLayerTub->GetInnerRadius()/(cm);
00661     thick  = goldLayerTub->GetOuterRadius()/(cm) - goldLayerTub->GetInnerRadius()/(cm);
00662     length = 2.0*goldLayerTub->GetZHalfLength()/(cm);
00663     z0     = 0.0;
00664     cout<<"gold layer: "<<" radius: "<<radius<<" thick:"<<thick<<" length: "<<length<<endl;
00665     RkFitCylinder goldLayerCylinder(&_BesBeamPipeMaterials[6], radius, thick, length , z0);
00666     _BesBeamPipeWalls.push_back(goldLayerCylinder);
00667     delete aMdcG4Geo; 
00668     delete aBesG4Geo;
00669 }//end of setBesFromGdml

void TRungeFitter::tof ( bool   ) 

Definition at line 119 of file TRungeFitter.cxx.

References _tof.

Referenced by fit().

00119                               {
00120   _tof = _in;
00121 }


Member Data Documentation

std::vector<RkFitMaterial> TRungeFitter::_BesBeamPipeMaterials

Definition at line 58 of file TRungeFitter.h.

Referenced by fit(), and setBesFromGdml().

std::vector<RkFitCylinder> TRungeFitter::_BesBeamPipeWalls

Definition at line 57 of file TRungeFitter.h.

Referenced by innerwall(), and setBesFromGdml().

int TRungeFitter::_propagation [private]

Definition at line 63 of file TRungeFitter.h.

Referenced by fit(), and propagation().

bool TRungeFitter::_sag [private]

Definition at line 62 of file TRungeFitter.h.

Referenced by fit(), and sag().

bool TRungeFitter::_tof [private]

Definition at line 64 of file TRungeFitter.h.

Referenced by fit(), and tof().

IMagneticFieldSvc* TRungeFitter::m_pmgnIMF [private]

Definition at line 61 of file TRungeFitter.h.

Referenced by fit(), and TRungeFitter().


Generated on Tue Nov 29 23:36:17 2016 for BOSS_7.0.2 by  doxygen 1.4.7