#include <MagneticFieldSvc.h>
Inheritance diagram for MagneticFieldSvc:
Public Member Functions | |
| virtual StatusCode | fieldVector (const HepGeom::Point3D< double > &xyz, HepGeom::Vector3D< double > &fvec) const =0 |
| virtual StatusCode | fieldVector (const HepPoint3D &xyz, HepVector3D &fvec) const |
| virtual StatusCode | fieldVector (const HepPoint3D &xyz, HepVector3D &fvec) const |
| virtual StatusCode | finalize () |
| Finalise the service. | |
| virtual StatusCode | finalize () |
| Finalise the service. | |
| virtual double | getReferField () |
| virtual double | getReferField () |
| void | handle (const Incident &) |
| void | handle (const Incident &) |
| virtual bool | ifRealField () const |
| virtual bool | ifRealField () const |
| void | init_params () |
| void | init_params () |
| virtual StatusCode | initialize () |
| Initialise the service (Inherited Service overrides). | |
| virtual StatusCode | initialize () |
| Initialise the service (Inherited Service overrides). | |
| virtual StatusCode | queryInterface (const InterfaceID &riid, void **ppvInterface) |
| virtual StatusCode | queryInterface (const InterfaceID &riid, void **ppvInterface) |
| virtual const InterfaceID & | type () const |
| Service type. | |
| virtual const InterfaceID & | type () const |
| Service type. | |
| virtual StatusCode | uniFieldVector (const HepGeom::Point3D< double > &xyz, HepGeom::Vector3D< double > &fvec) const =0 |
| virtual StatusCode | uniFieldVector (const HepPoint3D &xyz, HepVector3D &fvec) const |
| virtual StatusCode | uniFieldVector (const HepPoint3D &xyz, HepVector3D &fvec) const |
Static Public Member Functions | |
| const InterfaceID & | interfaceID () |
| Retrieve interface ID. | |
| const InterfaceID & | interfaceID () |
| Retrieve interface ID. | |
Protected Member Functions | |
| MagneticFieldSvc (const std::string &name, ISvcLocator *svc) | |
| MagneticFieldSvc (const std::string &name, ISvcLocator *svc) | |
| virtual | ~MagneticFieldSvc () |
| Virtual destructor. | |
| virtual | ~MagneticFieldSvc () |
| Virtual destructor. | |
Private Member Functions | |
| void | fieldGrid (const HepPoint3D &xyz, HepVector3D &fvec) const |
| Fills Q, the field vector. | |
| void | fieldGrid (const HepPoint3D &xyz, HepVector3D &fvec) const |
| Fills Q, the field vector. | |
| void | fieldGrid_TE (const HepPoint3D &xyz, HepVector3D &fvec) const |
| void | fieldGrid_TE (const HepPoint3D &xyz, HepVector3D &fvec) const |
| StatusCode | parseFile () |
| Reads the field map from file. | |
| StatusCode | parseFile () |
| Reads the field map from file. | |
| StatusCode | parseFile_TE () |
| Reads the field map from file. | |
| StatusCode | parseFile_TE () |
| Reads the field map from file. | |
Private Attributes | |
| FieldDBUtil::ConnectionDB * | m_connect_run |
| FieldDBUtil::ConnectionDB * | m_connect_run |
| double | m_Cur_SCQ1_55 |
| double | m_Cur_SCQ1_89 |
| double | m_Cur_SCQ2_10 |
| double | m_Dxyz [3] |
| Steps in x, y and z. | |
| double | m_Dxyz_TE [3] |
| Steps in x, y and z. | |
| IDataProviderSvc * | m_eventSvc |
| IDataProviderSvc * | m_eventSvc |
| std::string | m_filename |
| Magnetic field file name. | |
| std::string | m_filename_TE |
| Magnetic field file name. | |
| int | m_gridDistance |
| grid distance of field map | |
| bool | m_ifRealField |
| double | m_max_FL [3] |
| double | m_max_FL_TE [3] |
| double | m_min_FL [3] |
| double | m_min_FL_TE [3] |
| MucMagneticField * | m_Mucfield |
| MucMagneticField * | m_Mucfield |
| int | m_Nxyz [3] |
| Number of steps in x, y and z. | |
| int | m_Nxyz_TE [3] |
| Number of steps in x, y and z. | |
| int | m_outlevel |
| std::vector< double > | m_P |
| Grid position. | |
| std::vector< double > | m_P |
| Grid position. | |
| std::vector< double > | m_P_1 |
| Grid position. | |
| std::vector< double > | m_P_1 |
| Grid position. | |
| std::vector< double > | m_P_2 |
| Grid position. | |
| std::vector< double > | m_P_2 |
| Grid position. | |
| std::vector< double > | m_P_TE |
| Grid position. | |
| std::vector< double > | m_P_TE |
| Grid position. | |
| std::vector< double > | m_Q |
| Field vector. | |
| std::vector< double > | m_Q |
| Field vector. | |
| std::vector< double > | m_Q_1 |
| Field vector. | |
| std::vector< double > | m_Q_1 |
| Field vector. | |
| std::vector< double > | m_Q_2 |
| Field vector. | |
| std::vector< double > | m_Q_2 |
| Field vector. | |
| std::vector< double > | m_Q_TE |
| Field vector. | |
| std::vector< double > | m_Q_TE |
| Field vector. | |
| int | m_runmode |
| Run mode. | |
| bool | m_useDB |
| double | m_zfield |
| double | m_zOffSet |
| The z offset. | |
| double | m_zOffSet_TE |
| The z offset. | |
Friends | |
| class | SvcFactory<MagneticFieldSvc> |
| Allow SvcFactory to instantiate the service. | |
|
||||||||||||
|
Standard Constructor.
00038 : Service( name, svc )
00039 {
00040 declareProperty( "FieldMapFile", m_filename );
00041 declareProperty( "GridDistance", m_gridDistance = 5);
00042 declareProperty( "RunMode", m_runmode = 2);
00043 declareProperty( "IfRealField", m_ifRealField = true);
00044 declareProperty( "OutLevel", m_outlevel = 1);
00045
00046 declareProperty( "Cur_SCQ1_55", m_Cur_SCQ1_55 = 349.4);
00047 declareProperty( "Cur_SCQ1_89", m_Cur_SCQ1_89 = 426.2);
00048 declareProperty( "Cur_SCQ2_10", m_Cur_SCQ2_10 = 474.2);
00049
00050 declareProperty( "UseDBFlag", m_useDB = true);
00051
00052 m_Mucfield = new MucMagneticField();
00053 if(!m_Mucfield) cout<<"error: can not get MucMagneticField pointer"<<endl;
00054
00055 m_zfield = -1.0*tesla;
00056 }
|
|
|
Virtual destructor.
00061 {
00062 if(m_Mucfield) delete m_Mucfield;
00063 }
|
|
||||||||||||
|
Standard Constructor.
|
|
|
Virtual destructor.
|
|
||||||||||||
|
Fills Q, the field vector.
|
|
||||||||||||
|
Fills Q, the field vector. Linear interpolated field 01694 {
01695
01696 bf[0] = 0.0;
01697 bf[1] = 0.0;
01698 bf[2] = 0.0;
01699
01701 double z = r.z() - m_zOffSet;
01702 if( z < m_min_FL[2] || z > m_max_FL[2] ) return;
01703 double x = r.x();
01704 if( x < m_min_FL[0] || x > m_max_FL[0] ) return;
01705 double y = r.y();
01706 if( y < m_min_FL[1] || y > m_max_FL[1] ) return;
01707 int i = int( (x-m_min_FL[0])/m_Dxyz[0]);
01708 int j = int( (y-m_min_FL[1])/m_Dxyz[1] );
01709 int k = int( (z-m_min_FL[2])/m_Dxyz[2] );
01710
01711 int ijk000 = 3*( m_Nxyz[0]*( m_Nxyz[1]*k + j ) + i );
01712 int ijk001 = 3*( m_Nxyz[0]*( m_Nxyz[1]*(k+1) + j ) + i );
01713 int ijk010 = 3*( m_Nxyz[0]*( m_Nxyz[1]*k + j + 1 ) + i );
01714 int ijk011 = 3*( m_Nxyz[0]*( m_Nxyz[1]*(k+1) + j + 1) + i );
01715 int ijk100 = 3*( m_Nxyz[0]*( m_Nxyz[1]*k + j) + i + 1 );
01716 int ijk101 = 3*( m_Nxyz[0]*( m_Nxyz[1]*(k+1) + j) + i + 1 );
01717 int ijk110 = 3*( m_Nxyz[0]*( m_Nxyz[1]*k + j + 1) + i + 1 );
01718 int ijk111 = 3*( m_Nxyz[0]*( m_Nxyz[1]*(k+1) + j + 1 ) + i + 1 );
01719
01720 // auxiliary variables defined at the vertices of the cube that
01721 // contains the (x, y, z) point where the field is interpolated
01722 /* std::cout<<"x,y,z: "<<x<<","<<y<<","<<z<<std::endl;
01723 std::cout<<"point1(x,y,z): "<<m_P[ijk000]<<","<<m_P[ijk000+1]<<","<<m_P[ijk000+2]<<std::endl;
01724 std::cout<<"point2(x,y,z): "<<m_P[ijk001]<<","<<m_P[ijk001+1]<<","<<m_P[ijk001+2]<<std::endl;
01725 std::cout<<"point3(x,y,z): "<<m_P[ijk010]<<","<<m_P[ijk010+1]<<","<<m_P[ijk010+2]<<std::endl;
01726 std::cout<<"point4(x,y,z): "<<m_P[ijk011]<<","<<m_P[ijk011+1]<<","<<m_P[ijk011+2]<<std::endl;
01727 std::cout<<"point5(x,y,z): "<<m_P[ijk100]<<","<<m_P[ijk100+1]<<","<<m_P[ijk100+2]<<std::endl;
01728 std::cout<<"point6(x,y,z): "<<m_P[ijk101]<<","<<m_P[ijk101+1]<<","<<m_P[ijk101+2]<<std::endl;
01729 std::cout<<"point7(x,y,z): "<<m_P[ijk110]<<","<<m_P[ijk110+1]<<","<<m_P[ijk110+2]<<std::endl;
01730 std::cout<<"point8(x,y,z): "<<m_P[ijk111]<<","<<m_P[ijk111+1]<<","<<m_P[ijk111+2]<<std::endl;
01731
01732 if(std::fabs(m_P[ijk000]-x)>m_Dxyz[0]||std::fabs(m_P[ijk001]-x)>m_Dxyz[0]||std::fabs(m_P[ijk010]-x)>m_Dxyz[0]||std::fabs(m_P[ijk011]-x)>m_Dxyz[0]||std::fabs(m_P[ijk100]-x)>m_Dxyz[0]||std::fabs(m_P[ijk101]-x)>m_Dxyz[0]||std::fabs(m_P[ijk110]-x)>m_Dxyz[0]||std::fabs(m_P[ijk111]-x)>m_Dxyz[0])
01733 std::cout<<"FATALERRORX****************************"<<std::endl;
01734 if(std::fabs(m_P[ijk000+1]-y)>m_Dxyz[1]||std::fabs(m_P[ijk001+1]-y)>m_Dxyz[1]||std::fabs(m_P[ijk010+1]-y)>m_Dxyz[1]||std::fabs(m_P[ijk011+1]-y)>m_Dxyz[1]||std::fabs(m_P[ijk100+1]-y)>m_Dxyz[1]||std::fabs(m_P[ijk101+1]-y)>m_Dxyz[1]||std::fabs(m_P[ijk110+1]-y)>m_Dxyz[1]||std::fabs(m_P[ijk111+1]-y)>m_Dxyz[1])
01735 std::cout<<"FATALERRORY***************************"<<std::endl;
01736 if(std::fabs(m_P[ijk000+2]-z)>m_Dxyz[2]||std::fabs(m_P[ijk001+2]-z)>m_Dxyz[2]||std::fabs(m_P[ijk010+2]-z)>m_Dxyz[2]||std::fabs(m_P[ijk011+2]-z)>m_Dxyz[2]||std::fabs(m_P[ijk100+2]-z)>m_Dxyz[2]||std::fabs(m_P[ijk101+2]-z)>m_Dxyz[2]||std::fabs(m_P[ijk110+2]-z)>m_Dxyz[2]||std::fabs(m_P[ijk111+2]-z)>m_Dxyz[2])
01737 std::cout<<"FATALERRORZ****************************"<<std::endl; */
01738 double cx000 = m_Q[ ijk000 ];
01739 double cx001 = m_Q[ ijk001 ];
01740 double cx010 = m_Q[ ijk010 ];
01741 double cx011 = m_Q[ ijk011 ];
01742 double cx100 = m_Q[ ijk100 ];
01743 double cx101 = m_Q[ ijk101 ];
01744 double cx110 = m_Q[ ijk110 ];
01745 double cx111 = m_Q[ ijk111 ];
01746 double cy000 = m_Q[ ijk000+1 ];
01747 double cy001 = m_Q[ ijk001+1 ];
01748 double cy010 = m_Q[ ijk010+1 ];
01749 double cy011 = m_Q[ ijk011+1 ];
01750 double cy100 = m_Q[ ijk100+1 ];
01751 double cy101 = m_Q[ ijk101+1 ];
01752 double cy110 = m_Q[ ijk110+1 ];
01753 double cy111 = m_Q[ ijk111+1 ];
01754 double cz000 = m_Q[ ijk000+2 ];
01755 double cz001 = m_Q[ ijk001+2 ];
01756 double cz010 = m_Q[ ijk010+2 ];
01757 double cz011 = m_Q[ ijk011+2 ];
01758 double cz100 = m_Q[ ijk100+2 ];
01759 double cz101 = m_Q[ ijk101+2 ];
01760 double cz110 = m_Q[ ijk110+2 ];
01761 double cz111 = m_Q[ ijk111+2 ];
01762 double hx1 = ( x-m_min_FL[0]-i*m_Dxyz[0] )/m_Dxyz[0];
01763 double hy1 = ( y-m_min_FL[1]-j*m_Dxyz[1] )/m_Dxyz[1];
01764 double hz1 = ( z-m_min_FL[2]-k*m_Dxyz[2] )/m_Dxyz[2];
01765 double hx0 = 1.0-hx1;
01766 double hy0 = 1.0-hy1;
01767 double hz0 = 1.0-hz1;
01768 double h000 = hx0*hy0*hz0;
01769 if( fabs(h000) > 0.0 &&
01770 cx000 > 9.0e5 && cy000 > 9.0e5 && cz000 > 9.0e5) return;
01771
01772 double h001 = hx0*hy0*hz1;
01773 if( fabs(h001) > 0.0 &&
01774 cx001 > 9.0e5 && cy001 > 9.0e5 && cz001 > 9.0e5) return;
01775
01776 double h010 = hx0*hy1*hz0;
01777 if( fabs(h010) > 0.0 &&
01778 cx010 > 9.0e5 && cy010 > 9.0e5 && cz010 > 9.0e5) return;
01779
01780 double h011 = hx0*hy1*hz1;
01781 if( fabs(h011) > 0.0 &&
01782 cx011 > 9.0e5 && cy011 > 9.0e5 && cz011 > 9.0e5) return;
01783
01784 double h100 = hx1*hy0*hz0;
01785 if( fabs(h100) > 0.0 &&
01786 cx100 > 9.0e5 && cy100 > 9.0e5 && cz100 > 9.0e5) return;
01787
01788 double h101 = hx1*hy0*hz1;
01789 if( fabs(h101) > 0.0 &&
01790 cx101 > 9.0e5 && cy101 > 9.0e5 && cz101 > 9.0e5) return;
01791
01792 double h110 = hx1*hy1*hz0;
01793 if( fabs(h110) > 0.0 &&
01794 cx110 > 9.0e5 && cy110 > 9.0e5 && cz110 > 9.0e5) return;
01795
01796 double h111 = hx1*hy1*hz1;
01797 if( fabs(h111) > 0.0 &&
01798 cx111 > 9.0e5 && cy111 > 9.0e5 && cz111 > 9.0e5) return;
01799
01800 bf(0) = ( cx000*h000 + cx001*h001 + cx010*h010 + cx011*h011 +
01801 cx100*h100 + cx101*h101 + cx110*h110 + cx111*h111);
01802 bf(1) = ( cy000*h000 + cy001*h001 + cy010*h010 + cy011*h011 +
01803 cy100*h100 + cy101*h101 + cy110*h110 + cy111*h111 );
01804 bf(2) = ( cz000*h000 + cz001*h001 + cz010*h010 + cz011*h011 +
01805 cz100*h100 + cz101*h101 + cz110*h110 + cz111*h111 );
01806 return;
01807 }
|
|
||||||||||||
|
|
|
||||||||||||
|
Linear interpolated field 01813 {
01814
01815 bf[0] = 0.0;
01816 bf[1] = 0.0;
01817 bf[2] = 0.0;
01818
01820 double z = r.z() - m_zOffSet_TE;
01821 if( fabs(z) < m_min_FL_TE[2] || fabs(z) > m_max_FL_TE[2] ) return;
01822 double x = r.x();
01823 if( fabs(x) < m_min_FL_TE[0] || fabs(x) > m_max_FL_TE[0] ) return;
01824 double y = r.y();
01825 if( fabs(y) < m_min_FL_TE[1] || fabs(y) > m_max_FL_TE[1] ) return;
01826 int i = int( (fabs(x)-m_min_FL_TE[0])/m_Dxyz_TE[0]);
01827 int j = int( (fabs(y)-m_min_FL_TE[1])/m_Dxyz_TE[1] );
01828 int k = int( (fabs(z)-m_min_FL_TE[2])/m_Dxyz_TE[2] );
01829
01830 int ijk000 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*k + j ) + i );
01831 int ijk001 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*(k+1) + j ) + i );
01832 int ijk010 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*k + j + 1 ) + i );
01833 int ijk011 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*(k+1) + j + 1) + i );
01834 int ijk100 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*k + j) + i + 1 );
01835 int ijk101 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*(k+1) + j) + i + 1 );
01836 int ijk110 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*k + j + 1) + i + 1 );
01837 int ijk111 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*(k+1) + j + 1 ) + i + 1 );
01838
01839 // auxiliary variables defined at the vertices of the cube that
01840 // contains the (x, y, z) point where the field is interpolated
01841 /* std::cout<<"x,y,z: "<<x<<","<<y<<","<<z<<std::endl;
01842 std::cout<<"point1(x,y,z): "<<m_P[ijk000]<<","<<m_P[ijk000+1]<<","<<m_P[ijk000+2]<<std::endl;
01843 std::cout<<"point2(x,y,z): "<<m_P[ijk001]<<","<<m_P[ijk001+1]<<","<<m_P[ijk001+2]<<std::endl;
01844 std::cout<<"point3(x,y,z): "<<m_P[ijk010]<<","<<m_P[ijk010+1]<<","<<m_P[ijk010+2]<<std::endl;
01845 std::cout<<"point4(x,y,z): "<<m_P[ijk011]<<","<<m_P[ijk011+1]<<","<<m_P[ijk011+2]<<std::endl;
01846 std::cout<<"point5(x,y,z): "<<m_P[ijk100]<<","<<m_P[ijk100+1]<<","<<m_P[ijk100+2]<<std::endl;
01847 std::cout<<"point6(x,y,z): "<<m_P[ijk101]<<","<<m_P[ijk101+1]<<","<<m_P[ijk101+2]<<std::endl;
01848 std::cout<<"point7(x,y,z): "<<m_P[ijk110]<<","<<m_P[ijk110+1]<<","<<m_P[ijk110+2]<<std::endl;
01849 std::cout<<"point8(x,y,z): "<<m_P[ijk111]<<","<<m_P[ijk111+1]<<","<<m_P[ijk111+2]<<std::endl;
01850 */
01851 double cx000 = m_Q_TE[ ijk000 ];
01852 double cx001 = m_Q_TE[ ijk001 ];
01853 double cx010 = m_Q_TE[ ijk010 ];
01854 double cx011 = m_Q_TE[ ijk011 ];
01855 double cx100 = m_Q_TE[ ijk100 ];
01856 double cx101 = m_Q_TE[ ijk101 ];
01857 double cx110 = m_Q_TE[ ijk110 ];
01858 double cx111 = m_Q_TE[ ijk111 ];
01859 double cy000 = m_Q_TE[ ijk000+1 ];
01860 double cy001 = m_Q_TE[ ijk001+1 ];
01861 double cy010 = m_Q_TE[ ijk010+1 ];
01862 double cy011 = m_Q_TE[ ijk011+1 ];
01863 double cy100 = m_Q_TE[ ijk100+1 ];
01864 double cy101 = m_Q_TE[ ijk101+1 ];
01865 double cy110 = m_Q_TE[ ijk110+1 ];
01866 double cy111 = m_Q_TE[ ijk111+1 ];
01867 double cz000 = m_Q_TE[ ijk000+2 ];
01868 double cz001 = m_Q_TE[ ijk001+2 ];
01869 double cz010 = m_Q_TE[ ijk010+2 ];
01870 double cz011 = m_Q_TE[ ijk011+2 ];
01871 double cz100 = m_Q_TE[ ijk100+2 ];
01872 double cz101 = m_Q_TE[ ijk101+2 ];
01873 double cz110 = m_Q_TE[ ijk110+2 ];
01874 double cz111 = m_Q_TE[ ijk111+2 ];
01875 double hx1 = ( fabs(x)-m_min_FL_TE[0]-i*m_Dxyz_TE[0] )/m_Dxyz_TE[0];
01876 double hy1 = ( fabs(y)-m_min_FL_TE[1]-j*m_Dxyz_TE[1] )/m_Dxyz_TE[1];
01877 double hz1 = ( fabs(z)-m_min_FL_TE[2]-k*m_Dxyz_TE[2] )/m_Dxyz_TE[2];
01878 double hx0 = 1.0-hx1;
01879 double hy0 = 1.0-hy1;
01880 double hz0 = 1.0-hz1;
01881 double h000 = hx0*hy0*hz0;
01882 if( fabs(h000) > 0.0 &&
01883 cx000 > 9.0e5 && cy000 > 9.0e5 && cz000 > 9.0e5) return;
01884
01885 double h001 = hx0*hy0*hz1;
01886 if( fabs(h001) > 0.0 &&
01887 cx001 > 9.0e5 && cy001 > 9.0e5 && cz001 > 9.0e5) return;
01888
01889 double h010 = hx0*hy1*hz0;
01890 if( fabs(h010) > 0.0 &&
01891 cx010 > 9.0e5 && cy010 > 9.0e5 && cz010 > 9.0e5) return;
01892
01893 double h011 = hx0*hy1*hz1;
01894 if( fabs(h011) > 0.0 &&
01895 cx011 > 9.0e5 && cy011 > 9.0e5 && cz011 > 9.0e5) return;
01896
01897 double h100 = hx1*hy0*hz0;
01898 if( fabs(h100) > 0.0 &&
01899 cx100 > 9.0e5 && cy100 > 9.0e5 && cz100 > 9.0e5) return;
01900
01901 double h101 = hx1*hy0*hz1;
01902 if( fabs(h101) > 0.0 &&
01903 cx101 > 9.0e5 && cy101 > 9.0e5 && cz101 > 9.0e5) return;
01904
01905 double h110 = hx1*hy1*hz0;
01906 if( fabs(h110) > 0.0 &&
01907 cx110 > 9.0e5 && cy110 > 9.0e5 && cz110 > 9.0e5) return;
01908
01909 double h111 = hx1*hy1*hz1;
01910 if( fabs(h111) > 0.0 &&
01911 cx111 > 9.0e5 && cy111 > 9.0e5 && cz111 > 9.0e5) return;
01912
01913 bf(0) = ( cx000*h000 + cx001*h001 + cx010*h010 + cx011*h011 +
01914 cx100*h100 + cx101*h101 + cx110*h110 + cx111*h111);
01915 bf(1) = ( cy000*h000 + cy001*h001 + cy010*h010 + cy011*h011 +
01916 cy100*h100 + cy101*h101 + cy110*h110 + cy111*h111 );
01917 bf(2) = ( cz000*h000 + cz001*h001 + cz010*h010 + cz011*h011 +
01918 cz100*h100 + cz101*h101 + cz110*h110 + cz111*h111 );
01919
01920
01921 if( r.x() < 0. && r.y() >= 0. && r.z() > 0. ){
01922 bf(0) = -bf(0);
01923 }
01924 else if( r.x() <= 0. && r.y() < 0. && r.z() > 0. ){
01925 bf(0) = -bf(0);
01926 bf(1) = -bf(1);
01927 }
01928 else if( r.x() > 0. && r.y() < 0. && r.z() > 0. ){
01929 bf(1) = -bf(1);
01930 }
01931
01932 else if( r.x() >= 0. && r.y() > 0. && r.z() < 0. ){
01933 bf(0) = -bf(0);
01934 bf(1) = -bf(1);
01935 }
01936 else if( r.x() < 0. && r.y() >= 0. && r.z() < 0. ){
01937 bf(1) = -bf(1);
01938 }
01939 else if( r.x() <= 0. && r.y() < 0. && r.z() < 0. ){
01940 bf(0) = bf(0);
01941 bf(1) = bf(1);
01942 }
01943 else if( r.x() > 0. && r.y() <= 0. && r.z() < 0. ){
01944 bf(0) = -bf(0);
01945 }
01946
01947 return;
01948 }
|
|
||||||||||||
|
|
|
||||||||||||
|
IMagneticFieldSvc interface.
|
|
||||||||||||
|
IMagneticFieldSvc interface.
00737 {
00738 //reference frames defined by simulation and SCM are different. In simulation: x--south to north, y--down to up, but in SCM: x--north to south, y--down to up
00739 double new_x = -newr.x();
00740 double new_y = newr.y();
00741 double new_z = -newr.z();
00742 HepPoint3D r(new_x,new_y,new_z);
00743 HepVector3D b;
00744 b[0]=0.0*tesla;
00745 b[1]=0.0*tesla;
00746 b[2]=0.0*tesla;
00747 // This routine is now dummy. Was previously converting to/from CLHEP units
00748 if(-2.1*m<r.z() && r.z()<2.1*m && -1.8*m<r.x() && r.x()<1.8*m && -1.8*m<r.y() && r.y()<1.8*m)
00749 {
00750 if(-1.2*m<r.z()&&r.z()<1.2*m&&0*m<=std::sqrt(r.x()*r.x()+r.y()*r.y())&&std::sqrt(r.x()*r.x()+r.y()*r.y())<0.9*m)
00751 //if(std::abs(r.z())<1.2*m&&std::abs(r.x())<=0.9*m&&std::abs(r.y())<=0.9*m)
00752 {
00753 this->fieldGrid( r, b );
00754 }
00755 else
00756 {
00757 this->fieldGrid_TE( r, b );
00758 }
00759 }
00760 if((fabs(r.z())<=1970*mm && sqrt(r.x()*r.x()+r.y()*r.y()) >= 1740*mm) || (fabs(r.z())>=2050*mm))
00761 {
00762 HepPoint3D mr;
00763 HepVector3D tb;
00764 int part,layer,mat;
00765 double theta;
00766 bool ifbar = false;
00767 bool ifend = false;
00768 if(r.x()!=0.){
00769 theta = atan2(fabs(r.y()),fabs(r.x()));
00770 if(0<=theta&&theta<pi/8) {
00771 mr[0] = fabs(r.x()); mr[1] = -fabs(r.y()); mr[2] = fabs(r.z());
00772 if(mr[2]<=1970*mm&&1740*mm<=mr[0]&&mr[0]<=2620*mm){
00773 part = 1;
00774 if(1740*mm<=mr[0]&&mr[0]<1770*mm) { layer = 0; mat = 0; }
00775 if(1770*mm<=mr[0]&&mr[0]<1810*mm) { layer = 0; mat = 1; }
00776 if(1810*mm<=mr[0]&&mr[0]<1840*mm) { layer = 1; mat = 0; }
00777 if(1840*mm<=mr[0]&&mr[0]<1880*mm) { layer = 1; mat = 1; }
00778 if(1880*mm<=mr[0]&&mr[0]<1910*mm) { layer = 2; mat = 0; }
00779 if(1910*mm<=mr[0]&&mr[0]<1950*mm) { layer = 2; mat = 1; }
00780 if(1950*mm<=mr[0]&&mr[0]<1990*mm) { layer = 3; mat = 0; }
00781 if(1990*mm<=mr[0]&&mr[0]<2030*mm) { layer = 3; mat = 1; }
00782 if(2030*mm<=mr[0]&&mr[0]<2070*mm) { layer = 4; mat = 0; }
00783 if(2070*mm<=mr[0]&&mr[0]<2110*mm) { layer = 4; mat = 1; }
00784 if(2110*mm<=mr[0]&&mr[0]<2190*mm) { layer = 5; mat = 0; }
00785 if(2190*mm<=mr[0]&&mr[0]<2230*mm) { layer = 5; mat = 1; }
00786 if(2230*mm<=mr[0]&&mr[0]<2310*mm) { layer = 6; mat = 0; }
00787 if(2310*mm<=mr[0]&&mr[0]<2350*mm) { layer = 6; mat = 1; }
00788 if(2350*mm<=mr[0]&&mr[0]<2430*mm) { layer = 7; mat = 0; }
00789 if(2430*mm<=mr[0]&&mr[0]<2470*mm) { layer = 7; mat = 1; }
00790 if(2470*mm<=mr[0]&&mr[0]<=2620*mm) { layer = 8; mat = 0; }
00791 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00792 b[0] = tb[0];
00793 b[1] = -tb[1];
00794 b[2] = tb[2];
00795 ifbar = true;
00796 }
00797 if(2050*mm<=mr[2]&&mr[2]<2090*mm&&1034*mm<=mr[0]&&mr[0]<=2500*mm){
00798 part = 0; layer = 0; mat = 0;
00799 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00800 b[0] = tb[0];
00801 b[1] = -tb[1];
00802 b[2] = tb[2];
00803 ifend = true;
00804 }
00805 if(2090*mm<=mr[2]&&mr[2]<2130*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
00806 part = 0; layer = 0; mat = 1;
00807 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00808 b[0] = tb[0];
00809 b[1] = -tb[1];
00810 b[2] = tb[2];
00811 ifend = true;
00812 }
00813 if(2130*mm<=mr[2]&&mr[2]<2170*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
00814 part = 0; layer = 1; mat = 0;
00815 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00816 b[0] = tb[0];
00817 b[1] = -tb[1];
00818 b[2] = tb[2];
00819 ifend = true;
00820 }
00821 if(2170*mm<=mr[2]&&mr[2]<2210*mm&&1100*mm<=mr[0]&&mr[0]<=2500*mm) {
00822 part = 0; layer = 1; mat = 1;
00823 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00824 b[0] = tb[0];
00825 b[1] = -tb[1];
00826 b[2] = tb[2];
00827 ifend = true;
00828 }
00829 if(2210*mm<=mr[2]&&mr[2]<2240*mm&&1100*mm<mr[0]&&mr[0]<=2500*mm) {
00830 part = 0; layer = 2; mat = 0;
00831 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00832 b[0] = tb[0];
00833 b[1] = -tb[1];
00834 b[2] = tb[2];
00835 ifend = true;
00836 }
00837 if(2240*mm<=mr[2]&&mr[2]<2280*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
00838 part = 0; layer = 2; mat = 1;
00839 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00840 b[0] = tb[0];
00841 b[1] = -tb[1];
00842 b[2] = tb[2];
00843 ifend = true;
00844 }
00845 if(2280*mm<=mr[2]&&mr[2]<2310*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
00846 part = 0; layer = 3; mat = 0;
00847 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00848 b[0] = tb[0];
00849 b[1] = -tb[1];
00850 b[2] = tb[2];
00851 ifend = true;
00852 }
00853 if(2310*mm<=mr[2]&&mr[2]<2350*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
00854 part = 0; layer = 3; mat = 1;
00855 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00856 b[0] = tb[0];
00857 b[1] = -tb[1];
00858 b[2] = tb[2];
00859 ifend = true;
00860 }
00861 if(2350*mm<=mr[2]&&mr[2]<2380*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
00862 part = 0; layer = 4; mat = 0;
00863 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00864 b[0] = tb[0];
00865 b[1] = -tb[1];
00866 b[2] = tb[2];
00867 ifend = true;
00868 }
00869 if(2380*mm<=mr[2]&&mr[2]<2420*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
00870 part = 0; layer = 4; mat = 1;
00871 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00872 b[0] = tb[0];
00873 b[1] = -tb[1];
00874 b[2] = tb[2];
00875 ifend = true;
00876 }
00877 if(2420*mm<=mr[2]&&mr[2]<2470*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
00878 part = 0; layer = 5; mat = 0;
00879 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00880 b[0] = tb[0];
00881 b[1] = -tb[1];
00882 b[2] = tb[2];
00883 ifend = true;
00884 }
00885 if(2470*mm<=mr[2]&&mr[2]<2510*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
00886 part = 0; layer = 5; mat =1;
00887 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00888 b[0] = tb[0];
00889 b[1] = -tb[1];
00890 b[2] = tb[2];
00891 ifend = true;
00892 }
00893 if(2510*mm<=mr[2]&&mr[2]<2590*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
00894 part = 0; layer = 6; mat = 0;
00895 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00896 b[0] = tb[0];
00897 b[1] = -tb[1];
00898 b[2] = tb[2];
00899 ifend = true;
00900 }
00901 if(2590*mm<=mr[2]&&mr[2]<2630*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
00902 part = 0; layer = 6; mat = 1;
00903 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00904 b[0] = tb[0];
00905 b[1] = -tb[1];
00906 b[2] = tb[2];
00907 ifend = true;
00908 }
00909 if(2630*mm<=mr[2]&&mr[2]<2710*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
00910 part = 0; layer = 7; mat = 0;
00911 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00912 b[0] = tb[0];
00913 b[1] = -tb[1];
00914 b[2] = tb[2];
00915 ifend = true;
00916 }
00917 if(2710*mm<=mr[2]&&mr[2]<2750*mm&&1362*mm<=mr[0]&&mr[0]<=2500*mm) {
00918 part = 0; layer = 7; mat = 1;
00919 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00920 b[0] = tb[0];
00921 b[1] = -tb[1];
00922 b[2] = tb[2];
00923 ifend = true;
00924 }
00925 if(2750*mm<=mr[2]&&mr[2]<=2800*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
00926 part = 0; layer = 8; mat = 0;
00927 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00928 b[0] = tb[0];
00929 b[1] = -tb[1];
00930 b[2] = tb[2];
00931 ifend = true;
00932 }
00933 }
00934 if(pi/8<=theta&&theta<pi/4) {
00935 mr[0] = fabs(r.x())*cos(pi/4)+fabs(r.y())*sin(pi/4);
00936 mr[1] = -fabs(r.x())*sin(pi/4)+fabs(r.y())*cos(pi/4);
00937 mr[2] = fabs(r.z());
00938 if(mr[2]<=1970*mm&&1740*mm<=mr[0]&&mr[0]<=2620*mm) {
00939 part = 1;
00940 if(1740*mm<=mr[0]&&mr[0]<1770*mm) { layer = 0; mat = 0; }
00941 if(1770*mm<=mr[0]&&mr[0]<1810*mm) { layer = 0; mat = 1; }
00942 if(1810*mm<=mr[0]&&mr[0]<1840*mm) { layer = 1; mat = 0; }
00943 if(1840*mm<=mr[0]&&mr[0]<1880*mm) { layer = 1; mat = 1; }
00944 if(1880*mm<=mr[0]&&mr[0]<1910*mm) { layer = 2; mat = 0; }
00945 if(1910*mm<=mr[0]&&mr[0]<1950*mm) { layer = 2; mat = 1; }
00946 if(1950*mm<=mr[0]&&mr[0]<1990*mm) { layer = 3; mat = 0; }
00947 if(1990*mm<=mr[0]&&mr[0]<2030*mm) { layer = 3; mat = 1; }
00948 if(2030*mm<=mr[0]&&mr[0]<2070*mm) { layer = 4; mat = 0; }
00949 if(2070*mm<=mr[0]&&mr[0]<2110*mm) { layer = 4; mat = 1; }
00950 if(2110*mm<=mr[0]&&mr[0]<2190*mm) { layer = 5; mat = 0; }
00951 if(2190*mm<=mr[0]&&mr[0]<2230*mm) { layer = 5; mat = 1; }
00952 if(2230*mm<=mr[0]&&mr[0]<2310*mm) { layer = 6; mat = 0; }
00953 if(2310*mm<=mr[0]&&mr[0]<2350*mm) { layer = 6; mat = 1; }
00954 if(2350*mm<=mr[0]&&mr[0]<2430*mm) { layer = 7; mat = 0; }
00955 if(2430*mm<=mr[0]&&mr[0]<2470*mm) { layer = 7; mat = 1; }
00956 if(2470*mm<=mr[0]&&mr[0]<=2620*mm) { layer = 8; mat = 0; }
00957 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00958 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
00959 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
00960 b[2] = tb[2];
00961 ifbar = true;
00962 }
00963 if(2050*mm<=mr[2]&&mr[2]<2090*mm&&1034*mm<=mr[0]&&mr[0]<=2500*mm){
00964 part = 0; layer = 0; mat = 0;
00965 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00966 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
00967 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
00968 b[2] = tb[2];
00969 ifend = true;
00970 }
00971 if(2090*mm<=mr[2]&&mr[2]<2130*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
00972 part = 0; layer = 0; mat = 1;
00973 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00974 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
00975 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
00976 b[2] = tb[2];
00977 ifend = true;
00978 }
00979 if(2130*mm<=mr[2]&&mr[2]<2170*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
00980 part = 0; layer = 1; mat = 0;
00981 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00982 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
00983 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
00984 b[2] = tb[2];
00985 ifend = true;
00986 }
00987 if(2170*mm<=mr[2]&&mr[2]<2210*mm&&1100*mm<=mr[0]&&mr[0]<=2500*mm) {
00988 part = 0; layer = 1; mat = 1;
00989 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00990 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
00991 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
00992 b[2] = tb[2];
00993 ifend = true;
00994 }
00995 if(2210*mm<=mr[2]&&mr[2]<2240*mm&&1100*mm<mr[0]&&mr[0]<=2500*mm) {
00996 part = 0; layer = 2; mat = 0;
00997 m_Mucfield->getMucField(part,layer,mat,mr,tb);
00998 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
00999 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01000 b[2] = tb[2];
01001 ifend = true;
01002 }
01003 if(2240*mm<=mr[2]&&mr[2]<2280*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
01004 part = 0; layer = 2; mat = 1;
01005 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01006 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01007 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01008 b[2] = tb[2];
01009 ifend = true;
01010 }
01011 if(2280*mm<=mr[2]&&mr[2]<2310*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
01012 part = 0; layer = 3; mat = 0;
01013 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01014 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01015 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01016 b[2] = tb[2];
01017 ifend = true;
01018 }
01019 if(2310*mm<=mr[2]&&mr[2]<2350*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
01020 part = 0; layer = 3; mat = 1;
01021 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01022 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01023 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01024 b[2] = tb[2];
01025 ifend = true;
01026 }
01027 if(2350*mm<=mr[2]&&mr[2]<2380*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
01028 part = 0; layer = 4; mat = 0;
01029 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01030 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01031 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01032 b[2] = tb[2];
01033 ifend = true;
01034 }
01035 if(2380*mm<=mr[2]&&mr[2]<2420*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
01036 part = 0; layer = 4; mat = 1;
01037 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01038 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01039 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01040 b[2] = tb[2];
01041 ifend = true;
01042 }
01043 if(2420*mm<=mr[2]&&mr[2]<2470*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
01044 part = 0; layer = 5; mat = 0;
01045 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01046 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01047 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01048 b[2] = tb[2];
01049 ifend = true;
01050 }
01051 if(2470*mm<=mr[2]&&mr[2]<2510*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
01052 part = 0; layer = 5; mat =1;
01053 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01054 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01055 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01056 b[2] = tb[2];
01057 ifend = true;
01058 }
01059 if(2510*mm<=mr[2]&&mr[2]<2590*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
01060 part = 0; layer = 6; mat = 0;
01061 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01062 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01063 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01064 b[2] = tb[2];
01065 ifend = true;
01066 }
01067 if(2590*mm<=mr[2]&&mr[2]<2630*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01068 part = 0; layer = 6; mat = 1;
01069 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01070 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01071 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01072 b[2] = tb[2];
01073 ifend = true;
01074 }
01075 if(2630*mm<=mr[2]&&mr[2]<2710*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01076 part = 0; layer = 7; mat = 0;
01077 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01078 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01079 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01080 b[2] = tb[2];
01081 ifend = true;
01082 }
01083 if(2710*mm<=mr[2]&&mr[2]<2750*mm&&1362*mm<=mr[0]&&mr[0]<=2500*mm) {
01084 part = 0; layer = 7; mat = 1;
01085 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01086 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01087 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01088 b[2] = tb[2];
01089 ifend = true;
01090 }
01091 if(2750*mm<=mr[2]&&mr[2]<=2800*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01092 part = 0; layer = 8; mat = 0;
01093 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01094 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01095 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01096 b[2] = tb[2];
01097 ifend = true;
01098 }
01099 }
01100 if(pi/4<=theta&&theta<3*pi/8) {
01101 mr[0] = fabs(r.x())*cos(pi/4)+fabs(r.y())*sin(pi/4);
01102 mr[1] = fabs(r.x())*sin(pi/4)-fabs(r.y())*cos(pi/4);
01103 mr[2] = fabs(r.z());
01104 if(mr[2]<=1970*mm&&1740*mm<=mr[0]&&mr[0]<=2620*mm) {
01105 part = 1;
01106 if(1740*mm<=mr[0]&&mr[0]<1770*mm) { layer = 0; mat = 0; }
01107 if(1770*mm<=mr[0]&&mr[0]<1810*mm) { layer = 0; mat = 1; }
01108 if(1810*mm<=mr[0]&&mr[0]<1840*mm) { layer = 1; mat = 0; }
01109 if(1840*mm<=mr[0]&&mr[0]<1880*mm) { layer = 1; mat = 1; }
01110 if(1880*mm<=mr[0]&&mr[0]<1910*mm) { layer = 2; mat = 0; }
01111 if(1910*mm<=mr[0]&&mr[0]<1950*mm) { layer = 2; mat = 1; }
01112 if(1950*mm<=mr[0]&&mr[0]<1990*mm) { layer = 3; mat = 0; }
01113 if(1990*mm<=mr[0]&&mr[0]<2030*mm) { layer = 3; mat = 1; }
01114 if(2030*mm<=mr[0]&&mr[0]<2070*mm) { layer = 4; mat = 0; }
01115 if(2070*mm<=mr[0]&&mr[0]<2110*mm) { layer = 4; mat = 1; }
01116 if(2110*mm<=mr[0]&&mr[0]<2190*mm) { layer = 5; mat = 0; }
01117 if(2190*mm<=mr[0]&&mr[0]<2230*mm) { layer = 5; mat = 1; }
01118 if(2230*mm<=mr[0]&&mr[0]<2310*mm) { layer = 6; mat = 0; }
01119 if(2310*mm<=mr[0]&&mr[0]<2350*mm) { layer = 6; mat = 1; }
01120 if(2350*mm<=mr[0]&&mr[0]<2430*mm) { layer = 7; mat = 0; }
01121 if(2430*mm<=mr[0]&&mr[0]<2470*mm) { layer = 7; mat = 1; }
01122 if(2470*mm<=mr[0]&&mr[0]<=2620*mm) { layer = 8; mat = 0; }
01123 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01124 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01125 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01126 b[2] = tb[2];
01127 ifbar = true;
01128 }
01129 if(2050*mm<=mr[2]&&mr[2]<2090*mm&&1034*mm<=mr[0]&&mr[0]<=2500*mm){
01130 part = 0; layer = 0; mat = 0;
01131 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01132 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01133 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01134 b[2] = tb[2];
01135 ifend = true;
01136 }
01137 if(2090*mm<=mr[2]&&mr[2]<2130*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
01138 part = 0; layer = 0; mat = 1;
01139 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01140 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01141 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01142 b[2] = tb[2];
01143 ifend = true;
01144 }
01145 if(2130*mm<=mr[2]&&mr[2]<2170*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
01146 part = 0; layer = 1; mat = 0;
01147 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01148 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01149 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01150 b[2] = tb[2];
01151 ifend = true;
01152 }
01153 if(2170*mm<=mr[2]&&mr[2]<2210*mm&&1100*mm<=mr[0]&&mr[0]<=2500*mm) {
01154 part = 0; layer = 1; mat = 1;
01155 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01156 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01157 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01158 b[2] = tb[2];
01159 ifend = true;
01160 }
01161 if(2210*mm<=mr[2]&&mr[2]<2240*mm&&1100*mm<mr[0]&&mr[0]<=2500*mm) {
01162 part = 0; layer = 2; mat = 0;
01163 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01164 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01165 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01166 b[2] = tb[2];
01167 ifend = true;
01168 }
01169 if(2240*mm<=mr[2]&&mr[2]<2280*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
01170 part = 0; layer = 2; mat = 1;
01171 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01172 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01173 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01174 b[2] = tb[2];
01175 ifend = true;
01176 }
01177 if(2280*mm<=mr[2]&&mr[2]<2310*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
01178 part = 0; layer = 3; mat = 0;
01179 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01180 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01181 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01182 b[2] = tb[2];
01183 ifend = true;
01184 }
01185 if(2310*mm<=mr[2]&&mr[2]<2350*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
01186 part = 0; layer = 3; mat = 1;
01187 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01188 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01189 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01190 b[2] = tb[2];
01191 ifend = true;
01192 }
01193 if(2350*mm<=mr[2]&&mr[2]<2380*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
01194 part = 0; layer = 4; mat = 0;
01195 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01196 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01197 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01198 b[2] = tb[2];
01199 ifend = true;
01200 }
01201 if(2380*mm<=mr[2]&&mr[2]<2420*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
01202 part = 0; layer = 4; mat = 1;
01203 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01204 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01205 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01206 b[2] = tb[2];
01207 ifend = true;
01208 }
01209 if(2420*mm<=mr[2]&&mr[2]<2470*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
01210 part = 0; layer = 5; mat = 0;
01211 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01212 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01213 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01214 b[2] = tb[2];
01215 ifend = true;
01216 }
01217 if(2470*mm<=mr[2]&&mr[2]<2510*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
01218 part = 0; layer = 5; mat =1;
01219 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01220 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01221 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01222 b[2] = tb[2];
01223 ifend = true;
01224 }
01225 if(2510*mm<=mr[2]&&mr[2]<2590*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
01226 part = 0; layer = 6; mat = 0;
01227 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01228 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01229 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01230 b[2] = tb[2];
01231 ifend = true;
01232 }
01233 if(2590*mm<=mr[2]&&mr[2]<2630*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01234 part = 0; layer = 6; mat = 1;
01235 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01236 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01237 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01238 b[2] = tb[2];
01239 ifend = true;
01240 }
01241 if(2630*mm<=mr[2]&&mr[2]<2710*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01242 part = 0; layer = 7; mat = 0;
01243 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01244 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01245 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01246 b[2] = tb[2];
01247 ifend = true;
01248 }
01249 if(2710*mm<=mr[2]&&mr[2]<2750*mm&&1362*mm<=mr[0]&&mr[0]<=2500*mm) {
01250 part = 0; layer = 7; mat = 1;
01251 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01252 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01253 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01254 b[2] = tb[2];
01255 ifend = true;
01256 }
01257 if(2750*mm<=mr[2]&&mr[2]<=2800*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01258 part = 0; layer = 8; mat = 0;
01259 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01260 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01261 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01262 b[2] = tb[2];
01263 ifend = true;
01264 }
01265 }
01266 if(3*pi/8<=theta&&theta<pi/2) {
01267 mr[0] = fabs(r.y()); mr[1] = -fabs(r.x()); mr[2] = fabs(r.z());
01268 if(mr[2]<=1970*mm&&1740*mm<=mr[0]&&mr[0]<=2620*mm) {
01269 part = 1;
01270 if(1740*mm<=mr[0]&&mr[0]<1770*mm) { layer = 0; mat = 0; }
01271 if(1770*mm<=mr[0]&&mr[0]<1810*mm) { layer = 0; mat = 1; }
01272 if(1810*mm<=mr[0]&&mr[0]<1840*mm) { layer = 1; mat = 0; }
01273 if(1840*mm<=mr[0]&&mr[0]<1880*mm) { layer = 1; mat = 1; }
01274 if(1880*mm<=mr[0]&&mr[0]<1910*mm) { layer = 2; mat = 0; }
01275 if(1910*mm<=mr[0]&&mr[0]<1950*mm) { layer = 2; mat = 1; }
01276 if(1950*mm<=mr[0]&&mr[0]<1990*mm) { layer = 3; mat = 0; }
01277 if(1990*mm<=mr[0]&&mr[0]<2030*mm) { layer = 3; mat = 1; }
01278 if(2030*mm<=mr[0]&&mr[0]<2070*mm) { layer = 4; mat = 0; }
01279 if(2070*mm<=mr[0]&&mr[0]<2110*mm) { layer = 4; mat = 1; }
01280 if(2110*mm<=mr[0]&&mr[0]<2190*mm) { layer = 5; mat = 0; }
01281 if(2190*mm<=mr[0]&&mr[0]<2230*mm) { layer = 5; mat = 1; }
01282 if(2230*mm<=mr[0]&&mr[0]<2310*mm) { layer = 6; mat = 0; }
01283 if(2310*mm<=mr[0]&&mr[0]<2350*mm) { layer = 6; mat = 1; }
01284 if(2350*mm<=mr[0]&&mr[0]<2430*mm) { layer = 7; mat = 0; }
01285 if(2430*mm<=mr[0]&&mr[0]<2470*mm) { layer = 7; mat = 1; }
01286 if(2470*mm<=mr[0]&&mr[0]<=2620*mm) { layer = 8; mat = 0; }
01287 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01288 b[0] = -tb[1];
01289 b[1] = tb[0];
01290 b[2] = tb[2];
01291 ifbar = true;
01292 }
01293 if(2050*mm<=mr[2]&&mr[2]<2090*mm&&1034*mm<=mr[0]&&mr[0]<=2500*mm){
01294 part = 0; layer = 0; mat = 0;
01295 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01296 b[0] = -tb[1];
01297 b[1] = tb[0];
01298 b[2] = tb[2];
01299 ifend = true;
01300 }
01301 if(2090*mm<=mr[2]&&mr[2]<2130*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
01302 part = 0; layer = 0; mat = 1;
01303 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01304 b[0] = -tb[1];
01305 b[1] = tb[0];
01306 b[2] = tb[2];
01307 ifend = true;
01308 }
01309 if(2130*mm<=mr[2]&&mr[2]<2170*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
01310 part = 0; layer = 1; mat = 0;
01311 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01312 b[0] = -tb[1];
01313 b[1] = tb[0];
01314 b[2] = tb[2];
01315 ifend = true;
01316 }
01317 if(2170*mm<=mr[2]&&mr[2]<2210*mm&&1100*mm<=mr[0]&&mr[0]<=2500*mm) {
01318 part = 0; layer = 1; mat = 1;
01319 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01320 b[0] = -tb[1];
01321 b[1] = tb[0];
01322 b[2] = tb[2];
01323 ifend = true;
01324 }
01325 if(2210*mm<=mr[2]&&mr[2]<2240*mm&&1100*mm<mr[0]&&mr[0]<=2500*mm) {
01326 part = 0; layer = 2; mat = 0;
01327 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01328 b[0] = -tb[1];
01329 b[1] = tb[0];
01330 b[2] = tb[2];
01331 ifend = true;
01332 }
01333 if(2240*mm<=mr[2]&&mr[2]<2280*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
01334 part = 0; layer = 2; mat = 1;
01335 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01336 b[0] = -tb[1];
01337 b[1] = tb[0];
01338 b[2] = tb[2];
01339 ifend = true;
01340 }
01341 if(2280*mm<=mr[2]&&mr[2]<2310*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
01342 part = 0; layer = 3; mat = 0;
01343 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01344 b[0] = -tb[1];
01345 b[1] = tb[0];
01346 b[2] = tb[2];
01347 ifend = true;
01348 }
01349 if(2310*mm<=mr[2]&&mr[2]<2350*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
01350 part = 0; layer = 3; mat = 1;
01351 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01352 b[0] = -tb[1];
01353 b[1] = tb[0];
01354 b[2] = tb[2];
01355 ifend = true;
01356 }
01357 if(2350*mm<=mr[2]&&mr[2]<2380*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
01358 part = 0; layer = 4; mat = 0;
01359 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01360 b[0] = -tb[1];
01361 b[1] = tb[0];
01362 b[2] = tb[2];
01363 ifend = true;
01364 }
01365 if(2380*mm<=mr[2]&&mr[2]<2420*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
01366 part = 0; layer = 4; mat = 1;
01367 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01368 b[0] = -tb[1];
01369 b[1] = tb[0];
01370 b[2] = tb[2];
01371 ifend = true;
01372 }
01373 if(2420*mm<=mr[2]&&mr[2]<2470*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
01374 part = 0; layer = 5; mat = 0;
01375 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01376 b[0] = -tb[1];
01377 b[1] = tb[0];
01378 b[2] = tb[2];
01379 ifend = true;
01380 }
01381 if(2470*mm<=mr[2]&&mr[2]<2510*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
01382 part = 0; layer = 5; mat =1;
01383 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01384 b[0] = -tb[1];
01385 b[1] = tb[0];
01386 b[2] = tb[2];
01387 ifend = true;
01388 }
01389 if(2510*mm<=mr[2]&&mr[2]<2590*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
01390 part = 0; layer = 6; mat = 0;
01391 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01392 b[0] = -tb[1];
01393 b[1] = tb[0];
01394 b[2] = tb[2];
01395 ifend = true;
01396 }
01397 if(2590*mm<=mr[2]&&mr[2]<2630*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01398 part = 0; layer = 6; mat = 1;
01399 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01400 b[0] = -tb[1];
01401 b[1] = tb[0];
01402 b[2] = tb[2];
01403 ifend = true;
01404 }
01405 if(2630*mm<=mr[2]&&mr[2]<2710*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01406 part = 0; layer = 7; mat = 0;
01407 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01408 b[0] = -tb[1];
01409 b[1] = tb[0];
01410 b[2] = tb[2];
01411 ifend = true;
01412 }
01413 if(2710*mm<=mr[2]&&mr[2]<2750*mm&&1362*mm<=mr[0]&&mr[0]<=2500*mm) {
01414 part = 0; layer = 7; mat = 1;
01415 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01416 b[0] = -tb[1];
01417 b[1] = tb[0];
01418 b[2] = tb[2];
01419 ifend = true;
01420 }
01421 if(2750*mm<=mr[2]&&mr[2]<=2800*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01422 part = 0; layer = 8; mat = 0;
01423 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01424 b[0] = -tb[1];
01425 b[1] = tb[0];
01426 b[2] = tb[2];
01427 ifend = true;
01428 }
01429 }
01430 }
01431 if(r.x()==0.) {
01432 mr[0] = fabs(r.y()); mr[1] = -fabs(r.x()); mr[2] = fabs(r.z());
01433 if(mr[2]<=1970*mm&&1740*mm<=mr[0]&&mr[0]<=2620*mm) {
01434 part = 1;
01435 if(1740*mm<=mr[0]&&mr[0]<1770*mm) { layer = 0; mat = 0; }
01436 if(1770*mm<=mr[0]&&mr[0]<1810*mm) { layer = 0; mat = 1; }
01437 if(1810*mm<=mr[0]&&mr[0]<1840*mm) { layer = 1; mat = 0; }
01438 if(1840*mm<=mr[0]&&mr[0]<1880*mm) { layer = 1; mat = 1; }
01439 if(1880*mm<=mr[0]&&mr[0]<1910*mm) { layer = 2; mat = 0; }
01440 if(1910*mm<=mr[0]&&mr[0]<1950*mm) { layer = 2; mat = 1; }
01441 if(1950*mm<=mr[0]&&mr[0]<1990*mm) { layer = 3; mat = 0; }
01442 if(1990*mm<=mr[0]&&mr[0]<2030*mm) { layer = 3; mat = 1; }
01443 if(2030*mm<=mr[0]&&mr[0]<2070*mm) { layer = 4; mat = 0; }
01444 if(2070*mm<=mr[0]&&mr[0]<2110*mm) { layer = 4; mat = 1; }
01445 if(2110*mm<=mr[0]&&mr[0]<2190*mm) { layer = 5; mat = 0; }
01446 if(2190*mm<=mr[0]&&mr[0]<2230*mm) { layer = 5; mat = 1; }
01447 if(2230*mm<=mr[0]&&mr[0]<2310*mm) { layer = 6; mat = 0; }
01448 if(2310*mm<=mr[0]&&mr[0]<2350*mm) { layer = 6; mat = 1; }
01449 if(2350*mm<=mr[0]&&mr[0]<2430*mm) { layer = 7; mat = 0; }
01450 if(2430*mm<=mr[0]&&mr[0]<2470*mm) { layer = 7; mat = 1; }
01451 if(2470*mm<=mr[0]&&mr[0]<=2620*mm) { layer = 8; mat = 0; }
01452 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01453 b[0] = -tb[1];
01454 b[1] = tb[0];
01455 b[2] = tb[2];
01456 ifbar = true;
01457 }
01458 if(2050*mm<=mr[2]&&mr[2]<2090*mm&&1034*mm<=mr[0]&&mr[0]<=2500*mm){
01459 part = 0; layer = 0; mat = 0;
01460 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01461 b[0] = -tb[1];
01462 b[1] = tb[0];
01463 b[2] = tb[2];
01464 ifend = true;
01465 }
01466 if(2090*mm<=mr[2]&&mr[2]<2130*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
01467 part = 0; layer = 0; mat = 1;
01468 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01469 b[0] = -tb[1];
01470 b[1] = tb[0];
01471 b[2] = tb[2];
01472 ifend = true;
01473 }
01474 if(2130*mm<=mr[2]&&mr[2]<2170*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
01475 part = 0; layer = 1; mat = 0;
01476 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01477 b[0] = -tb[1];
01478 b[1] = tb[0];
01479 b[2] = tb[2];
01480 ifend = true;
01481 }
01482 if(2170*mm<=mr[2]&&mr[2]<2210*mm&&1100*mm<=mr[0]&&mr[0]<=2500*mm) {
01483 part = 0; layer = 1; mat = 1;
01484 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01485 b[0] = -tb[1];
01486 b[1] = tb[0];
01487 b[2] = tb[2];
01488 ifend = true;
01489 }
01490 if(2210*mm<=mr[2]&&mr[2]<2240*mm&&1100*mm<mr[0]&&mr[0]<=2500*mm) {
01491 part = 0; layer = 2; mat = 0;
01492 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01493 b[0] = -tb[1];
01494 b[1] = tb[0];
01495 b[2] = tb[2];
01496 ifend = true;
01497 }
01498 if(2240*mm<=mr[2]&&mr[2]<2280*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
01499 part = 0; layer = 2; mat = 1;
01500 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01501 b[0] = -tb[1];
01502 b[1] = tb[0];
01503 b[2] = tb[2];
01504 ifend = true;
01505 }
01506 if(2280*mm<=mr[2]&&mr[2]<2310*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
01507 part = 0; layer = 3; mat = 0;
01508 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01509 b[0] = -tb[1];
01510 b[1] = tb[0];
01511 b[2] = tb[2];
01512 ifend = true;
01513 }
01514 if(2310*mm<=mr[2]&&mr[2]<2350*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
01515 part = 0; layer = 3; mat = 1;
01516 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01517 b[0] = -tb[1];
01518 b[1] = tb[0];
01519 b[2] = tb[2];
01520 ifend = true;
01521 }
01522 if(2350*mm<=mr[2]&&mr[2]<2380*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
01523 part = 0; layer = 4; mat = 0;
01524 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01525 b[0] = -tb[1];
01526 b[1] = tb[0];
01527 b[2] = tb[2];
01528 ifend = true;
01529 }
01530 if(2380*mm<=mr[2]&&mr[2]<2420*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
01531 part = 0; layer = 4; mat = 1;
01532 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01533 b[0] = -tb[1];
01534 b[1] = tb[0];
01535 b[2] = tb[2];
01536 ifend = true;
01537 }
01538 if(2420*mm<=mr[2]&&mr[2]<2470*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
01539 part = 0; layer = 5; mat = 0;
01540 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01541 b[0] = -tb[1];
01542 b[1] = tb[0];
01543 b[2] = tb[2];
01544 ifend = true;
01545 }
01546 if(2470*mm<=mr[2]&&mr[2]<2510*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
01547 part = 0; layer = 5; mat =1;
01548 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01549 b[0] = -tb[1];
01550 b[1] = tb[0];
01551 b[2] = tb[2];
01552 ifend = true;
01553 }
01554 if(2510*mm<=mr[2]&&mr[2]<2590*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
01555 part = 0; layer = 6; mat = 0;
01556 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01557 b[0] = -tb[1];
01558 b[1] = tb[0];
01559 b[2] = tb[2];
01560 ifend = true;
01561 }
01562 if(2590*mm<=mr[2]&&mr[2]<2630*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01563 part = 0; layer = 6; mat = 1;
01564 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01565 b[0] = -tb[1];
01566 b[1] = tb[0];
01567 b[2] = tb[2];
01568 ifend = true;
01569 }
01570 if(2630*mm<=mr[2]&&mr[2]<2710*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01571 part = 0; layer = 7; mat = 0;
01572 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01573 b[0] = -tb[1];
01574 b[1] = tb[0];
01575 b[2] = tb[2];
01576 ifend = true;
01577 }
01578 if(2710*mm<=mr[2]&&mr[2]<2750*mm&&1362*mm<=mr[0]&&mr[0]<=2500*mm) {
01579 part = 0; layer = 7; mat = 1;
01580 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01581 b[0] = -tb[1];
01582 b[1] = tb[0];
01583 b[2] = tb[2];
01584 ifend = true;
01585 }
01586 if(2750*mm<=mr[2]&&mr[2]<=2800*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01587 part = 0; layer = 8; mat = 0;
01588 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01589 b[0] = -tb[1];
01590 b[1] = tb[0];
01591 b[2] = tb[2];
01592 ifend = true;
01593 }
01594 }
01595 if(ifbar==true||ifend==true) {
01596 if( r.x() < 0. && r.y() >= 0. && r.z() > 0. ){
01597 b[0] = -b[0];
01598 }
01599 else if( r.x() <= 0. && r.y() < 0. && r.z() > 0. ){
01600 b[0] = -b[0];
01601 b[1] = -b[1];
01602 }
01603 else if( r.x() > 0. && r.y() < 0. && r.z() > 0. ){
01604 b[1] = -b[1];
01605 }
01606 else if( r.x() >= 0. && r.y() > 0. && r.z() < 0. ){
01607 b[0] = -b[0];
01608 b[1] = -b[1];
01609 }
01610 else if( r.x() < 0. && r.y() >= 0. && r.z() < 0. ){
01611 b[1] = -b[1];
01612 }
01613 else if( r.x() <= 0. && r.y() < 0. && r.z() < 0. ){
01614 b[0] = b[0];
01615 b[1] = b[1];
01616 }
01617 else if( r.x() > 0. && r.y() <= 0. && r.z() < 0. ){
01618 b[0] = -b[0];
01619 }
01620 }
01621 }
01622
01623 //reference frames defined by simulation and SCM are different.
01624 newb[0] = -b[0];
01625 newb[1] = b[1];
01626 newb[2] = -b[2];
01627 /* if(-1.2*m<r.z()&&r.z()<1.2*m&&0*m<=std::sqrt(r.x()*r.x()+r.y()*r.y())&&std::sqrt(r.x()*r.x()+r.y()*r.y())<0.9*m) {
01628 return StatusCode::SUCCESS;
01629 }
01630 else {
01631 newb[0] = newb[0] - 0.10*tesla;
01632 newb[1] = newb[1] + 0.10*tesla;
01633 newb[2] = newb[2] - 0.10*tesla;
01634 }*/
01635
01636 //newb[0] = b[0];
01637 //newb[1] = b[1];
01638 //newb[2] = b[2];
01639
01640 return StatusCode::SUCCESS;
01641 }
|
|
|
Finalise the service.
|
|
|
Finalise the service.
00445 {
00446 MsgStream log( msgSvc(), name() );
00447 //if(m_connect_run) delete m_connect_run;
00448 StatusCode status = Service::finalize();
00449
00450 if ( status.isSuccess() )
00451 log << MSG::INFO << "Service finalized successfully" << endreq;
00452 return status;
00453 }
|
|
|
Implements IMagneticFieldSvc. |
|
|
Implements IMagneticFieldSvc. 01678 {
01679 if(m_useDB == false) {
01680 if(m_runmode == 8 || m_runmode == 7) m_zfield = -0.9*tesla;
01681 else m_zfield = -1.0*tesla;
01682 }
01683 return m_zfield;
01684 }
|
|
|
|
|
|
00476 {
00477 MsgStream log( messageService(), name() );
00478 log << MSG::DEBUG << "handle: " << inc.type() << endreq;
00479 if ( inc.type() != "NewRun" ){
00480 return;
00481 }
00482 log << MSG::DEBUG << "Begin New Run" << endreq;
00483 if(m_useDB == true) {
00484 init_params();
00485 }
00486 }
|
|
|
Implements IMagneticFieldSvc. |
|
|
Implements IMagneticFieldSvc. 01686 {
01687 return m_ifRealField;
01688 }
|
|
|
|
|
|
00181 {
00182 MsgStream log(msgSvc(), name());
00183
00184 m_Q.clear();
00185 m_P.clear();
00186 m_Q_1.clear();
00187 m_P_1.clear();
00188 m_Q_2.clear();
00189 m_P_2.clear();
00190 m_P_TE.clear();
00191 m_Q_TE.clear();
00192
00193 if(m_gridDistance == 5) {
00194 m_Q.reserve(201250);
00195 m_P.reserve(201250);
00196 m_Q_1.reserve(201250);
00197 m_P_1.reserve(201250);
00198 m_Q_2.reserve(201250);
00199 m_P_2.reserve(201250);
00200 m_Q_TE.reserve(176608);
00201 m_P_TE.reserve(176608);
00202 }
00203 if(m_gridDistance == 10) {
00204 m_Q.reserve(27082);
00205 m_P.reserve(27082);
00206 m_Q_1.reserve(27082);
00207 m_P_1.reserve(27082);
00208 m_Q_2.reserve(27082);
00209 m_P_2.reserve(27082);
00210 m_Q_TE.reserve(23833);
00211 m_P_TE.reserve(23833);
00212 }
00213
00214 int runNo;
00215 SmartDataPtr<Event::EventHeader> evt(m_eventSvc,"/Event/EventHeader");
00216 if( evt ){
00217 runNo = evt -> runNumber();
00218 log << MSG::INFO <<"The runNumber of current event is "<<runNo<<endreq;
00219 }
00220 else {
00221 log << MSG::ERROR <<"Can not get EventHeader!"<<endreq;
00222 }
00223
00224 using FieldDBUtil::ConnectionDB;
00225 std::vector<double> current;
00226 current.clear();
00227 ConnectionDB::eRet e = m_connect_run->getReadSC_MagnetInfo(current,std::abs(runNo));
00228
00229 int ssm_curr = int (current[0]);
00230 double scql_curr = current[1];
00231 double scqr_curr = current[2];
00232 log << MSG::INFO << "SSM current: " << current[0] << " SCQL current: " << current[1] << " SCQR current: " << current[2] << " in Run " << runNo << endreq;
00233
00234 //int ssm_curr = 3369;
00235 //double scql_curr = 426.2;
00236 //double scqr_curr = 426.2;
00237 //for the energy less than 1.89GeV
00238 if((ssm_curr >= 3367) && (ssm_curr <= 3370) && ((scql_curr+scqr_curr)/2 < m_Cur_SCQ1_89)) {
00239 m_zfield = -1.0*tesla;
00240 if(getenv("MAGNETICFIELDROOT") != NULL) {
00241 m_filename = std::string(getenv( "MAGNETICFIELDROOT" ));
00242 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode2.dat");
00243 }
00244 log << MSG::INFO << "Select field map: " << m_filename << endreq;
00245 StatusCode status = parseFile();
00246 if ( !status.isSuccess() ) {
00247 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
00248 }
00249 m_Q_1 = m_Q;
00250 m_P_1 = m_P;
00251
00252 m_Q.clear();
00253 m_P.clear();
00254 if(m_gridDistance == 5) {
00255 m_Q.reserve(201250);
00256 m_P.reserve(201250);
00257 }
00258 if(m_gridDistance == 10) {
00259 m_Q.reserve(27082);
00260 m_P.reserve(27082);
00261 }
00262
00263 if(getenv("MAGNETICFIELDROOT") != NULL) {
00264 m_filename = std::string(getenv( "MAGNETICFIELDROOT" ));
00265 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode3.dat");
00266 }
00267 log << MSG::INFO << "Select field map: " << m_filename << endreq;
00268 status = parseFile();
00269 if ( !status.isSuccess() ) {
00270 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
00271 }
00272 m_Q_2 = m_Q;
00273 m_P_2 = m_P;
00274
00275 m_Q.clear();
00276 m_P.clear();
00277 if(m_gridDistance == 5) {
00278 m_Q.reserve(201250);
00279 m_P.reserve(201250);
00280 }
00281 if(m_gridDistance == 10) {
00282 m_Q.reserve(27082);
00283 m_P.reserve(27082);
00284 }
00285 //check
00286 if(m_Q_2.size() != m_Q_1.size()) log << MSG::FATAL << "The two field maps used in this run are wrong!!!" << endreq;
00287
00288 for(int iQ = 0; iQ < m_Q_1.size(); iQ++) {
00289 double fieldvalue = (m_Q_1[iQ] - m_Q_2[iQ])/(m_Cur_SCQ1_55 - m_Cur_SCQ1_89)*((scql_curr+scqr_curr)/2 - m_Cur_SCQ1_89) + m_Q_2[iQ];
00290 m_Q.push_back(fieldvalue);
00291 m_P.push_back(m_P_2[iQ]);
00292 }
00293 }
00294 //for the energy greater than 1.89GeV
00295 if((ssm_curr >= 3367) && (ssm_curr <= 3370) && ((scql_curr+scqr_curr)/2 >= m_Cur_SCQ1_89)) {
00296 m_zfield = -1.0*tesla;
00297 if(getenv("MAGNETICFIELDROOT") != NULL) {
00298 m_filename = std::string(getenv( "MAGNETICFIELDROOT" ));
00299 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode3.dat");
00300 }
00301 log << MSG::INFO << "Select field map: " << m_filename << endreq;
00302 StatusCode status = parseFile();
00303 if ( !status.isSuccess() ) {
00304 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
00305 }
00306 m_Q_1 = m_Q;
00307 m_P_1 = m_P;
00308
00309 m_Q.clear();
00310 m_P.clear();
00311 if(m_gridDistance == 5) {
00312 m_Q.reserve(201250);
00313 m_P.reserve(201250);
00314 }
00315 if(m_gridDistance == 10) {
00316 m_Q.reserve(27082);
00317 m_P.reserve(27082);
00318 }
00319
00320 if(getenv("MAGNETICFIELDROOT") != NULL) {
00321 m_filename = std::string(getenv( "MAGNETICFIELDROOT" ));
00322 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode4.dat");
00323 }
00324 log << MSG::INFO << "Select field map: " << m_filename << endreq;
00325 status = parseFile();
00326 if ( !status.isSuccess() ) {
00327 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
00328 }
00329 m_Q_2 = m_Q;
00330 m_P_2 = m_P;
00331
00332 m_Q.clear();
00333 m_P.clear();
00334 if(m_gridDistance == 5) {
00335 m_Q.reserve(201250);
00336 m_P.reserve(201250);
00337 }
00338 if(m_gridDistance == 10) {
00339 m_Q.reserve(27082);
00340 m_P.reserve(27082);
00341 }
00342 //check
00343 if(m_Q_2.size() != m_Q_1.size()) log << MSG::FATAL << "The two field maps used in this run are wrong!!!" << endreq;
00344
00345 for(int iQ = 0; iQ < m_Q_1.size(); iQ++) {
00346 double fieldvalue = (m_Q_1[iQ] - m_Q_2[iQ])/(m_Cur_SCQ1_89 - m_Cur_SCQ2_10)*((scql_curr+scqr_curr)/2 - m_Cur_SCQ2_10) + m_Q_2[iQ];
00347 m_Q.push_back(fieldvalue);
00348 m_P.push_back(m_P_2[iQ]);
00349 }
00350 }
00351 if((ssm_curr >= 3033) && (ssm_curr <= 3035)) {
00352 m_zfield = -0.9*tesla;
00353 if(getenv("MAGNETICFIELDROOT") != NULL) {
00354 m_filename = std::string(getenv( "MAGNETICFIELDROOT" ));
00355 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode7.dat");
00356 }
00357 log << MSG::INFO << "Select field map: " << m_filename << endreq;
00358 StatusCode status = parseFile();
00359 if ( !status.isSuccess() ) {
00360 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
00361 }
00362 m_Q_1 = m_Q;
00363 m_P_1 = m_P;
00364
00365 m_Q.clear();
00366 m_P.clear();
00367 if(m_gridDistance == 5) {
00368 m_Q.reserve(201250);
00369 m_P.reserve(201250);
00370 }
00371 if(m_gridDistance == 10) {
00372 m_Q.reserve(27082);
00373 m_P.reserve(27082);
00374 }
00375
00376 if(getenv("MAGNETICFIELDROOT") != NULL) {
00377 m_filename = std::string(getenv( "MAGNETICFIELDROOT" ));
00378 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode8.dat");
00379 }
00380 log << MSG::INFO << "Select field map: " << m_filename << endreq;
00381 status = parseFile();
00382 if ( !status.isSuccess() ) {
00383 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
00384 }
00385 m_Q_2 = m_Q;
00386 m_P_2 = m_P;
00387
00388 m_Q.clear();
00389 m_P.clear();
00390 if(m_gridDistance == 5) {
00391 m_Q.reserve(201250);
00392 m_P.reserve(201250);
00393 }
00394 if(m_gridDistance == 10) {
00395 m_Q.reserve(27082);
00396 m_P.reserve(27082);
00397 }
00398 //check
00399 if(m_Q_2.size() != m_Q_1.size()) log << MSG::FATAL << "The two field maps used in this run are wrong!!!" << endreq;
00400
00401 for(int iQ = 0; iQ < m_Q_1.size(); iQ++) {
00402 double fieldvalue = (m_Q_1[iQ] - m_Q_2[iQ])/(m_Cur_SCQ1_55 - m_Cur_SCQ1_89)*((scql_curr+scqr_curr)/2 - m_Cur_SCQ1_89) + m_Q_2[iQ];
00403 m_Q.push_back(fieldvalue);
00404 m_P.push_back(m_P_2[iQ]);
00405 }
00406 }
00407
00408 if((!((ssm_curr >= 3367) && (ssm_curr <= 3370)) && !((ssm_curr >= 3033) && (ssm_curr <= 3035))) || scql_curr == 0 || scqr_curr == 0) {
00409 log << MSG::FATAL << "The current of run " << runNo << " is abnormal in database, please check it! " << endreq;
00410 log << MSG::FATAL << "The current of SSM is " << ssm_curr << ", SCQR is " << scqr_curr << ", SCQL is " << scql_curr << endreq;
00411 exit(1);
00412 }
00413 if(m_Q.size() == 0) {
00414 log << MSG::FATAL << "This size of field map vector is ZERO, please check the current of run " << runNo << " in database!" << endreq;
00415 exit(1);
00416 }
00417
00418 m_filename_TE = std::string(getenv( "MAGNETICFIELDROOT" ));
00419 if(m_gridDistance == 10) m_filename_TE += std::string( "/dat/TEMap10cm.dat");
00420 if(m_gridDistance == 5) m_filename_TE += std::string( "/dat/TEMap5cm.dat");
00421 log << MSG::INFO << "Select field map: " << m_filename_TE << endreq;
00422 StatusCode status = parseFile_TE();
00423 if ( !status.isSuccess() ) {
00424 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
00425 }
00426
00427 for (int iC = 0; iC< 2; ++iC ){
00428 m_min_FL[iC] = -90.0*cm;
00429 m_max_FL[iC] = m_min_FL[iC]+( m_Nxyz[iC]-1 )* m_Dxyz[iC];
00430 //for tof and emc
00431 m_min_FL_TE[iC] = 0.0*cm;
00432 m_max_FL_TE[iC] = m_min_FL_TE[iC]+( m_Nxyz_TE[iC]-1 )* m_Dxyz_TE[iC];
00433 } // iC
00434 m_min_FL[2] = -120.0*cm;
00435 m_max_FL[2] = m_min_FL[2]+( m_Nxyz[2]-1 )* m_Dxyz[2];
00436 //for tof and emc
00437 m_min_FL_TE[2] = 0.0*cm;
00438 m_max_FL_TE[2] = m_min_FL_TE[2]+( m_Nxyz_TE[2]-1 )* m_Dxyz_TE[2];
00439 }
|
|
|
Initialise the service (Inherited Service overrides).
|
|
|
Initialise the service (Inherited Service overrides).
00069 {
00070 MsgStream log(msgSvc(), name());
00071 StatusCode status = Service::initialize();
00072 // Get the references to the services that are needed by the ApplicationMgr itself
00073 IIncidentSvc* incsvc;
00074 status = service("IncidentSvc", incsvc);
00075 int priority = 100;
00076 if( status.isSuccess() ){
00077 incsvc -> addListener(this, "NewRun", priority);
00078 }
00079
00080 status = serviceLocator()->service("EventDataSvc", m_eventSvc, true);
00081 if (status.isFailure() ) {
00082 log << MSG::ERROR << "Unable to find EventDataSvc " << endreq;
00083 return status;
00084 }
00085
00086 if(m_useDB == false) {
00087 if(m_gridDistance == 5) {
00088 m_Q.reserve(201250);
00089 m_P.reserve(201250);
00090 m_Q_TE.reserve(176608);
00091 m_P_TE.reserve(176608);
00092 }
00093 if(m_gridDistance == 10) {
00094 m_Q.reserve(27082);
00095 m_P.reserve(27082);
00096 m_Q_TE.reserve(23833);
00097 m_P_TE.reserve(23833);
00098 }
00099
00100 if(getenv("MAGNETICFIELDROOT") != NULL) {
00101 m_filename = std::string(getenv( "MAGNETICFIELDROOT" ));
00102 m_filename_TE = std::string(getenv( "MAGNETICFIELDROOT" ));
00103 if(m_gridDistance == 10) {
00104 m_filename_TE += std::string("/dat/TEMap10cm.dat");
00105 if(m_runmode == 2)
00106 m_filename += std::string( "/dat/2008-04-03BESIII-field-Mode2.dat");
00107 else if(m_runmode == 4)
00108 m_filename += std::string( "/dat/2008-04-03BESIII-field-Mode4.dat");
00109 else if(m_runmode == 6)
00110 m_filename += std::string( "/dat/2008-04-03BESIII-field-Mode6.dat");
00111 else if(m_runmode == 7)
00112 m_filename += std::string( "/dat/2008-04-03BESIII-field-Mode7.dat");
00113 else {
00114 m_filename += std::string( "/dat/2007-08-28BESIII-field.dat");
00115 std::cout<<"WARNING: YOU ARE USING OLD FIELD MAP!!!!"<<std::endl;
00116 }
00117 }
00118 if(m_gridDistance == 5) {
00119 m_filename_TE += std::string("/dat/TEMap5cm.dat");
00120 if(m_runmode == 1)
00121 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode1.dat");
00122 else if(m_runmode == 2)
00123 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode2.dat");
00124 else if(m_runmode == 3)
00125 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode3.dat");
00126 else if(m_runmode == 4)
00127 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode4.dat");
00128 else if(m_runmode == 5)
00129 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode5.dat");
00130 else if(m_runmode == 6)
00131 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode6.dat");
00132 else if(m_runmode == 7)
00133 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode7.dat");
00134 else if(m_runmode == 8)
00135 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode8.dat");
00136 else {
00137 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode2.dat");
00138 std::cout<<"WARNING: NO RUN MODE, YOU ARE USING DEFAULT FIELD MAP!!!!"<<std::endl;
00139 }
00140 }
00141 cout<<"*** Read field map: "<<m_filename<<endl;
00142 cout<<"*** Read field map: "<<m_filename_TE<<endl;
00143 }
00144 else {
00145 std::cerr<<"Couldn't find MAGNETICFIELDROOT"<<std::endl;
00146 }
00147
00148 status = parseFile();
00149 status = parseFile_TE();
00150 if ( status.isSuccess() ) {
00151 log << MSG::DEBUG << "Magnetic field parsed successfully" << endreq;
00152
00153 for (int iC = 0; iC< 2; ++iC ){
00154 m_min_FL[iC] = -90.0*cm;
00155 m_max_FL[iC] = m_min_FL[iC]+( m_Nxyz[iC]-1 )* m_Dxyz[iC];
00156 //for tof and emc
00157 m_min_FL_TE[iC] = 0.0*cm;
00158 m_max_FL_TE[iC] = m_min_FL_TE[iC]+( m_Nxyz_TE[iC]-1 )* m_Dxyz_TE[iC];
00159 } // iC
00160 m_min_FL[2] = -120.0*cm;
00161 m_max_FL[2] = m_min_FL[2]+( m_Nxyz[2]-1 )* m_Dxyz[2];
00162 //for tof and emc
00163 m_min_FL_TE[2] = 0.0*cm;
00164 m_max_FL_TE[2] = m_min_FL_TE[2]+( m_Nxyz_TE[2]-1 )* m_Dxyz_TE[2];
00165 }
00166 else {
00167 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
00168 }
00169 }
00170 else {
00171 m_connect_run = new FieldDBUtil::ConnectionDB();
00172 if (!m_connect_run) {
00173 log << MSG::ERROR << "Could not open connection to database" << endreq;
00174 }
00175 }
00176
00177 return status;
00178 }
|
|
|
Retrieve interface ID.
00033 { return IID_IMagneticFieldSvc; }
|
|
|
Retrieve interface ID.
00033 { return IID_IMagneticFieldSvc; }
|
|
|
Reads the field map from file.
|
|
|
Reads the field map from file.
00492 {
00493 StatusCode sc = StatusCode::FAILURE;
00494
00495 MsgStream log( msgSvc(), name() );
00496 char line[ 255 ];
00497 std::ifstream infile( m_filename.c_str() );
00498
00499 if ( infile ) {
00500 sc = StatusCode::SUCCESS;
00501
00502 // Skip the header till PARAMETER
00503 do{
00504 infile.getline( line, 255 );
00505 } while( line[0] != 'P' );
00506
00507 // Get the PARAMETER
00508 std::string sPar[2];
00509 char* token = strtok( line, " " );
00510 int ip = 0;
00511 do{
00512 if ( token ) { sPar[ip] = token; token = strtok( NULL, " " );}
00513 else continue;
00514 ip++;
00515 } while ( token != NULL );
00516 long int npar = atoi( sPar[1].c_str() );
00517
00518 // Skip the header till GEOMETRY
00519 do{
00520 infile.getline( line, 255 );
00521 } while( line[0] != 'G' );
00522
00523 // Skip any comment before GEOMETRY
00524 do{
00525 infile.getline( line, 255 );
00526 } while( line[0] != '#' );
00527
00528 // Get the GEOMETRY
00529 infile.getline( line, 255 );
00530 std::string sGeom[7];
00531 token = strtok( line, " " );
00532 int ig = 0;
00533 do{
00534 if ( token ) { sGeom[ig] = token; token = strtok( NULL, " " );}
00535 else continue;
00536 ig++;
00537 } while (token != NULL);
00538
00539 // Grid dimensions are given in cm in CDF file. Convert to CLHEP units
00540 m_Dxyz[0] = atof( sGeom[0].c_str() ) * cm;
00541 m_Dxyz[1] = atof( sGeom[1].c_str() ) * cm;
00542 m_Dxyz[2] = atof( sGeom[2].c_str() ) * cm;
00543 m_Nxyz[0] = atoi( sGeom[3].c_str() );
00544 m_Nxyz[1] = atoi( sGeom[4].c_str() );
00545 m_Nxyz[2] = atoi( sGeom[5].c_str() );
00546 m_zOffSet = atof( sGeom[6].c_str() ) * cm;
00547 // Number of lines with data to be read
00548 long int nlines = ( npar - 7 ) / 3;
00549
00550 // Check number of lines with data read in the loop
00551 long int ncheck = 0;
00552
00553 // Skip comments and fill a vector of magnetic components for the
00554 // x, y and z positions given in GEOMETRY
00555
00556 while( infile ) {
00557 // parse each line of the file,
00558 // comment lines begin with '#' in the cdf file
00559 infile.getline( line, 255 );
00560 if ( line[0] == '#' ) continue;
00561 std::string gridx, gridy, gridz, sFx, sFy, sFz;
00562 char* token = strtok( line, " " );
00563 if ( token ) { gridx = token; token = strtok( NULL, " " );} else continue;
00564 if ( token ) { gridy = token; token = strtok( NULL, " " );} else continue;
00565 if ( token ) { gridz = token; token = strtok( NULL, " " );} else continue;
00566 if ( token ) { sFx = token; token = strtok( NULL, " " );} else continue;
00567 if ( token ) { sFy = token; token = strtok( NULL, " " );} else continue;
00568 if ( token ) { sFz = token; token = strtok( NULL, " " );} else continue;
00569 if ( token != NULL ) continue;
00570 //Grid position
00571 double gx = atof( gridx.c_str() ) * m;
00572 double gy = atof( gridy.c_str() ) * m;
00573 double gz = atof( gridz.c_str() ) * m;
00574 // Field values are given in tesla in CDF file. Convert to CLHEP units
00575 double fx = atof( sFx.c_str() ) * tesla;
00576 double fy = atof( sFy.c_str() ) * tesla;
00577 double fz = atof( sFz.c_str() ) * tesla;
00578 //for debug
00579 if(m_outlevel == 0) {
00580 log << MSG::DEBUG << "grid x, y, z is "<< gx << ", " << gy << ", " << gz << endreq;
00581 log << MSG::DEBUG << "field x, y, z is "<< fx << ", " << fy << ", " << fz << endreq;
00582 }
00583 //Fill the postion to the vector
00584 m_P.push_back( gx );
00585 m_P.push_back( gy );
00586 m_P.push_back( gz );
00587 // Add the magnetic field components of each point to
00588 // sequentialy in a vector
00589 m_Q.push_back( fx );
00590 m_Q.push_back( fy );
00591 m_Q.push_back( fz );
00592 // counts after reading and filling to match the number of lines
00593 ncheck++;
00594 }
00595 infile.close();
00596 if ( nlines != ncheck) {
00597 log << MSG::ERROR << "nline is " << nlines << "; ncheck is " << ncheck << " Number of points in field map does not match!"
00598 << endreq;
00599 return StatusCode::FAILURE;
00600 }
00601 }
00602 else {
00603 log << MSG::ERROR << "Unable to open magnetic field file : "
00604 << m_filename << endreq;
00605 }
00606
00607 return sc;
00608 }
|
|
|
Reads the field map from file.
|
|
|
Reads the field map from file.
00615 {
00616 StatusCode sc = StatusCode::FAILURE;
00617
00618 MsgStream log( msgSvc(), name() );
00619 char line[ 255 ];
00620 std::ifstream infile( m_filename_TE.c_str() );
00621
00622 if ( infile ) {
00623 sc = StatusCode::SUCCESS;
00624
00625 // Skip the header till PARAMETER
00626 do{
00627 infile.getline( line, 255 );
00628 } while( line[0] != 'P' );
00629
00630 // Get the PARAMETER
00631 std::string sPar[2];
00632 char* token = strtok( line, " " );
00633 int ip = 0;
00634 do{
00635 if ( token ) { sPar[ip] = token; token = strtok( NULL, " " );}
00636 else continue;
00637 ip++;
00638 } while ( token != NULL );
00639 long int npar = atoi( sPar[1].c_str() );
00640
00641 // Skip the header till GEOMETRY
00642 do{
00643 infile.getline( line, 255 );
00644 } while( line[0] != 'G' );
00645
00646 // Skip any comment before GEOMETRY
00647 do{
00648 infile.getline( line, 255 );
00649 } while( line[0] != '#' );
00650
00651 // Get the GEOMETRY
00652 infile.getline( line, 255 );
00653 std::string sGeom[7];
00654 token = strtok( line, " " );
00655 int ig = 0;
00656 do{
00657 if ( token ) { sGeom[ig] = token; token = strtok( NULL, " " );}
00658 else continue;
00659 ig++;
00660 } while (token != NULL);
00661
00662 // Grid dimensions are given in cm in CDF file. Convert to CLHEP units
00663 m_Dxyz_TE[0] = atof( sGeom[0].c_str() ) * cm;
00664 m_Dxyz_TE[1] = atof( sGeom[1].c_str() ) * cm;
00665 m_Dxyz_TE[2] = atof( sGeom[2].c_str() ) * cm;
00666 m_Nxyz_TE[0] = atoi( sGeom[3].c_str() );
00667 m_Nxyz_TE[1] = atoi( sGeom[4].c_str() );
00668 m_Nxyz_TE[2] = atoi( sGeom[5].c_str() );
00669 m_zOffSet_TE = atof( sGeom[6].c_str() ) * cm;
00670 // Number of lines with data to be read
00671 long int nlines = ( npar - 7 ) / 3;
00672
00673 // Check number of lines with data read in the loop
00674 long int ncheck = 0;
00675
00676 // Skip comments and fill a vector of magnetic components for the
00677 // x, y and z positions given in GEOMETRY
00678
00679 while( infile ) {
00680 // parse each line of the file,
00681 // comment lines begin with '#' in the cdf file
00682 infile.getline( line, 255 );
00683 if ( line[0] == '#' ) continue;
00684 std::string gridx, gridy, gridz, sFx, sFy, sFz;
00685 char* token = strtok( line, " " );
00686 if ( token ) { gridx = token; token = strtok( NULL, " " );} else continue;
00687 if ( token ) { gridy = token; token = strtok( NULL, " " );} else continue;
00688 if ( token ) { gridz = token; token = strtok( NULL, " " );} else continue;
00689 if ( token ) { sFx = token; token = strtok( NULL, " " );} else continue;
00690 if ( token ) { sFy = token; token = strtok( NULL, " " );} else continue;
00691 if ( token ) { sFz = token; token = strtok( NULL, " " );} else continue;
00692 if ( token != NULL ) continue;
00693 //Grid position
00694 double gx = atof( gridx.c_str() ) * m;
00695 double gy = atof( gridy.c_str() ) * m;
00696 double gz = atof( gridz.c_str() ) * m;
00697 // Field values are given in tesla in CDF file. Convert to CLHEP units
00698 double fx = atof( sFx.c_str() ) * tesla;
00699 double fy = atof( sFy.c_str() ) * tesla;
00700 double fz = atof( sFz.c_str() ) * tesla;
00701 //for debug
00702 if(m_outlevel == 0) {
00703 log << MSG::DEBUG << "grid x, y, z is "<< gx << ", " << gy << ", " << gz << endreq;
00704 log << MSG::DEBUG << "field x, y, z is "<< fx << ", " << fy << ", " << fz << endreq;
00705 }
00706 //Fill the postion to the vector
00707 m_P_TE.push_back( gx );
00708 m_P_TE.push_back( gy );
00709 m_P_TE.push_back( gz );
00710 // Add the magnetic field components of each point to
00711 // sequentialy in a vector
00712 m_Q_TE.push_back( fx );
00713 m_Q_TE.push_back( fy );
00714 m_Q_TE.push_back( fz );
00715 // counts after reading and filling to match the number of lines
00716 ncheck++;
00717 }
00718 infile.close();
00719 if ( nlines != ncheck) {
00720 log << MSG::ERROR << "nline is " << nlines << "; ncheck is " << ncheck << " Number of points in field map does not match!"
00721 << endreq;
00722 return StatusCode::FAILURE;
00723 }
00724 }
00725 else {
00726 log << MSG::ERROR << "Unable to open magnetic field file : "
00727 << m_filename_TE << endreq;
00728 }
00729
00730 return sc;
00731 }
|
|
||||||||||||
|
Query the available interfaces.
|
|
||||||||||||
|
Query the available interfaces.
00460 {
00461 StatusCode sc = StatusCode::FAILURE;
00462 if ( ppvInterface ) {
00463 *ppvInterface = 0;
00464
00465 if ( riid == IID_IMagneticFieldSvc ) {
00466 *ppvInterface = static_cast<IMagneticFieldSvc*>(this);
00467 sc = StatusCode::SUCCESS;
00468 addRef();
00469 }
00470 else
00471 sc = Service::queryInterface( riid, ppvInterface );
00472 }
00473 return sc;
00474 }
|
|
|
Service type.
00064 { return IID_IMagneticFieldSvc; };
|
|
|
Service type.
00064 { return IID_IMagneticFieldSvc; };
|
|
||||||||||||
|
|
|
||||||||||||
|
|
|
||||||||||||
|
01645 {
01646 if(m_runmode == 8 || m_runmode == 7) {
01647 if(-1.5*m<=r.z()&&r.z()<=1.5*m&&0*m<=std::sqrt(r.x()*r.x()+r.y()*r.y())&&std::sqrt(r.x()*r.x()+r.y()*r.y())<=1.5*m)
01648 {
01649 b[0]=0.0;
01650 b[1]=0.0;
01651 b[2]=-0.9*tesla;
01652 }
01653 else
01654 {
01655 b[0]=0.0;
01656 b[1]=0.0;
01657 b[2]=0.0;
01658 }
01659 }
01660 else {
01661 if(-1.5*m<=r.z()&&r.z()<=1.5*m&&0*m<=std::sqrt(r.x()*r.x()+r.y()*r.y())&&std::sqrt(r.x()*r.x()+r.y()*r.y())<=1.5*m)
01662 {
01663 b[0]=0.0;
01664 b[1]=0.0;
01665 b[2]=-1.0*tesla;
01666 }
01667 else
01668 {
01669 b[0]=0.0;
01670 b[1]=0.0;
01671 b[2]=0.0;
01672 }
01673 }
01674 return StatusCode::SUCCESS;
01675 }
|
|
|
Allow SvcFactory to instantiate the service.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Steps in x, y and z.
|
|
|
Steps in x, y and z.
|
|
|
|
|
|
|
|
|
Magnetic field file name.
|
|
|
Magnetic field file name.
|
|
|
grid distance of field map
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Number of steps in x, y and z.
|
|
|
Number of steps in x, y and z.
|
|
|
|
|
|
Grid position.
|
|
|
Grid position.
|
|
|
Grid position.
|
|
|
Grid position.
|
|
|
Grid position.
|
|
|
Grid position.
|
|
|
Grid position.
|
|
|
Grid position.
|
|
|
Field vector.
|
|
|
Field vector.
|
|
|
Field vector.
|
|
|
Field vector.
|
|
|
Field vector.
|
|
|
Field vector.
|
|
|
Field vector.
|
|
|
Field vector.
|
|
|
Run mode.
|
|
|
|
|
|
|
|
|
The z offset.
|
|
|
The z offset.
|
1.3.9.1