00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifdef TRKRECO_DEBUG_DETAIL
00014 #ifndef TRKRECO_DEBUG
00015 #define TRKRECO_DEBUG
00016 #endif
00017 #endif
00018 #include <float.h>
00019 #include <iostream>
00020 #include "TrkReco/TRungeFitter.h"
00021 #include "TrkReco/TRunge.h"
00022 #include "GaudiKernel/Service.h"
00023 #include "GaudiKernel/ISvcLocator.h"
00024 #include "GaudiKernel/SvcFactory.h"
00025 #define HEP_SHORT_NAMES
00026
00027
00028 #include "MdcTables/MdcTables.h"
00029 #include "CLHEP/Matrix/Vector.h"
00030 #include "CLHEP/Matrix/SymMatrix.h"
00031 #include "CLHEP/Matrix/Matrix.h"
00032 #include "TrkReco/TMLink.h"
00033 #include "TrkReco/TTrackBase.h"
00034 #include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
00035 #include "MdcCalibFunSvc/MdcCalibFunSvc.h"
00036 #include "GaudiKernel/StatusCode.h"
00037 #include "GaudiKernel/IInterface.h"
00038 #include "GaudiKernel/Kernel.h"
00039 #include "GaudiKernel/Service.h"
00040 #include "GaudiKernel/ISvcLocator.h"
00041 #include "GaudiKernel/SvcFactory.h"
00042 #include "GaudiKernel/IDataProviderSvc.h"
00043 #include "GaudiKernel/Bootstrap.h"
00044 #include "GaudiKernel/MsgStream.h"
00045 #include "GaudiKernel/SmartDataPtr.h"
00046 #include "GaudiKernel/IMessageSvc.h"
00047 #include "MdcTables/MdcTables.h"
00048 #include "CLHEP/Matrix/Vector.h"
00049 #include "CLHEP/Matrix/SymMatrix.h"
00050 #include "CLHEP/Matrix/Matrix.h"
00051 #include "TrkReco/TMLink.h"
00052 #include "TrkReco/TTrackBase.h"
00053
00054 #include "GaudiKernel/ISvcLocator.h"
00055 #include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
00056 #include "MdcCalibFunSvc/MdcCalibFunSvc.h"
00057 #include "GaudiKernel/Service.h"
00058 #include "GaudiKernel/ISvcLocator.h"
00059 #include "GaudiKernel/SvcFactory.h"
00060 #include "AIDA/IHistogram1D.h"
00061
00062 #include "AIDA/IHistogram2D.h"
00063 #include "G4Geo/MdcG4Geo.h"
00064 #include "G4Geo/BesG4Geo.h"
00065 #include "G4Material.hh"
00066 #include "G4Tubs.hh"
00067
00068 #include <time.h>
00069
00070 using CLHEP::HepVector;
00071 using CLHEP::HepMatrix;
00072 using CLHEP::HepSymMatrix;
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090 TRungeFitter::TRungeFitter(const std::string& name)
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 }
00099 TRungeFitter::TRungeFitter(const std::string& name,
00100 bool m_sag,int m_prop,bool m_tof)
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 }
00109
00110 TRungeFitter::~TRungeFitter() {
00111 }
00112
00113 void TRungeFitter::sag(bool _in){
00114 _sag = _in;
00115 }
00116 void TRungeFitter::propagation(int _in){
00117 _propagation = _in;
00118 }
00119 void TRungeFitter::tof(bool _in){
00120 _tof = _in;
00121 }
00122 int TRungeFitter::fit(TTrackBase& tb) const{
00123 return fit(tb,0,-1);
00124 }
00125
00126 int TRungeFitter::fit(TTrackBase& tb, float t0Offset, int layer) const{
00127
00128 IMdcCalibFunSvc* l_mdcCalFunSvc;
00129 Gaudi::svcLocator() -> service("MdcCalibFunSvc", l_mdcCalFunSvc);
00130
00131
00132 if(tb.objectType() != Runge) return TFitUnavailable;
00133 TRunge& t = (TRunge&) tb;
00134
00135 if(t.fitted()&&layer==-1) return TFitAlreadyFitted;
00136
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
00152
00153
00154
00155
00156 const HepPoint3D pivot_bak = t.pivot();
00157 t.pivot(ORIGIN);
00158
00159
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
00182 Vector ea(5);
00183
00184 float tof;
00185 HepVector3D p;
00186 unsigned i;
00187
00188 double distance;
00189 double eDistance;
00190
00191
00192 const double ea_init=0.000001;
00193 ea[0]=ea_init;
00194 ea[1]=ea_init;
00195 ea[2]=ea_init;
00196 ea[3]=ea_init;
00197 ea[4]=ea_init;
00198
00199
00200
00201
00202
00203
00204 Vector def_a(t.a());
00205 Vector ta(def_a);
00206
00207
00208 unsigned nTrial = 0;
00209 while(nTrial < 100){
00210
00211
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
00219 for(unsigned j=0;j<5;j++){
00220 ta=def_a;
00221 ta[j]+=ea[j];
00222 t.a(ta);
00223
00224 i=0;
00225
00226 while(TMLink * l = cores[i++]){
00227
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
00233 }
00234
00235
00236
00237 }
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256 t.a(def_a);
00257
00258
00259
00260 i=0;
00261 while(TMLink * l = cores[i++]){
00262 const TMDCWireHit& h = * l->hit();
00263
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
00275 if ((t0Offset == 0.) && (_propagation==0) && (! _tof)) {
00276
00277 distance = l->drift(leftRight);
00278 eDistance = h.dDrift(leftRight);
00279 }else{
00280
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;
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
00332 const double d0 = (onTrack-onWire).mag();
00333
00334
00335 double dDistance = d0 - distance;
00336
00337
00338 for(int j=0;j<5;j++){
00339 dDda[j]=(d[i-1][j]-d0)/ea[j];
00340
00341 }
00342
00343
00344 dchi2da += (dDistance / eDistance2) * dDda;
00345 d2chi2d2a += vT_times_v(dDda) / eDistance2;
00346 double pChi2 = dDistance * dDistance / eDistance2;
00347 chi2 += pChi2;
00348
00349
00350
00351
00352
00353
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 }
00366
00367
00368 double change = chi2Old - chi2;
00369 if (0 <= change && change < 0.01) break;
00370 if (change < 0.) {
00371 factor *= 100;
00372 a += da;
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
00392 return err;
00393 }
00394 alpha2=alpha;
00395 for(unsigned j=0;j<5;j++) alpha2[j][j]*=(1+factor);
00396
00397 da = solve(alpha2,beta);
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407 a -= da;
00408 t.a(a);
00409
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
00418
00419
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
00427 double y_[6]={0,0,0,0,0,0};
00428 t.SetFirst(y_);
00429 int lmass=0;
00430 innerwall(t,lmass,y_);
00431 int nsign=a[2]/abs(a[2]);
00432
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
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
00443 t._fitted = true;
00444 }else{
00445 err = TFitFailed;
00446 }
00447
00448 t._ndf = nCores - dim;
00449 t._chi2 = chi2;
00450
00451
00452 t.pivot(pivot_bak);
00453
00454 delete [] d;
00455
00456
00457 return err;
00458 }
00459 void TRungeFitter::setBesFromGdml(void){
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
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 }
00670
00671 void TRungeFitter::innerwall(TRunge& track, int l_mass,double y[6])const{
00672 for(int i = 0; i < _BesBeamPipeWalls.size(); i++) {
00673 _BesBeamPipeWalls[i].updateTrack(track,y);
00674 }
00675 }