#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) |
const std::string & | name (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 | |
IMagneticFieldSvc * | m_pmgnIMF |
bool | _sag |
int | _propagation |
bool | _tof |
Definition at line 34 of file TRungeFitter.h.
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] |
void TRungeFitter::dump | ( | const std::string & | message = std::string("") , |
|
const std::string & | prefix = std::string("") | |||
) | const |
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 | ) |
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 }
std::vector<RkFitMaterial> TRungeFitter::_BesBeamPipeMaterials |
std::vector<RkFitCylinder> TRungeFitter::_BesBeamPipeWalls |
int TRungeFitter::_propagation [private] |
bool TRungeFitter::_sag [private] |
bool TRungeFitter::_tof [private] |
IMagneticFieldSvc* TRungeFitter::m_pmgnIMF [private] |