00001
00002 #ifndef BEAN
00003 #include "GaudiKernel/AlgFactory.h"
00004 #include "GaudiKernel/SvcFactory.h"
00005 #include "GaudiKernel/ISvcLocator.h"
00006 #include "GaudiKernel/MsgStream.h"
00007 #include "GaudiKernel/Incident.h"
00008 #include "GaudiKernel/IIncidentSvc.h"
00009 #include "GaudiKernel/IDataProviderSvc.h"
00010 #include "GaudiKernel/DataObject.h"
00011 #include "GaudiKernel/SmartDataPtr.h"
00012
00013 #include "MagneticField/IMagneticFieldSvc.h"
00014 #endif
00015
00016 #include "MagneticField/MagneticFieldSvc.h"
00017 #include "MagneticField/MucMagneticField.h"
00018
00019 #include "CLHEP/Units/PhysicalConstants.h"
00020
00021 #ifndef BEAN
00022 #include "EventModel/EventModel.h"
00023 #include "EventModel/EventHeader.h"
00024 #endif
00025
00026 #include <cstring>
00027 #include <cstdlib>
00028 #include <fstream>
00029 #include <iostream>
00030 using namespace std;
00031 using namespace CLHEP;
00032
00038 #ifndef BEAN
00039
00040
00041
00042
00043
00044
00045
00046
00047 MagneticFieldSvc::MagneticFieldSvc( const std::string& name,
00048 ISvcLocator* svc ) : Service( name, svc )
00049 {
00050 declareProperty( "TurnOffField", m_turnOffField = false);
00051 declareProperty( "FieldMapFile", m_filename );
00052 declareProperty( "GridDistance", m_gridDistance = 5);
00053 declareProperty( "RunMode", m_runmode = 2);
00054 declareProperty( "IfRealField", m_ifRealField = true);
00055 declareProperty( "OutLevel", m_outlevel = 1);
00056 declareProperty( "Scale", m_scale = 1.0);
00057 declareProperty( "UniField", m_uniField = false);
00058
00059 declareProperty( "Cur_SCQ1_55", m_Cur_SCQ1_55 = 349.4);
00060 declareProperty( "Cur_SCQ1_89", m_Cur_SCQ1_89 = 426.2);
00061 declareProperty( "Cur_SCQ2_10", m_Cur_SCQ2_10 = 474.2);
00062
00063 declareProperty( "UseDBFlag", m_useDB = true);
00064
00065 m_Mucfield = new MucMagneticField();
00066 if(!m_Mucfield) cout<<"error: can not get MucMagneticField pointer"<<endl;
00067
00068 m_zfield = -1.0*tesla;
00069
00070 if(getenv("MAGNETICFIELDROOT") != NULL) {
00071 path = std::string(getenv( "MAGNETICFIELDROOT" ));
00072 } else {
00073 std::cerr<<"Couldn't find MAGNETICFIELDROOT"<<std::endl;
00074 }
00075
00076 }
00077
00078 #else
00079
00080 MagneticFieldSvc* MagneticFieldSvc::m_field = 0;
00081
00082
00083
00084
00085 MagneticFieldSvc::MagneticFieldSvc()
00086 {
00087
00088 m_gridDistance = 5;
00089 m_runmode = 2;
00090 m_ifRealField = true;
00091 m_outlevel = 1;
00092
00093 m_Cur_SCQ1_55 = 349.4;
00094 m_Cur_SCQ1_89 = 426.2;
00095 m_Cur_SCQ2_10 = 474.2;
00096
00097 m_useDB = true;
00098
00099 path = "";
00100 m_Mucfield = 0;
00101
00102 m_zfield = -1.0*tesla;
00103 }
00104 #endif
00105
00106
00107
00108
00109 MagneticFieldSvc::~MagneticFieldSvc()
00110 {
00111 if(m_Mucfield) delete m_Mucfield;
00112 }
00113
00114
00115
00116
00117 #ifdef BEAN
00118 bool MagneticFieldSvc::init_mucMagneticField()
00119 {
00120
00121 if( m_Mucfield ) {
00122 if( m_Mucfield->getPath() == path ) return true;
00123 delete m_Mucfield;
00124 }
00125
00126 m_Mucfield = new(std::nothrow) MucMagneticField(path);
00127 if( !m_Mucfield ) {
00128 cout<<"error: can not get MucMagneticField pointer"<<endl;
00129 return false;
00130 }
00131
00132 return true;
00133 }
00134 #endif
00135
00136 StatusCode MagneticFieldSvc::initialize()
00137 {
00138 #ifndef BEAN
00139 MsgStream log(msgSvc(), name());
00140 StatusCode status = Service::initialize();
00141
00142 IIncidentSvc* incsvc;
00143 status = service("IncidentSvc", incsvc);
00144 int priority = 100;
00145 if( status.isSuccess() ){
00146 incsvc -> addListener(this, "NewRun", priority);
00147 }
00148
00149 status = serviceLocator()->service("EventDataSvc", m_eventSvc, true);
00150 if (status.isFailure() ) {
00151 log << MSG::ERROR << "Unable to find EventDataSvc " << endreq;
00152 return status;
00153 }
00154 #else
00155 if( !init_mucMagneticField() ) return false;
00156 StatusCode status = true;
00157 #endif
00158
00159 if(m_useDB == false) {
00160 if(m_gridDistance == 5) {
00161 m_Q.reserve(201250);
00162 m_P.reserve(201250);
00163 m_Q_TE.reserve(176608);
00164 m_P_TE.reserve(176608);
00165 }
00166 if(m_gridDistance == 10) {
00167 m_Q.reserve(27082);
00168 m_P.reserve(27082);
00169 m_Q_TE.reserve(23833);
00170 m_P_TE.reserve(23833);
00171 }
00172
00173 m_filename = path;
00174 m_filename_TE = path;
00175 if(m_gridDistance == 10) {
00176 m_filename_TE += std::string("/dat/TEMap10cm.dat");
00177 if(m_runmode == 2)
00178 m_filename += std::string( "/dat/2008-04-03BESIII-field-Mode2.dat");
00179 else if(m_runmode == 4)
00180 m_filename += std::string( "/dat/2008-04-03BESIII-field-Mode4.dat");
00181 else if(m_runmode == 6)
00182 m_filename += std::string( "/dat/2008-04-03BESIII-field-Mode6.dat");
00183 else if(m_runmode == 7)
00184 m_filename += std::string( "/dat/2008-04-03BESIII-field-Mode7.dat");
00185 else {
00186 m_filename += std::string( "/dat/2007-08-28BESIII-field.dat");
00187 std::cout<<"WARNING: YOU ARE USING OLD FIELD MAP!!!!"<<std::endl;
00188 }
00189 }
00190 if(m_gridDistance == 5) {
00191 m_filename_TE += std::string("/dat/TEMap5cm.dat");
00192 if(m_runmode == 1)
00193 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode1.dat");
00194 else if(m_runmode == 2)
00195 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode2.dat");
00196 else if(m_runmode == 3)
00197 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode3.dat");
00198 else if(m_runmode == 4)
00199 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode4.dat");
00200 else if(m_runmode == 5)
00201 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode5.dat");
00202 else if(m_runmode == 6)
00203 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode6.dat");
00204 else if(m_runmode == 7)
00205 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode7.dat");
00206 else if(m_runmode == 8)
00207 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode8.dat");
00208 else {
00209 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode2.dat");
00210 std::cout<<"WARNING: NO RUN MODE, YOU ARE USING DEFAULT FIELD MAP!!!!"<<std::endl;
00211 }
00212 }
00213 cout<<"*** Read field map: "<<m_filename<<endl;
00214 cout<<"*** Read field map: "<<m_filename_TE<<endl;
00215
00216 status = parseFile();
00217 status = parseFile_TE();
00218
00219 #ifndef BEAN
00220 if ( status.isSuccess() ) {
00221 log << MSG::DEBUG << "Magnetic field parsed successfully" << endreq;
00222 #else
00223 if ( status ) {
00224 cout << "Magnetic field parsed successfully" << endl;
00225 #endif
00226
00227 for (int iC = 0; iC< 2; ++iC ){
00228 m_min_FL[iC] = -90.0*cm;
00229 m_max_FL[iC] = m_min_FL[iC]+( m_Nxyz[iC]-1 )* m_Dxyz[iC];
00230
00231 m_min_FL_TE[iC] = 0.0*cm;
00232 m_max_FL_TE[iC] = m_min_FL_TE[iC]+( m_Nxyz_TE[iC]-1 )* m_Dxyz_TE[iC];
00233 }
00234 m_min_FL[2] = -120.0*cm;
00235 m_max_FL[2] = m_min_FL[2]+( m_Nxyz[2]-1 )* m_Dxyz[2];
00236
00237 m_min_FL_TE[2] = 0.0*cm;
00238 m_max_FL_TE[2] = m_min_FL_TE[2]+( m_Nxyz_TE[2]-1 )* m_Dxyz_TE[2];
00239 }
00240 else {
00241 #ifndef BEAN
00242 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
00243 #else
00244 cout << "Magnetic field parse failled" << endl;
00245 #endif
00246 }
00247 }
00248 else {
00249 m_connect_run = new FieldDBUtil::ConnectionDB();
00250 #ifndef BEAN
00251 if (!m_connect_run) {
00252 log << MSG::ERROR << "Could not open connection to database" << endreq;
00253 }
00254 #endif
00255 }
00256
00257 return status;
00258 }
00259
00260 #ifndef BEAN
00261 void MagneticFieldSvc::init_params()
00262 {
00263 MsgStream log(msgSvc(), name());
00264 #else
00265 void MagneticFieldSvc::init_params(int run)
00266 {
00267 if( !init_mucMagneticField() ) {
00268 cerr << " STOP! " << endl;
00269 exit(1);
00270 }
00271 #endif
00272
00273 m_Q.clear();
00274 m_P.clear();
00275 m_Q_1.clear();
00276 m_P_1.clear();
00277 m_Q_2.clear();
00278 m_P_2.clear();
00279 m_P_TE.clear();
00280 m_Q_TE.clear();
00281
00282 if(m_gridDistance == 5) {
00283 m_Q.reserve(201250);
00284 m_P.reserve(201250);
00285 m_Q_1.reserve(201250);
00286 m_P_1.reserve(201250);
00287 m_Q_2.reserve(201250);
00288 m_P_2.reserve(201250);
00289 m_Q_TE.reserve(176608);
00290 m_P_TE.reserve(176608);
00291 }
00292 if(m_gridDistance == 10) {
00293 m_Q.reserve(27082);
00294 m_P.reserve(27082);
00295 m_Q_1.reserve(27082);
00296 m_P_1.reserve(27082);
00297 m_Q_2.reserve(27082);
00298 m_P_2.reserve(27082);
00299 m_Q_TE.reserve(23833);
00300 m_P_TE.reserve(23833);
00301 }
00302
00303 #ifndef BEAN
00304 int runNo;
00305 SmartDataPtr<Event::EventHeader> evt(m_eventSvc,"/Event/EventHeader");
00306 if( evt ){
00307 runNo = evt -> runNumber();
00308 log << MSG::INFO <<"The runNumber of current event is "<<runNo<<endreq;
00309 }
00310 else {
00311 log << MSG::ERROR <<"Can not get EventHeader!"<<endreq;
00312 }
00313 #else
00314 int runNo = run;
00315 #endif
00316
00317 using FieldDBUtil::ConnectionDB;
00318 std::vector<double> current(3,0.);
00319 ConnectionDB::eRet e = m_connect_run->getReadSC_MagnetInfo(current,std::abs(runNo));
00320
00321 std::vector<double> beamEnergy;
00322 beamEnergy.clear();
00323 e = m_connect_run->getBeamEnergy(beamEnergy, std::abs(runNo));
00324 char BPR_PRB[200];
00325 char BER_PRB[200];
00326
00327 sprintf(BPR_PRB,"%f ",beamEnergy[0]);
00328 sprintf(BER_PRB,"%f ",beamEnergy[1]);
00329
00330 setenv("BEPCII_INFO.BPR_PRB", BPR_PRB, 1);
00331 setenv("BEPCII_INFO.BER_PRB", BER_PRB, 1);
00332
00333 int ssm_curr = int (current[0]);
00334 double scql_curr = current[1];
00335 double scqr_curr = current[2];
00336
00337 #ifndef BEAN
00338 log << MSG::INFO << "SSM current: " << current[0] << " SCQL current: " << current[1] << " SCQR current: " << current[2] << " in Run " << runNo << endreq;
00339 #else
00340 cout << "SSM current: " << current[0] << " SCQL current: " << current[1]
00341 << " SCQR current: " << current[2] << " in Run " << runNo << endl;
00342 #endif
00343
00344
00345
00346
00347
00348 if((ssm_curr >= 3367) && (ssm_curr <= 3370) && ((scql_curr+scqr_curr)/2 < m_Cur_SCQ1_89)) {
00349 m_zfield = -1.0*tesla;
00350 m_filename = path;
00351 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode2.dat");
00352 #ifndef BEAN
00353 log << MSG::INFO << "Select field map: " << m_filename << endreq;
00354 #else
00355 cout << "Select field map: " << m_filename << endl;
00356 #endif
00357
00358 StatusCode status = parseFile();
00359
00360 #ifndef BEAN
00361 if ( !status.isSuccess() ) {
00362 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
00363 }
00364 #else
00365 if ( !status ) {
00366 cout << "Magnetic field parse failled" << endl;
00367 }
00368 #endif
00369 m_Q_1 = m_Q;
00370 m_P_1 = m_P;
00371
00372 m_Q.clear();
00373 m_P.clear();
00374 if(m_gridDistance == 5) {
00375 m_Q.reserve(201250);
00376 m_P.reserve(201250);
00377 }
00378 if(m_gridDistance == 10) {
00379 m_Q.reserve(27082);
00380 m_P.reserve(27082);
00381 }
00382
00383 m_filename = path;
00384 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode3.dat");
00385 #ifndef BEAN
00386 log << MSG::INFO << "Select field map: " << m_filename << endreq;
00387 #else
00388 cout << "Select field map: " << m_filename << endl;
00389 #endif
00390
00391 status = parseFile();
00392
00393 #ifndef BEAN
00394 if ( !status.isSuccess() ) {
00395 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
00396 }
00397 #else
00398 if ( !status ) {
00399 cout << "Magnetic field parse failled" << endl;
00400 }
00401 #endif
00402 m_Q_2 = m_Q;
00403 m_P_2 = m_P;
00404
00405 m_Q.clear();
00406 m_P.clear();
00407 if(m_gridDistance == 5) {
00408 m_Q.reserve(201250);
00409 m_P.reserve(201250);
00410 }
00411 if(m_gridDistance == 10) {
00412 m_Q.reserve(27082);
00413 m_P.reserve(27082);
00414 }
00415
00416 if(m_Q_2.size() != m_Q_1.size()) {
00417 #ifndef BEAN
00418 log << MSG::FATAL << "The two field maps used in this run are wrong!!!" << endreq;
00419 #else
00420 cout << "The two field maps used in this run are wrong!!!" << endl;
00421 #endif
00422 exit(1);
00423 }
00424
00425 for(unsigned int iQ = 0; iQ < m_Q_1.size(); iQ++) {
00426 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];
00427 m_Q.push_back(fieldvalue);
00428 m_P.push_back(m_P_2[iQ]);
00429 }
00430 }
00431
00432 if((ssm_curr >= 3367) && (ssm_curr <= 3370) && ((scql_curr+scqr_curr)/2 >= m_Cur_SCQ1_89)) {
00433 m_zfield = -1.0*tesla;
00434 m_filename = path;
00435 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode3.dat");
00436 #ifndef BEAN
00437 log << MSG::INFO << "Select field map: " << m_filename << endreq;
00438 #else
00439 cout << "Select field map: " << m_filename << endl;
00440 #endif
00441
00442 StatusCode status = parseFile();
00443
00444 #ifndef BEAN
00445 if ( !status.isSuccess() ) {
00446 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
00447 }
00448 #else
00449 if ( !status ) {
00450 cout << "Magnetic field parse failled" << endl;
00451 }
00452 #endif
00453 m_Q_1 = m_Q;
00454 m_P_1 = m_P;
00455
00456 m_Q.clear();
00457 m_P.clear();
00458 if(m_gridDistance == 5) {
00459 m_Q.reserve(201250);
00460 m_P.reserve(201250);
00461 }
00462 if(m_gridDistance == 10) {
00463 m_Q.reserve(27082);
00464 m_P.reserve(27082);
00465 }
00466
00467 m_filename = path;
00468 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode4.dat");
00469 #ifndef BEAN
00470 log << MSG::INFO << "Select field map: " << m_filename << endreq;
00471 #else
00472 cout << "Select field map: " << m_filename << endl;
00473 #endif
00474
00475 status = parseFile();
00476
00477 #ifndef BEAN
00478 if ( !status.isSuccess() ) {
00479 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
00480 }
00481 #else
00482 if ( !status ) {
00483 cout << "Magnetic field parse failled" << endl;
00484 }
00485 #endif
00486 m_Q_2 = m_Q;
00487 m_P_2 = m_P;
00488
00489 m_Q.clear();
00490 m_P.clear();
00491 if(m_gridDistance == 5) {
00492 m_Q.reserve(201250);
00493 m_P.reserve(201250);
00494 }
00495 if(m_gridDistance == 10) {
00496 m_Q.reserve(27082);
00497 m_P.reserve(27082);
00498 }
00499
00500 if(m_Q_2.size() != m_Q_1.size()) {
00501 #ifndef BEAN
00502 log << MSG::FATAL << "The two field maps used in this run are wrong!!!" << endreq;
00503 #else
00504 cout << "The two field maps used in this run are wrong!!!" << endl;
00505 #endif
00506 exit(1);
00507 }
00508
00509 for(unsigned int iQ = 0; iQ < m_Q_1.size(); iQ++) {
00510 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];
00511 m_Q.push_back(fieldvalue);
00512 m_P.push_back(m_P_2[iQ]);
00513 }
00514 }
00515 if((ssm_curr >= 3033) && (ssm_curr <= 3035)) {
00516 m_zfield = -0.9*tesla;
00517 m_filename = path;
00518 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode7.dat");
00519 #ifndef BEAN
00520 log << MSG::INFO << "Select field map: " << m_filename << endreq;
00521 #else
00522 cout << "Select field map: " << m_filename << endl;
00523 #endif
00524
00525 StatusCode status = parseFile();
00526
00527 #ifndef BEAN
00528 if ( !status.isSuccess() ) {
00529 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
00530 }
00531 #else
00532 if ( !status ) {
00533 cout << "Magnetic field parse failled" << endl;
00534 }
00535 #endif
00536 m_Q_1 = m_Q;
00537 m_P_1 = m_P;
00538
00539 m_Q.clear();
00540 m_P.clear();
00541 if(m_gridDistance == 5) {
00542 m_Q.reserve(201250);
00543 m_P.reserve(201250);
00544 }
00545 if(m_gridDistance == 10) {
00546 m_Q.reserve(27082);
00547 m_P.reserve(27082);
00548 }
00549
00550 m_filename = path;
00551 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode8.dat");
00552 #ifndef BEAN
00553 log << MSG::INFO << "Select field map: " << m_filename << endreq;
00554 #else
00555 cout << "Select field map: " << m_filename << endl;
00556 #endif
00557
00558 status = parseFile();
00559
00560 #ifndef BEAN
00561 if ( !status.isSuccess() ) {
00562 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
00563 }
00564 #else
00565 if ( !status ) {
00566 cout << "Magnetic field parse failled" << endl;
00567 }
00568 #endif
00569 m_Q_2 = m_Q;
00570 m_P_2 = m_P;
00571
00572 m_Q.clear();
00573 m_P.clear();
00574 if(m_gridDistance == 5) {
00575 m_Q.reserve(201250);
00576 m_P.reserve(201250);
00577 }
00578 if(m_gridDistance == 10) {
00579 m_Q.reserve(27082);
00580 m_P.reserve(27082);
00581 }
00582
00583 if(m_Q_2.size() != m_Q_1.size()) {
00584 #ifndef BEAN
00585 log << MSG::FATAL << "The two field maps used in this run are wrong!!!" << endreq;
00586 #else
00587 cout << "The two field maps used in this run are wrong!!!" << endl;
00588 #endif
00589 exit(1);
00590 }
00591
00592 for(unsigned int iQ = 0; iQ < m_Q_1.size(); iQ++) {
00593 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];
00594 m_Q.push_back(fieldvalue);
00595 m_P.push_back(m_P_2[iQ]);
00596 }
00597 }
00598
00599 if((!((ssm_curr >= 3367) && (ssm_curr <= 3370)) && !((ssm_curr >= 3033) && (ssm_curr <= 3035))) || scql_curr == 0 || scqr_curr == 0) {
00600 #ifndef BEAN
00601 log << MSG::FATAL << "The current of run " << runNo << " is abnormal in database, please check it! " << endreq;
00602 log << MSG::FATAL << "The current of SSM is " << ssm_curr << ", SCQR is " << scqr_curr << ", SCQL is " << scql_curr << endreq;
00603 #else
00604 cout << "The current of run " << runNo
00605 << " is abnormal in database, please check it! " << endl;
00606 cout << "The current of SSM is " << ssm_curr
00607 << ", SCQR is " << scqr_curr << ", SCQL is " << scql_curr << endl;
00608 #endif
00609 exit(1);
00610 }
00611
00612 if(m_Q.size() == 0) {
00613 #ifndef BEAN
00614 log << MSG::FATAL << "This size of field map vector is ZERO, please check the current of run " << runNo << " in database!" << endreq;
00615 #else
00616 cout << "This size of field map vector is ZERO,"
00617 << " please check the current of run " << runNo
00618 << " in database!" << endl;
00619 #endif
00620 exit(1);
00621 }
00622
00623 m_filename_TE = path;
00624 if(m_gridDistance == 10) m_filename_TE += std::string( "/dat/TEMap10cm.dat");
00625 if(m_gridDistance == 5) m_filename_TE += std::string( "/dat/TEMap5cm.dat");
00626 #ifndef BEAN
00627 log << MSG::INFO << "Select field map: " << m_filename_TE << endreq;
00628 #else
00629 cout << "Select field map: " << m_filename_TE << endl;
00630 #endif
00631
00632 StatusCode status = parseFile_TE();
00633
00634 #ifndef BEAN
00635 if ( !status.isSuccess() ) {
00636 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
00637 }
00638 #else
00639 if ( !status ) {
00640 cout << "Magnetic field parse failled" << endl;
00641 }
00642 #endif
00643
00644 for (int iC = 0; iC< 2; ++iC ){
00645 m_min_FL[iC] = -90.0*cm;
00646 m_max_FL[iC] = m_min_FL[iC]+( m_Nxyz[iC]-1 )* m_Dxyz[iC];
00647
00648 m_min_FL_TE[iC] = 0.0*cm;
00649 m_max_FL_TE[iC] = m_min_FL_TE[iC]+( m_Nxyz_TE[iC]-1 )* m_Dxyz_TE[iC];
00650 }
00651 m_min_FL[2] = -120.0*cm;
00652 m_max_FL[2] = m_min_FL[2]+( m_Nxyz[2]-1 )* m_Dxyz[2];
00653
00654 m_min_FL_TE[2] = 0.0*cm;
00655 m_max_FL_TE[2] = m_min_FL_TE[2]+( m_Nxyz_TE[2]-1 )* m_Dxyz_TE[2];
00656 }
00657
00658 #ifndef BEAN
00659
00660
00661
00662 StatusCode MagneticFieldSvc::finalize()
00663 {
00664 MsgStream log( msgSvc(), name() );
00665
00666 StatusCode status = Service::finalize();
00667
00668 if ( status.isSuccess() )
00669 log << MSG::INFO << "Service finalized successfully" << endreq;
00670 return status;
00671 }
00672
00673
00674
00675
00676 StatusCode MagneticFieldSvc::queryInterface( const InterfaceID& riid,
00677 void** ppvInterface )
00678 {
00679 StatusCode sc = StatusCode::FAILURE;
00680 if ( ppvInterface ) {
00681 *ppvInterface = 0;
00682
00683 if ( riid == IID_IMagneticFieldSvc ) {
00684 *ppvInterface = static_cast<IMagneticFieldSvc*>(this);
00685 sc = StatusCode::SUCCESS;
00686 addRef();
00687 }
00688 else
00689 sc = Service::queryInterface( riid, ppvInterface );
00690 }
00691 return sc;
00692 }
00693
00694 void MagneticFieldSvc::handle(const Incident& inc) {
00695 MsgStream log( messageService(), name() );
00696 log << MSG::DEBUG << "handle: " << inc.type() << endreq;
00697 if ( inc.type() != "NewRun" ){
00698 return;
00699 }
00700 log << MSG::DEBUG << "Begin New Run" << endreq;
00701 if(m_useDB == true) {
00702 init_params();
00703 }
00704 }
00705
00706 #else
00707 void MagneticFieldSvc::handle(int new_run) {
00708 static int save_run = 0;
00709 if( new_run == save_run ) return;
00710
00711 cout << "Begin New Run " << new_run << endl;
00712 save_run = new_run;
00713 if(m_useDB == true) {
00714 init_params(new_run);
00715 }
00716 }
00717 #endif
00718
00719
00720
00721
00722
00723 StatusCode MagneticFieldSvc::parseFile() {
00724 #ifndef BEAN
00725 StatusCode sc = StatusCode::FAILURE;
00726
00727 MsgStream log( msgSvc(), name() );
00728 #else
00729 StatusCode sc = false;
00730 #endif
00731
00732 char line[ 255 ];
00733 std::ifstream infile( m_filename.c_str() );
00734
00735 if ( infile ) {
00736 #ifndef BEAN
00737 sc = StatusCode::SUCCESS;
00738 #else
00739 sc = true;
00740 #endif
00741
00742
00743 do{
00744 infile.getline( line, 255 );
00745 } while( line[0] != 'P' );
00746
00747
00748 std::string sPar[2];
00749 char* token = strtok( line, " " );
00750 int ip = 0;
00751 do{
00752 if ( token ) { sPar[ip] = token; token = strtok( NULL, " " );}
00753 else continue;
00754 ip++;
00755 } while ( token != NULL );
00756 long int npar = atoi( sPar[1].c_str() );
00757
00758
00759 do{
00760 infile.getline( line, 255 );
00761 } while( line[0] != 'G' );
00762
00763
00764 do{
00765 infile.getline( line, 255 );
00766 } while( line[0] != '#' );
00767
00768
00769 infile.getline( line, 255 );
00770 std::string sGeom[7];
00771 token = strtok( line, " " );
00772 int ig = 0;
00773 do{
00774 if ( token ) { sGeom[ig] = token; token = strtok( NULL, " " );}
00775 else continue;
00776 ig++;
00777 } while (token != NULL);
00778
00779
00780 m_Dxyz[0] = atof( sGeom[0].c_str() ) * cm;
00781 m_Dxyz[1] = atof( sGeom[1].c_str() ) * cm;
00782 m_Dxyz[2] = atof( sGeom[2].c_str() ) * cm;
00783 m_Nxyz[0] = atoi( sGeom[3].c_str() );
00784 m_Nxyz[1] = atoi( sGeom[4].c_str() );
00785 m_Nxyz[2] = atoi( sGeom[5].c_str() );
00786 m_zOffSet = atof( sGeom[6].c_str() ) * cm;
00787
00788 long int nlines = ( npar - 7 ) / 3;
00789
00790
00791 long int ncheck = 0;
00792
00793
00794
00795
00796 while( infile ) {
00797
00798
00799 infile.getline( line, 255 );
00800 if ( line[0] == '#' ) continue;
00801 std::string gridx, gridy, gridz, sFx, sFy, sFz;
00802 char* token = strtok( line, " " );
00803 if ( token ) { gridx = token; token = strtok( NULL, " " );} else continue;
00804 if ( token ) { gridy = token; token = strtok( NULL, " " );} else continue;
00805 if ( token ) { gridz = token; token = strtok( NULL, " " );} else continue;
00806 if ( token ) { sFx = token; token = strtok( NULL, " " );} else continue;
00807 if ( token ) { sFy = token; token = strtok( NULL, " " );} else continue;
00808 if ( token ) { sFz = token; token = strtok( NULL, " " );} else continue;
00809 if ( token != NULL ) continue;
00810
00811 double gx = atof( gridx.c_str() ) * m;
00812 double gy = atof( gridy.c_str() ) * m;
00813 double gz = atof( gridz.c_str() ) * m;
00814
00815 double fx = atof( sFx.c_str() ) * tesla;
00816 double fy = atof( sFy.c_str() ) * tesla;
00817 double fz = atof( sFz.c_str() ) * tesla;
00818
00819 if(m_outlevel == 0) {
00820 #ifndef BEAN
00821 log << MSG::DEBUG << "grid x, y, z is "<< gx << ", " << gy << ", " << gz << endreq;
00822 log << MSG::DEBUG << "field x, y, z is "<< fx << ", " << fy << ", " << fz << endreq;
00823 #else
00824 cout << "grid x, y, z is "<< gx << ", " << gy << ", " << gz << endl;
00825 cout << "field x, y, z is "<< fx << ", " << fy << ", " << fz << endl;
00826 #endif
00827 }
00828
00829 m_P.push_back( gx );
00830 m_P.push_back( gy );
00831 m_P.push_back( gz );
00832
00833
00834 m_Q.push_back( fx );
00835 m_Q.push_back( fy );
00836 m_Q.push_back( fz );
00837
00838 ncheck++;
00839 }
00840 infile.close();
00841 if ( nlines != ncheck) {
00842 #ifndef BEAN
00843 log << MSG::ERROR << "nline is " << nlines << "; ncheck is " << ncheck << " Number of points in field map does not match!"
00844 << endreq;
00845 return StatusCode::FAILURE;
00846 #else
00847 cout << "nline is " << nlines << "; ncheck is " << ncheck << " Number of points in field map does not match!"
00848 << endl;
00849 return false;
00850 #endif
00851 }
00852 }
00853 else {
00854 #ifndef BEAN
00855 log << MSG::ERROR << "Unable to open magnetic field file : "
00856 << m_filename << endreq;
00857 #else
00858 cout << "Unable to open magnetic field file : "
00859 << m_filename << endl;
00860 #endif
00861 }
00862
00863 return sc;
00864 }
00865
00866
00867
00868
00869
00870
00871 StatusCode MagneticFieldSvc::parseFile_TE() {
00872 #ifndef BEAN
00873 StatusCode sc = StatusCode::FAILURE;
00874
00875 MsgStream log( msgSvc(), name() );
00876 #else
00877 StatusCode sc = false;
00878 #endif
00879
00880 char line[ 255 ];
00881 std::ifstream infile( m_filename_TE.c_str() );
00882
00883 if ( infile ) {
00884 #ifndef BEAN
00885 sc = StatusCode::SUCCESS;
00886 #else
00887 sc = true;
00888 #endif
00889
00890
00891 do{
00892 infile.getline( line, 255 );
00893 } while( line[0] != 'P' );
00894
00895
00896 std::string sPar[2];
00897 char* token = strtok( line, " " );
00898 int ip = 0;
00899 do{
00900 if ( token ) { sPar[ip] = token; token = strtok( NULL, " " );}
00901 else continue;
00902 ip++;
00903 } while ( token != NULL );
00904 long int npar = atoi( sPar[1].c_str() );
00905
00906
00907 do{
00908 infile.getline( line, 255 );
00909 } while( line[0] != 'G' );
00910
00911
00912 do{
00913 infile.getline( line, 255 );
00914 } while( line[0] != '#' );
00915
00916
00917 infile.getline( line, 255 );
00918 std::string sGeom[7];
00919 token = strtok( line, " " );
00920 int ig = 0;
00921 do{
00922 if ( token ) { sGeom[ig] = token; token = strtok( NULL, " " );}
00923 else continue;
00924 ig++;
00925 } while (token != NULL);
00926
00927
00928 m_Dxyz_TE[0] = atof( sGeom[0].c_str() ) * cm;
00929 m_Dxyz_TE[1] = atof( sGeom[1].c_str() ) * cm;
00930 m_Dxyz_TE[2] = atof( sGeom[2].c_str() ) * cm;
00931 m_Nxyz_TE[0] = atoi( sGeom[3].c_str() );
00932 m_Nxyz_TE[1] = atoi( sGeom[4].c_str() );
00933 m_Nxyz_TE[2] = atoi( sGeom[5].c_str() );
00934 m_zOffSet_TE = atof( sGeom[6].c_str() ) * cm;
00935
00936 long int nlines = ( npar - 7 ) / 3;
00937
00938
00939 long int ncheck = 0;
00940
00941
00942
00943
00944 while( infile ) {
00945
00946
00947 infile.getline( line, 255 );
00948 if ( line[0] == '#' ) continue;
00949 std::string gridx, gridy, gridz, sFx, sFy, sFz;
00950 char* token = strtok( line, " " );
00951 if ( token ) { gridx = token; token = strtok( NULL, " " );} else continue;
00952 if ( token ) { gridy = token; token = strtok( NULL, " " );} else continue;
00953 if ( token ) { gridz = token; token = strtok( NULL, " " );} else continue;
00954 if ( token ) { sFx = token; token = strtok( NULL, " " );} else continue;
00955 if ( token ) { sFy = token; token = strtok( NULL, " " );} else continue;
00956 if ( token ) { sFz = token; token = strtok( NULL, " " );} else continue;
00957 if ( token != NULL ) continue;
00958
00959 double gx = atof( gridx.c_str() ) * m;
00960 double gy = atof( gridy.c_str() ) * m;
00961 double gz = atof( gridz.c_str() ) * m;
00962
00963 double fx = atof( sFx.c_str() ) * tesla;
00964 double fy = atof( sFy.c_str() ) * tesla;
00965 double fz = atof( sFz.c_str() ) * tesla;
00966
00967 if(m_outlevel == 0) {
00968 #ifndef BEAN
00969 log << MSG::DEBUG << "grid x, y, z is "<< gx << ", " << gy << ", " << gz << endreq;
00970 log << MSG::DEBUG << "field x, y, z is "<< fx << ", " << fy << ", " << fz << endreq;
00971 #else
00972 cout << "grid x, y, z is "<< gx << ", " << gy << ", " << gz << endl;
00973 cout << "field x, y, z is "<< fx << ", " << fy << ", " << fz << endl;
00974 #endif
00975 }
00976
00977 m_P_TE.push_back( gx );
00978 m_P_TE.push_back( gy );
00979 m_P_TE.push_back( gz );
00980
00981
00982 m_Q_TE.push_back( fx );
00983 m_Q_TE.push_back( fy );
00984 m_Q_TE.push_back( fz );
00985
00986 ncheck++;
00987 }
00988 infile.close();
00989 if ( nlines != ncheck) {
00990 #ifndef BEAN
00991 log << MSG::ERROR << "nline is " << nlines << "; ncheck is " << ncheck << " Number of points in field map does not match!"
00992 << endreq;
00993 return StatusCode::FAILURE;
00994 #else
00995 cout << "nline is " << nlines << "; ncheck is " << ncheck << " Number of points in field map does not match!"
00996 << endl;
00997 return false;
00998 #endif
00999 }
01000 }
01001 else {
01002 #ifndef BEAN
01003 log << MSG::ERROR << "Unable to open magnetic field file : "
01004 << m_filename_TE << endreq;
01005 #else
01006 cout << "Unable to open magnetic field file : "
01007 << m_filename_TE << endl;
01008 #endif
01009 }
01010
01011 return sc;
01012 }
01013
01014
01015
01016 StatusCode MagneticFieldSvc::fieldVector(const HepPoint3D& newr,
01017 HepVector3D& newb) const
01018 {
01019
01020 if(m_turnOffField == true) {
01021 newb[0] = 0.;
01022 newb[1] = 0.;
01023 newb[2] = 0.;
01024 #ifndef BEAN
01025 return StatusCode::SUCCESS;
01026 #else
01027 return true;
01028 #endif
01029 }
01030
01031
01032 if(m_uniField) {
01033 uniFieldVector(newr,newb);
01034 #ifndef BEAN
01035 return StatusCode::SUCCESS;
01036 #else
01037 return true;
01038 #endif
01039 }
01040
01041
01042
01043 double new_x = -newr.x();
01044 double new_y = newr.y();
01045 double new_z = -newr.z();
01046 HepPoint3D r(new_x,new_y,new_z);
01047 HepVector3D b;
01048 b[0]=0.0*tesla;
01049 b[1]=0.0*tesla;
01050 b[2]=0.0*tesla;
01051
01052 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)
01053 {
01054 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)
01055
01056 {
01057 this->fieldGrid( r, b );
01058 }
01059 else
01060 {
01061 this->fieldGrid_TE( r, b );
01062 }
01063 }
01064 if((fabs(r.z())<=1970*mm && sqrt(r.x()*r.x()+r.y()*r.y()) >= 1740*mm) || (fabs(r.z())>=2050*mm))
01065 {
01066 HepPoint3D mr;
01067 HepVector3D tb;
01068 int part = 0, layer = 0, mat = 0;
01069 double theta;
01070 bool ifbar = false;
01071 bool ifend = false;
01072 if(r.x()!=0.){
01073 theta = atan2(fabs(r.y()),fabs(r.x()));
01074 if(0<=theta&&theta<pi/8) {
01075 mr[0] = fabs(r.x()); mr[1] = -fabs(r.y()); mr[2] = fabs(r.z());
01076 if(mr[2]<=1970*mm&&1740*mm<=mr[0]&&mr[0]<=2620*mm){
01077 part = 1;
01078 if(1740*mm<=mr[0]&&mr[0]<1770*mm) { layer = 0; mat = 0; }
01079 if(1770*mm<=mr[0]&&mr[0]<1810*mm) { layer = 0; mat = 1; }
01080 if(1810*mm<=mr[0]&&mr[0]<1840*mm) { layer = 1; mat = 0; }
01081 if(1840*mm<=mr[0]&&mr[0]<1880*mm) { layer = 1; mat = 1; }
01082 if(1880*mm<=mr[0]&&mr[0]<1910*mm) { layer = 2; mat = 0; }
01083 if(1910*mm<=mr[0]&&mr[0]<1950*mm) { layer = 2; mat = 1; }
01084 if(1950*mm<=mr[0]&&mr[0]<1990*mm) { layer = 3; mat = 0; }
01085 if(1990*mm<=mr[0]&&mr[0]<2030*mm) { layer = 3; mat = 1; }
01086 if(2030*mm<=mr[0]&&mr[0]<2070*mm) { layer = 4; mat = 0; }
01087 if(2070*mm<=mr[0]&&mr[0]<2110*mm) { layer = 4; mat = 1; }
01088 if(2110*mm<=mr[0]&&mr[0]<2190*mm) { layer = 5; mat = 0; }
01089 if(2190*mm<=mr[0]&&mr[0]<2230*mm) { layer = 5; mat = 1; }
01090 if(2230*mm<=mr[0]&&mr[0]<2310*mm) { layer = 6; mat = 0; }
01091 if(2310*mm<=mr[0]&&mr[0]<2350*mm) { layer = 6; mat = 1; }
01092 if(2350*mm<=mr[0]&&mr[0]<2430*mm) { layer = 7; mat = 0; }
01093 if(2430*mm<=mr[0]&&mr[0]<2470*mm) { layer = 7; mat = 1; }
01094 if(2470*mm<=mr[0]&&mr[0]<=2620*mm) { layer = 8; mat = 0; }
01095 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01096 b[0] = tb[0];
01097 b[1] = -tb[1];
01098 b[2] = tb[2];
01099 ifbar = true;
01100 }
01101 if(2050*mm<=mr[2]&&mr[2]<2090*mm&&1034*mm<=mr[0]&&mr[0]<=2500*mm){
01102 part = 0; layer = 0; mat = 0;
01103 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01104 b[0] = tb[0];
01105 b[1] = -tb[1];
01106 b[2] = tb[2];
01107 ifend = true;
01108 }
01109 if(2090*mm<=mr[2]&&mr[2]<2130*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
01110 part = 0; layer = 0; mat = 1;
01111 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01112 b[0] = tb[0];
01113 b[1] = -tb[1];
01114 b[2] = tb[2];
01115 ifend = true;
01116 }
01117 if(2130*mm<=mr[2]&&mr[2]<2170*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
01118 part = 0; layer = 1; mat = 0;
01119 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01120 b[0] = tb[0];
01121 b[1] = -tb[1];
01122 b[2] = tb[2];
01123 ifend = true;
01124 }
01125 if(2170*mm<=mr[2]&&mr[2]<2210*mm&&1100*mm<=mr[0]&&mr[0]<=2500*mm) {
01126 part = 0; layer = 1; mat = 1;
01127 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01128 b[0] = tb[0];
01129 b[1] = -tb[1];
01130 b[2] = tb[2];
01131 ifend = true;
01132 }
01133 if(2210*mm<=mr[2]&&mr[2]<2240*mm&&1100*mm<mr[0]&&mr[0]<=2500*mm) {
01134 part = 0; layer = 2; mat = 0;
01135 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01136 b[0] = tb[0];
01137 b[1] = -tb[1];
01138 b[2] = tb[2];
01139 ifend = true;
01140 }
01141 if(2240*mm<=mr[2]&&mr[2]<2280*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
01142 part = 0; layer = 2; mat = 1;
01143 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01144 b[0] = tb[0];
01145 b[1] = -tb[1];
01146 b[2] = tb[2];
01147 ifend = true;
01148 }
01149 if(2280*mm<=mr[2]&&mr[2]<2310*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
01150 part = 0; layer = 3; mat = 0;
01151 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01152 b[0] = tb[0];
01153 b[1] = -tb[1];
01154 b[2] = tb[2];
01155 ifend = true;
01156 }
01157 if(2310*mm<=mr[2]&&mr[2]<2350*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
01158 part = 0; layer = 3; mat = 1;
01159 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01160 b[0] = tb[0];
01161 b[1] = -tb[1];
01162 b[2] = tb[2];
01163 ifend = true;
01164 }
01165 if(2350*mm<=mr[2]&&mr[2]<2380*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
01166 part = 0; layer = 4; mat = 0;
01167 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01168 b[0] = tb[0];
01169 b[1] = -tb[1];
01170 b[2] = tb[2];
01171 ifend = true;
01172 }
01173 if(2380*mm<=mr[2]&&mr[2]<2420*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
01174 part = 0; layer = 4; mat = 1;
01175 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01176 b[0] = tb[0];
01177 b[1] = -tb[1];
01178 b[2] = tb[2];
01179 ifend = true;
01180 }
01181 if(2420*mm<=mr[2]&&mr[2]<2470*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
01182 part = 0; layer = 5; mat = 0;
01183 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01184 b[0] = tb[0];
01185 b[1] = -tb[1];
01186 b[2] = tb[2];
01187 ifend = true;
01188 }
01189 if(2470*mm<=mr[2]&&mr[2]<2510*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
01190 part = 0; layer = 5; mat =1;
01191 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01192 b[0] = tb[0];
01193 b[1] = -tb[1];
01194 b[2] = tb[2];
01195 ifend = true;
01196 }
01197 if(2510*mm<=mr[2]&&mr[2]<2590*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
01198 part = 0; layer = 6; mat = 0;
01199 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01200 b[0] = tb[0];
01201 b[1] = -tb[1];
01202 b[2] = tb[2];
01203 ifend = true;
01204 }
01205 if(2590*mm<=mr[2]&&mr[2]<2630*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01206 part = 0; layer = 6; mat = 1;
01207 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01208 b[0] = tb[0];
01209 b[1] = -tb[1];
01210 b[2] = tb[2];
01211 ifend = true;
01212 }
01213 if(2630*mm<=mr[2]&&mr[2]<2710*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01214 part = 0; layer = 7; mat = 0;
01215 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01216 b[0] = tb[0];
01217 b[1] = -tb[1];
01218 b[2] = tb[2];
01219 ifend = true;
01220 }
01221 if(2710*mm<=mr[2]&&mr[2]<2750*mm&&1362*mm<=mr[0]&&mr[0]<=2500*mm) {
01222 part = 0; layer = 7; mat = 1;
01223 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01224 b[0] = tb[0];
01225 b[1] = -tb[1];
01226 b[2] = tb[2];
01227 ifend = true;
01228 }
01229 if(2750*mm<=mr[2]&&mr[2]<=2800*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01230 part = 0; layer = 8; mat = 0;
01231 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01232 b[0] = tb[0];
01233 b[1] = -tb[1];
01234 b[2] = tb[2];
01235 ifend = true;
01236 }
01237 }
01238 if(pi/8<=theta&&theta<pi/4) {
01239 mr[0] = fabs(r.x())*cos(pi/4)+fabs(r.y())*sin(pi/4);
01240 mr[1] = -fabs(r.x())*sin(pi/4)+fabs(r.y())*cos(pi/4);
01241 mr[2] = fabs(r.z());
01242 if(mr[2]<=1970*mm&&1740*mm<=mr[0]&&mr[0]<=2620*mm) {
01243 part = 1;
01244 if(1740*mm<=mr[0]&&mr[0]<1770*mm) { layer = 0; mat = 0; }
01245 if(1770*mm<=mr[0]&&mr[0]<1810*mm) { layer = 0; mat = 1; }
01246 if(1810*mm<=mr[0]&&mr[0]<1840*mm) { layer = 1; mat = 0; }
01247 if(1840*mm<=mr[0]&&mr[0]<1880*mm) { layer = 1; mat = 1; }
01248 if(1880*mm<=mr[0]&&mr[0]<1910*mm) { layer = 2; mat = 0; }
01249 if(1910*mm<=mr[0]&&mr[0]<1950*mm) { layer = 2; mat = 1; }
01250 if(1950*mm<=mr[0]&&mr[0]<1990*mm) { layer = 3; mat = 0; }
01251 if(1990*mm<=mr[0]&&mr[0]<2030*mm) { layer = 3; mat = 1; }
01252 if(2030*mm<=mr[0]&&mr[0]<2070*mm) { layer = 4; mat = 0; }
01253 if(2070*mm<=mr[0]&&mr[0]<2110*mm) { layer = 4; mat = 1; }
01254 if(2110*mm<=mr[0]&&mr[0]<2190*mm) { layer = 5; mat = 0; }
01255 if(2190*mm<=mr[0]&&mr[0]<2230*mm) { layer = 5; mat = 1; }
01256 if(2230*mm<=mr[0]&&mr[0]<2310*mm) { layer = 6; mat = 0; }
01257 if(2310*mm<=mr[0]&&mr[0]<2350*mm) { layer = 6; mat = 1; }
01258 if(2350*mm<=mr[0]&&mr[0]<2430*mm) { layer = 7; mat = 0; }
01259 if(2430*mm<=mr[0]&&mr[0]<2470*mm) { layer = 7; mat = 1; }
01260 if(2470*mm<=mr[0]&&mr[0]<=2620*mm) { layer = 8; mat = 0; }
01261 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01262 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01263 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01264 b[2] = tb[2];
01265 ifbar = true;
01266 }
01267 if(2050*mm<=mr[2]&&mr[2]<2090*mm&&1034*mm<=mr[0]&&mr[0]<=2500*mm){
01268 part = 0; layer = 0; mat = 0;
01269 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01270 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01271 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01272 b[2] = tb[2];
01273 ifend = true;
01274 }
01275 if(2090*mm<=mr[2]&&mr[2]<2130*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
01276 part = 0; layer = 0; mat = 1;
01277 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01278 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01279 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01280 b[2] = tb[2];
01281 ifend = true;
01282 }
01283 if(2130*mm<=mr[2]&&mr[2]<2170*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
01284 part = 0; layer = 1; mat = 0;
01285 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01286 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01287 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01288 b[2] = tb[2];
01289 ifend = true;
01290 }
01291 if(2170*mm<=mr[2]&&mr[2]<2210*mm&&1100*mm<=mr[0]&&mr[0]<=2500*mm) {
01292 part = 0; layer = 1; mat = 1;
01293 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01294 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01295 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01296 b[2] = tb[2];
01297 ifend = true;
01298 }
01299 if(2210*mm<=mr[2]&&mr[2]<2240*mm&&1100*mm<mr[0]&&mr[0]<=2500*mm) {
01300 part = 0; layer = 2; mat = 0;
01301 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01302 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01303 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01304 b[2] = tb[2];
01305 ifend = true;
01306 }
01307 if(2240*mm<=mr[2]&&mr[2]<2280*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
01308 part = 0; layer = 2; mat = 1;
01309 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01310 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01311 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01312 b[2] = tb[2];
01313 ifend = true;
01314 }
01315 if(2280*mm<=mr[2]&&mr[2]<2310*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
01316 part = 0; layer = 3; mat = 0;
01317 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01318 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01319 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01320 b[2] = tb[2];
01321 ifend = true;
01322 }
01323 if(2310*mm<=mr[2]&&mr[2]<2350*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
01324 part = 0; layer = 3; mat = 1;
01325 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01326 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01327 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01328 b[2] = tb[2];
01329 ifend = true;
01330 }
01331 if(2350*mm<=mr[2]&&mr[2]<2380*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
01332 part = 0; layer = 4; mat = 0;
01333 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01334 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01335 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01336 b[2] = tb[2];
01337 ifend = true;
01338 }
01339 if(2380*mm<=mr[2]&&mr[2]<2420*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
01340 part = 0; layer = 4; mat = 1;
01341 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01342 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01343 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01344 b[2] = tb[2];
01345 ifend = true;
01346 }
01347 if(2420*mm<=mr[2]&&mr[2]<2470*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
01348 part = 0; layer = 5; mat = 0;
01349 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01350 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01351 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01352 b[2] = tb[2];
01353 ifend = true;
01354 }
01355 if(2470*mm<=mr[2]&&mr[2]<2510*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
01356 part = 0; layer = 5; mat =1;
01357 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01358 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01359 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01360 b[2] = tb[2];
01361 ifend = true;
01362 }
01363 if(2510*mm<=mr[2]&&mr[2]<2590*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
01364 part = 0; layer = 6; mat = 0;
01365 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01366 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01367 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01368 b[2] = tb[2];
01369 ifend = true;
01370 }
01371 if(2590*mm<=mr[2]&&mr[2]<2630*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01372 part = 0; layer = 6; mat = 1;
01373 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01374 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01375 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01376 b[2] = tb[2];
01377 ifend = true;
01378 }
01379 if(2630*mm<=mr[2]&&mr[2]<2710*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01380 part = 0; layer = 7; mat = 0;
01381 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01382 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01383 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01384 b[2] = tb[2];
01385 ifend = true;
01386 }
01387 if(2710*mm<=mr[2]&&mr[2]<2750*mm&&1362*mm<=mr[0]&&mr[0]<=2500*mm) {
01388 part = 0; layer = 7; mat = 1;
01389 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01390 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01391 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01392 b[2] = tb[2];
01393 ifend = true;
01394 }
01395 if(2750*mm<=mr[2]&&mr[2]<=2800*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01396 part = 0; layer = 8; mat = 0;
01397 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01398 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
01399 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
01400 b[2] = tb[2];
01401 ifend = true;
01402 }
01403 }
01404 if(pi/4<=theta&&theta<3*pi/8) {
01405 mr[0] = fabs(r.x())*cos(pi/4)+fabs(r.y())*sin(pi/4);
01406 mr[1] = fabs(r.x())*sin(pi/4)-fabs(r.y())*cos(pi/4);
01407 mr[2] = fabs(r.z());
01408 if(mr[2]<=1970*mm&&1740*mm<=mr[0]&&mr[0]<=2620*mm) {
01409 part = 1;
01410 if(1740*mm<=mr[0]&&mr[0]<1770*mm) { layer = 0; mat = 0; }
01411 if(1770*mm<=mr[0]&&mr[0]<1810*mm) { layer = 0; mat = 1; }
01412 if(1810*mm<=mr[0]&&mr[0]<1840*mm) { layer = 1; mat = 0; }
01413 if(1840*mm<=mr[0]&&mr[0]<1880*mm) { layer = 1; mat = 1; }
01414 if(1880*mm<=mr[0]&&mr[0]<1910*mm) { layer = 2; mat = 0; }
01415 if(1910*mm<=mr[0]&&mr[0]<1950*mm) { layer = 2; mat = 1; }
01416 if(1950*mm<=mr[0]&&mr[0]<1990*mm) { layer = 3; mat = 0; }
01417 if(1990*mm<=mr[0]&&mr[0]<2030*mm) { layer = 3; mat = 1; }
01418 if(2030*mm<=mr[0]&&mr[0]<2070*mm) { layer = 4; mat = 0; }
01419 if(2070*mm<=mr[0]&&mr[0]<2110*mm) { layer = 4; mat = 1; }
01420 if(2110*mm<=mr[0]&&mr[0]<2190*mm) { layer = 5; mat = 0; }
01421 if(2190*mm<=mr[0]&&mr[0]<2230*mm) { layer = 5; mat = 1; }
01422 if(2230*mm<=mr[0]&&mr[0]<2310*mm) { layer = 6; mat = 0; }
01423 if(2310*mm<=mr[0]&&mr[0]<2350*mm) { layer = 6; mat = 1; }
01424 if(2350*mm<=mr[0]&&mr[0]<2430*mm) { layer = 7; mat = 0; }
01425 if(2430*mm<=mr[0]&&mr[0]<2470*mm) { layer = 7; mat = 1; }
01426 if(2470*mm<=mr[0]&&mr[0]<=2620*mm) { layer = 8; mat = 0; }
01427 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01428 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01429 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01430 b[2] = tb[2];
01431 ifbar = true;
01432 }
01433 if(2050*mm<=mr[2]&&mr[2]<2090*mm&&1034*mm<=mr[0]&&mr[0]<=2500*mm){
01434 part = 0; layer = 0; mat = 0;
01435 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01436 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01437 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01438 b[2] = tb[2];
01439 ifend = true;
01440 }
01441 if(2090*mm<=mr[2]&&mr[2]<2130*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
01442 part = 0; layer = 0; mat = 1;
01443 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01444 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01445 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01446 b[2] = tb[2];
01447 ifend = true;
01448 }
01449 if(2130*mm<=mr[2]&&mr[2]<2170*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
01450 part = 0; layer = 1; mat = 0;
01451 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01452 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01453 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01454 b[2] = tb[2];
01455 ifend = true;
01456 }
01457 if(2170*mm<=mr[2]&&mr[2]<2210*mm&&1100*mm<=mr[0]&&mr[0]<=2500*mm) {
01458 part = 0; layer = 1; mat = 1;
01459 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01460 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01461 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01462 b[2] = tb[2];
01463 ifend = true;
01464 }
01465 if(2210*mm<=mr[2]&&mr[2]<2240*mm&&1100*mm<mr[0]&&mr[0]<=2500*mm) {
01466 part = 0; layer = 2; mat = 0;
01467 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01468 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01469 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01470 b[2] = tb[2];
01471 ifend = true;
01472 }
01473 if(2240*mm<=mr[2]&&mr[2]<2280*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
01474 part = 0; layer = 2; mat = 1;
01475 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01476 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01477 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01478 b[2] = tb[2];
01479 ifend = true;
01480 }
01481 if(2280*mm<=mr[2]&&mr[2]<2310*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
01482 part = 0; layer = 3; mat = 0;
01483 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01484 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01485 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01486 b[2] = tb[2];
01487 ifend = true;
01488 }
01489 if(2310*mm<=mr[2]&&mr[2]<2350*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
01490 part = 0; layer = 3; mat = 1;
01491 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01492 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01493 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01494 b[2] = tb[2];
01495 ifend = true;
01496 }
01497 if(2350*mm<=mr[2]&&mr[2]<2380*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
01498 part = 0; layer = 4; mat = 0;
01499 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01500 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01501 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01502 b[2] = tb[2];
01503 ifend = true;
01504 }
01505 if(2380*mm<=mr[2]&&mr[2]<2420*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
01506 part = 0; layer = 4; mat = 1;
01507 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01508 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01509 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01510 b[2] = tb[2];
01511 ifend = true;
01512 }
01513 if(2420*mm<=mr[2]&&mr[2]<2470*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
01514 part = 0; layer = 5; mat = 0;
01515 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01516 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01517 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01518 b[2] = tb[2];
01519 ifend = true;
01520 }
01521 if(2470*mm<=mr[2]&&mr[2]<2510*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
01522 part = 0; layer = 5; mat =1;
01523 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01524 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01525 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01526 b[2] = tb[2];
01527 ifend = true;
01528 }
01529 if(2510*mm<=mr[2]&&mr[2]<2590*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
01530 part = 0; layer = 6; mat = 0;
01531 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01532 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01533 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01534 b[2] = tb[2];
01535 ifend = true;
01536 }
01537 if(2590*mm<=mr[2]&&mr[2]<2630*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01538 part = 0; layer = 6; mat = 1;
01539 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01540 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01541 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01542 b[2] = tb[2];
01543 ifend = true;
01544 }
01545 if(2630*mm<=mr[2]&&mr[2]<2710*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01546 part = 0; layer = 7; mat = 0;
01547 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01548 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01549 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01550 b[2] = tb[2];
01551 ifend = true;
01552 }
01553 if(2710*mm<=mr[2]&&mr[2]<2750*mm&&1362*mm<=mr[0]&&mr[0]<=2500*mm) {
01554 part = 0; layer = 7; mat = 1;
01555 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01556 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01557 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01558 b[2] = tb[2];
01559 ifend = true;
01560 }
01561 if(2750*mm<=mr[2]&&mr[2]<=2800*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01562 part = 0; layer = 8; mat = 0;
01563 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01564 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
01565 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
01566 b[2] = tb[2];
01567 ifend = true;
01568 }
01569 }
01570 if(3*pi/8<=theta&&theta<pi/2) {
01571 mr[0] = fabs(r.y()); mr[1] = -fabs(r.x()); mr[2] = fabs(r.z());
01572 if(mr[2]<=1970*mm&&1740*mm<=mr[0]&&mr[0]<=2620*mm) {
01573 part = 1;
01574 if(1740*mm<=mr[0]&&mr[0]<1770*mm) { layer = 0; mat = 0; }
01575 if(1770*mm<=mr[0]&&mr[0]<1810*mm) { layer = 0; mat = 1; }
01576 if(1810*mm<=mr[0]&&mr[0]<1840*mm) { layer = 1; mat = 0; }
01577 if(1840*mm<=mr[0]&&mr[0]<1880*mm) { layer = 1; mat = 1; }
01578 if(1880*mm<=mr[0]&&mr[0]<1910*mm) { layer = 2; mat = 0; }
01579 if(1910*mm<=mr[0]&&mr[0]<1950*mm) { layer = 2; mat = 1; }
01580 if(1950*mm<=mr[0]&&mr[0]<1990*mm) { layer = 3; mat = 0; }
01581 if(1990*mm<=mr[0]&&mr[0]<2030*mm) { layer = 3; mat = 1; }
01582 if(2030*mm<=mr[0]&&mr[0]<2070*mm) { layer = 4; mat = 0; }
01583 if(2070*mm<=mr[0]&&mr[0]<2110*mm) { layer = 4; mat = 1; }
01584 if(2110*mm<=mr[0]&&mr[0]<2190*mm) { layer = 5; mat = 0; }
01585 if(2190*mm<=mr[0]&&mr[0]<2230*mm) { layer = 5; mat = 1; }
01586 if(2230*mm<=mr[0]&&mr[0]<2310*mm) { layer = 6; mat = 0; }
01587 if(2310*mm<=mr[0]&&mr[0]<2350*mm) { layer = 6; mat = 1; }
01588 if(2350*mm<=mr[0]&&mr[0]<2430*mm) { layer = 7; mat = 0; }
01589 if(2430*mm<=mr[0]&&mr[0]<2470*mm) { layer = 7; mat = 1; }
01590 if(2470*mm<=mr[0]&&mr[0]<=2620*mm) { layer = 8; mat = 0; }
01591 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01592 b[0] = -tb[1];
01593 b[1] = tb[0];
01594 b[2] = tb[2];
01595 ifbar = true;
01596 }
01597 if(2050*mm<=mr[2]&&mr[2]<2090*mm&&1034*mm<=mr[0]&&mr[0]<=2500*mm){
01598 part = 0; layer = 0; mat = 0;
01599 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01600 b[0] = -tb[1];
01601 b[1] = tb[0];
01602 b[2] = tb[2];
01603 ifend = true;
01604 }
01605 if(2090*mm<=mr[2]&&mr[2]<2130*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
01606 part = 0; layer = 0; mat = 1;
01607 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01608 b[0] = -tb[1];
01609 b[1] = tb[0];
01610 b[2] = tb[2];
01611 ifend = true;
01612 }
01613 if(2130*mm<=mr[2]&&mr[2]<2170*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
01614 part = 0; layer = 1; mat = 0;
01615 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01616 b[0] = -tb[1];
01617 b[1] = tb[0];
01618 b[2] = tb[2];
01619 ifend = true;
01620 }
01621 if(2170*mm<=mr[2]&&mr[2]<2210*mm&&1100*mm<=mr[0]&&mr[0]<=2500*mm) {
01622 part = 0; layer = 1; mat = 1;
01623 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01624 b[0] = -tb[1];
01625 b[1] = tb[0];
01626 b[2] = tb[2];
01627 ifend = true;
01628 }
01629 if(2210*mm<=mr[2]&&mr[2]<2240*mm&&1100*mm<mr[0]&&mr[0]<=2500*mm) {
01630 part = 0; layer = 2; mat = 0;
01631 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01632 b[0] = -tb[1];
01633 b[1] = tb[0];
01634 b[2] = tb[2];
01635 ifend = true;
01636 }
01637 if(2240*mm<=mr[2]&&mr[2]<2280*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
01638 part = 0; layer = 2; mat = 1;
01639 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01640 b[0] = -tb[1];
01641 b[1] = tb[0];
01642 b[2] = tb[2];
01643 ifend = true;
01644 }
01645 if(2280*mm<=mr[2]&&mr[2]<2310*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
01646 part = 0; layer = 3; mat = 0;
01647 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01648 b[0] = -tb[1];
01649 b[1] = tb[0];
01650 b[2] = tb[2];
01651 ifend = true;
01652 }
01653 if(2310*mm<=mr[2]&&mr[2]<2350*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
01654 part = 0; layer = 3; mat = 1;
01655 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01656 b[0] = -tb[1];
01657 b[1] = tb[0];
01658 b[2] = tb[2];
01659 ifend = true;
01660 }
01661 if(2350*mm<=mr[2]&&mr[2]<2380*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
01662 part = 0; layer = 4; mat = 0;
01663 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01664 b[0] = -tb[1];
01665 b[1] = tb[0];
01666 b[2] = tb[2];
01667 ifend = true;
01668 }
01669 if(2380*mm<=mr[2]&&mr[2]<2420*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
01670 part = 0; layer = 4; mat = 1;
01671 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01672 b[0] = -tb[1];
01673 b[1] = tb[0];
01674 b[2] = tb[2];
01675 ifend = true;
01676 }
01677 if(2420*mm<=mr[2]&&mr[2]<2470*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
01678 part = 0; layer = 5; mat = 0;
01679 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01680 b[0] = -tb[1];
01681 b[1] = tb[0];
01682 b[2] = tb[2];
01683 ifend = true;
01684 }
01685 if(2470*mm<=mr[2]&&mr[2]<2510*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
01686 part = 0; layer = 5; mat =1;
01687 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01688 b[0] = -tb[1];
01689 b[1] = tb[0];
01690 b[2] = tb[2];
01691 ifend = true;
01692 }
01693 if(2510*mm<=mr[2]&&mr[2]<2590*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
01694 part = 0; layer = 6; mat = 0;
01695 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01696 b[0] = -tb[1];
01697 b[1] = tb[0];
01698 b[2] = tb[2];
01699 ifend = true;
01700 }
01701 if(2590*mm<=mr[2]&&mr[2]<2630*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01702 part = 0; layer = 6; mat = 1;
01703 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01704 b[0] = -tb[1];
01705 b[1] = tb[0];
01706 b[2] = tb[2];
01707 ifend = true;
01708 }
01709 if(2630*mm<=mr[2]&&mr[2]<2710*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01710 part = 0; layer = 7; mat = 0;
01711 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01712 b[0] = -tb[1];
01713 b[1] = tb[0];
01714 b[2] = tb[2];
01715 ifend = true;
01716 }
01717 if(2710*mm<=mr[2]&&mr[2]<2750*mm&&1362*mm<=mr[0]&&mr[0]<=2500*mm) {
01718 part = 0; layer = 7; mat = 1;
01719 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01720 b[0] = -tb[1];
01721 b[1] = tb[0];
01722 b[2] = tb[2];
01723 ifend = true;
01724 }
01725 if(2750*mm<=mr[2]&&mr[2]<=2800*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01726 part = 0; layer = 8; mat = 0;
01727 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01728 b[0] = -tb[1];
01729 b[1] = tb[0];
01730 b[2] = tb[2];
01731 ifend = true;
01732 }
01733 }
01734 }
01735 if(r.x()==0.) {
01736 mr[0] = fabs(r.y()); mr[1] = -fabs(r.x()); mr[2] = fabs(r.z());
01737 if(mr[2]<=1970*mm&&1740*mm<=mr[0]&&mr[0]<=2620*mm) {
01738 part = 1;
01739 if(1740*mm<=mr[0]&&mr[0]<1770*mm) { layer = 0; mat = 0; }
01740 if(1770*mm<=mr[0]&&mr[0]<1810*mm) { layer = 0; mat = 1; }
01741 if(1810*mm<=mr[0]&&mr[0]<1840*mm) { layer = 1; mat = 0; }
01742 if(1840*mm<=mr[0]&&mr[0]<1880*mm) { layer = 1; mat = 1; }
01743 if(1880*mm<=mr[0]&&mr[0]<1910*mm) { layer = 2; mat = 0; }
01744 if(1910*mm<=mr[0]&&mr[0]<1950*mm) { layer = 2; mat = 1; }
01745 if(1950*mm<=mr[0]&&mr[0]<1990*mm) { layer = 3; mat = 0; }
01746 if(1990*mm<=mr[0]&&mr[0]<2030*mm) { layer = 3; mat = 1; }
01747 if(2030*mm<=mr[0]&&mr[0]<2070*mm) { layer = 4; mat = 0; }
01748 if(2070*mm<=mr[0]&&mr[0]<2110*mm) { layer = 4; mat = 1; }
01749 if(2110*mm<=mr[0]&&mr[0]<2190*mm) { layer = 5; mat = 0; }
01750 if(2190*mm<=mr[0]&&mr[0]<2230*mm) { layer = 5; mat = 1; }
01751 if(2230*mm<=mr[0]&&mr[0]<2310*mm) { layer = 6; mat = 0; }
01752 if(2310*mm<=mr[0]&&mr[0]<2350*mm) { layer = 6; mat = 1; }
01753 if(2350*mm<=mr[0]&&mr[0]<2430*mm) { layer = 7; mat = 0; }
01754 if(2430*mm<=mr[0]&&mr[0]<2470*mm) { layer = 7; mat = 1; }
01755 if(2470*mm<=mr[0]&&mr[0]<=2620*mm) { layer = 8; mat = 0; }
01756 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01757 b[0] = -tb[1];
01758 b[1] = tb[0];
01759 b[2] = tb[2];
01760 ifbar = true;
01761 }
01762 if(2050*mm<=mr[2]&&mr[2]<2090*mm&&1034*mm<=mr[0]&&mr[0]<=2500*mm){
01763 part = 0; layer = 0; mat = 0;
01764 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01765 b[0] = -tb[1];
01766 b[1] = tb[0];
01767 b[2] = tb[2];
01768 ifend = true;
01769 }
01770 if(2090*mm<=mr[2]&&mr[2]<2130*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
01771 part = 0; layer = 0; mat = 1;
01772 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01773 b[0] = -tb[1];
01774 b[1] = tb[0];
01775 b[2] = tb[2];
01776 ifend = true;
01777 }
01778 if(2130*mm<=mr[2]&&mr[2]<2170*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
01779 part = 0; layer = 1; mat = 0;
01780 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01781 b[0] = -tb[1];
01782 b[1] = tb[0];
01783 b[2] = tb[2];
01784 ifend = true;
01785 }
01786 if(2170*mm<=mr[2]&&mr[2]<2210*mm&&1100*mm<=mr[0]&&mr[0]<=2500*mm) {
01787 part = 0; layer = 1; mat = 1;
01788 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01789 b[0] = -tb[1];
01790 b[1] = tb[0];
01791 b[2] = tb[2];
01792 ifend = true;
01793 }
01794 if(2210*mm<=mr[2]&&mr[2]<2240*mm&&1100*mm<mr[0]&&mr[0]<=2500*mm) {
01795 part = 0; layer = 2; mat = 0;
01796 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01797 b[0] = -tb[1];
01798 b[1] = tb[0];
01799 b[2] = tb[2];
01800 ifend = true;
01801 }
01802 if(2240*mm<=mr[2]&&mr[2]<2280*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
01803 part = 0; layer = 2; mat = 1;
01804 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01805 b[0] = -tb[1];
01806 b[1] = tb[0];
01807 b[2] = tb[2];
01808 ifend = true;
01809 }
01810 if(2280*mm<=mr[2]&&mr[2]<2310*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
01811 part = 0; layer = 3; mat = 0;
01812 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01813 b[0] = -tb[1];
01814 b[1] = tb[0];
01815 b[2] = tb[2];
01816 ifend = true;
01817 }
01818 if(2310*mm<=mr[2]&&mr[2]<2350*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
01819 part = 0; layer = 3; mat = 1;
01820 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01821 b[0] = -tb[1];
01822 b[1] = tb[0];
01823 b[2] = tb[2];
01824 ifend = true;
01825 }
01826 if(2350*mm<=mr[2]&&mr[2]<2380*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
01827 part = 0; layer = 4; mat = 0;
01828 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01829 b[0] = -tb[1];
01830 b[1] = tb[0];
01831 b[2] = tb[2];
01832 ifend = true;
01833 }
01834 if(2380*mm<=mr[2]&&mr[2]<2420*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
01835 part = 0; layer = 4; mat = 1;
01836 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01837 b[0] = -tb[1];
01838 b[1] = tb[0];
01839 b[2] = tb[2];
01840 ifend = true;
01841 }
01842 if(2420*mm<=mr[2]&&mr[2]<2470*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
01843 part = 0; layer = 5; mat = 0;
01844 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01845 b[0] = -tb[1];
01846 b[1] = tb[0];
01847 b[2] = tb[2];
01848 ifend = true;
01849 }
01850 if(2470*mm<=mr[2]&&mr[2]<2510*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
01851 part = 0; layer = 5; mat =1;
01852 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01853 b[0] = -tb[1];
01854 b[1] = tb[0];
01855 b[2] = tb[2];
01856 ifend = true;
01857 }
01858 if(2510*mm<=mr[2]&&mr[2]<2590*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
01859 part = 0; layer = 6; mat = 0;
01860 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01861 b[0] = -tb[1];
01862 b[1] = tb[0];
01863 b[2] = tb[2];
01864 ifend = true;
01865 }
01866 if(2590*mm<=mr[2]&&mr[2]<2630*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01867 part = 0; layer = 6; mat = 1;
01868 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01869 b[0] = -tb[1];
01870 b[1] = tb[0];
01871 b[2] = tb[2];
01872 ifend = true;
01873 }
01874 if(2630*mm<=mr[2]&&mr[2]<2710*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01875 part = 0; layer = 7; mat = 0;
01876 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01877 b[0] = -tb[1];
01878 b[1] = tb[0];
01879 b[2] = tb[2];
01880 ifend = true;
01881 }
01882 if(2710*mm<=mr[2]&&mr[2]<2750*mm&&1362*mm<=mr[0]&&mr[0]<=2500*mm) {
01883 part = 0; layer = 7; mat = 1;
01884 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01885 b[0] = -tb[1];
01886 b[1] = tb[0];
01887 b[2] = tb[2];
01888 ifend = true;
01889 }
01890 if(2750*mm<=mr[2]&&mr[2]<=2800*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
01891 part = 0; layer = 8; mat = 0;
01892 m_Mucfield->getMucField(part,layer,mat,mr,tb);
01893 b[0] = -tb[1];
01894 b[1] = tb[0];
01895 b[2] = tb[2];
01896 ifend = true;
01897 }
01898 }
01899 if(ifbar==true||ifend==true) {
01900 if( r.x() < 0. && r.y() >= 0. && r.z() > 0. ){
01901 b[0] = -b[0];
01902 }
01903 else if( r.x() <= 0. && r.y() < 0. && r.z() > 0. ){
01904 b[0] = -b[0];
01905 b[1] = -b[1];
01906 }
01907 else if( r.x() > 0. && r.y() < 0. && r.z() > 0. ){
01908 b[1] = -b[1];
01909 }
01910 else if( r.x() >= 0. && r.y() > 0. && r.z() < 0. ){
01911 b[0] = -b[0];
01912 b[1] = -b[1];
01913 }
01914 else if( r.x() < 0. && r.y() >= 0. && r.z() < 0. ){
01915 b[1] = -b[1];
01916 }
01917 else if( r.x() <= 0. && r.y() < 0. && r.z() < 0. ){
01918 b[0] = b[0];
01919 b[1] = b[1];
01920 }
01921 else if( r.x() > 0. && r.y() <= 0. && r.z() < 0. ){
01922 b[0] = -b[0];
01923 }
01924 }
01925 }
01926
01927
01928 newb[0] = -b[0] * m_scale;
01929 newb[1] = b[1] * m_scale;
01930 newb[2] = -b[2] * m_scale;
01931
01932
01933
01934
01935
01936
01937
01938
01939
01940
01941
01942
01943
01944 #ifndef BEAN
01945 return StatusCode::SUCCESS;
01946 #else
01947 return true;
01948 #endif
01949 }
01950
01951 StatusCode MagneticFieldSvc::uniFieldVector(const HepPoint3D& r,
01952 HepVector3D& b) const
01953 {
01954 if(m_runmode == 8 || m_runmode == 7) {
01955 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)
01956 {
01957 b[0]=0.0;
01958 b[1]=0.0;
01959 b[2]=-0.9*tesla;
01960 }
01961 else
01962 {
01963 b[0]=0.0;
01964 b[1]=0.0;
01965 b[2]=0.0;
01966 }
01967 }
01968 else {
01969 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)
01970 {
01971 b[0]=0.0;
01972 b[1]=0.0;
01973 b[2]=-1.0*tesla;
01974 }
01975 else
01976 {
01977 b[0]=0.0;
01978 b[1]=0.0;
01979 b[2]=0.0;
01980 }
01981 }
01982
01983 if(m_turnOffField == true) {
01984 b[0] = 0.;
01985 b[1] = 0.;
01986 b[2] = 0.;
01987 }
01988
01989 b[0] *= m_scale;
01990 b[1] *= m_scale;
01991 b[2] *= m_scale;
01992 #ifndef BEAN
01993 return StatusCode::SUCCESS;
01994 #else
01995 return true;
01996 #endif
01997 }
01998
01999 double MagneticFieldSvc::getReferField()
02000 {
02001 if(m_useDB == false) {
02002 if(m_runmode == 8 || m_runmode == 7) m_zfield = -0.9*tesla;
02003 else m_zfield = -1.0*tesla;
02004 }
02005
02006 if(m_turnOffField == true) {
02007 m_zfield = 0.;
02008 }
02009 return m_zfield * m_scale;
02010 }
02011
02012 bool MagneticFieldSvc::ifRealField() const {
02013 return m_ifRealField;
02014 }
02015
02016
02017
02018
02019 void MagneticFieldSvc::fieldGrid (const HepPoint3D& r,
02020 HepVector3D& bf ) const {
02021
02022 bf[0] = 0.0;
02023 bf[1] = 0.0;
02024 bf[2] = 0.0;
02025
02027 double z = r.z() - m_zOffSet;
02028 if( z < m_min_FL[2] || z > m_max_FL[2] ) return;
02029 double x = r.x();
02030 if( x < m_min_FL[0] || x > m_max_FL[0] ) return;
02031 double y = r.y();
02032 if( y < m_min_FL[1] || y > m_max_FL[1] ) return;
02033 int i = int( (x-m_min_FL[0])/m_Dxyz[0]);
02034 int j = int( (y-m_min_FL[1])/m_Dxyz[1] );
02035 int k = int( (z-m_min_FL[2])/m_Dxyz[2] );
02036
02037 int ijk000 = 3*( m_Nxyz[0]*( m_Nxyz[1]*k + j ) + i );
02038 int ijk001 = 3*( m_Nxyz[0]*( m_Nxyz[1]*(k+1) + j ) + i );
02039 int ijk010 = 3*( m_Nxyz[0]*( m_Nxyz[1]*k + j + 1 ) + i );
02040 int ijk011 = 3*( m_Nxyz[0]*( m_Nxyz[1]*(k+1) + j + 1) + i );
02041 int ijk100 = 3*( m_Nxyz[0]*( m_Nxyz[1]*k + j) + i + 1 );
02042 int ijk101 = 3*( m_Nxyz[0]*( m_Nxyz[1]*(k+1) + j) + i + 1 );
02043 int ijk110 = 3*( m_Nxyz[0]*( m_Nxyz[1]*k + j + 1) + i + 1 );
02044 int ijk111 = 3*( m_Nxyz[0]*( m_Nxyz[1]*(k+1) + j + 1 ) + i + 1 );
02045
02046
02047
02048
02049
02050
02051
02052
02053
02054
02055
02056
02057
02058
02059
02060
02061
02062
02063
02064 double cx000 = m_Q[ ijk000 ];
02065 double cx001 = m_Q[ ijk001 ];
02066 double cx010 = m_Q[ ijk010 ];
02067 double cx011 = m_Q[ ijk011 ];
02068 double cx100 = m_Q[ ijk100 ];
02069 double cx101 = m_Q[ ijk101 ];
02070 double cx110 = m_Q[ ijk110 ];
02071 double cx111 = m_Q[ ijk111 ];
02072 double cy000 = m_Q[ ijk000+1 ];
02073 double cy001 = m_Q[ ijk001+1 ];
02074 double cy010 = m_Q[ ijk010+1 ];
02075 double cy011 = m_Q[ ijk011+1 ];
02076 double cy100 = m_Q[ ijk100+1 ];
02077 double cy101 = m_Q[ ijk101+1 ];
02078 double cy110 = m_Q[ ijk110+1 ];
02079 double cy111 = m_Q[ ijk111+1 ];
02080 double cz000 = m_Q[ ijk000+2 ];
02081 double cz001 = m_Q[ ijk001+2 ];
02082 double cz010 = m_Q[ ijk010+2 ];
02083 double cz011 = m_Q[ ijk011+2 ];
02084 double cz100 = m_Q[ ijk100+2 ];
02085 double cz101 = m_Q[ ijk101+2 ];
02086 double cz110 = m_Q[ ijk110+2 ];
02087 double cz111 = m_Q[ ijk111+2 ];
02088 double hx1 = ( x-m_min_FL[0]-i*m_Dxyz[0] )/m_Dxyz[0];
02089 double hy1 = ( y-m_min_FL[1]-j*m_Dxyz[1] )/m_Dxyz[1];
02090 double hz1 = ( z-m_min_FL[2]-k*m_Dxyz[2] )/m_Dxyz[2];
02091 double hx0 = 1.0-hx1;
02092 double hy0 = 1.0-hy1;
02093 double hz0 = 1.0-hz1;
02094 double h000 = hx0*hy0*hz0;
02095 if( fabs(h000) > 0.0 &&
02096 cx000 > 9.0e5 && cy000 > 9.0e5 && cz000 > 9.0e5) return;
02097
02098 double h001 = hx0*hy0*hz1;
02099 if( fabs(h001) > 0.0 &&
02100 cx001 > 9.0e5 && cy001 > 9.0e5 && cz001 > 9.0e5) return;
02101
02102 double h010 = hx0*hy1*hz0;
02103 if( fabs(h010) > 0.0 &&
02104 cx010 > 9.0e5 && cy010 > 9.0e5 && cz010 > 9.0e5) return;
02105
02106 double h011 = hx0*hy1*hz1;
02107 if( fabs(h011) > 0.0 &&
02108 cx011 > 9.0e5 && cy011 > 9.0e5 && cz011 > 9.0e5) return;
02109
02110 double h100 = hx1*hy0*hz0;
02111 if( fabs(h100) > 0.0 &&
02112 cx100 > 9.0e5 && cy100 > 9.0e5 && cz100 > 9.0e5) return;
02113
02114 double h101 = hx1*hy0*hz1;
02115 if( fabs(h101) > 0.0 &&
02116 cx101 > 9.0e5 && cy101 > 9.0e5 && cz101 > 9.0e5) return;
02117
02118 double h110 = hx1*hy1*hz0;
02119 if( fabs(h110) > 0.0 &&
02120 cx110 > 9.0e5 && cy110 > 9.0e5 && cz110 > 9.0e5) return;
02121
02122 double h111 = hx1*hy1*hz1;
02123 if( fabs(h111) > 0.0 &&
02124 cx111 > 9.0e5 && cy111 > 9.0e5 && cz111 > 9.0e5) return;
02125
02126 bf(0) = ( cx000*h000 + cx001*h001 + cx010*h010 + cx011*h011 +
02127 cx100*h100 + cx101*h101 + cx110*h110 + cx111*h111);
02128 bf(1) = ( cy000*h000 + cy001*h001 + cy010*h010 + cy011*h011 +
02129 cy100*h100 + cy101*h101 + cy110*h110 + cy111*h111 );
02130 bf(2) = ( cz000*h000 + cz001*h001 + cz010*h010 + cz011*h011 +
02131 cz100*h100 + cz101*h101 + cz110*h110 + cz111*h111 );
02132 return;
02133 }
02134
02135
02136
02137
02138 void MagneticFieldSvc::fieldGrid_TE (const HepPoint3D& r,
02139 HepVector3D& bf ) const {
02140
02141 bf[0] = 0.0;
02142 bf[1] = 0.0;
02143 bf[2] = 0.0;
02144
02146 double z = r.z() - m_zOffSet_TE;
02147 if( fabs(z) < m_min_FL_TE[2] || fabs(z) > m_max_FL_TE[2] ) return;
02148 double x = r.x();
02149 if( fabs(x) < m_min_FL_TE[0] || fabs(x) > m_max_FL_TE[0] ) return;
02150 double y = r.y();
02151 if( fabs(y) < m_min_FL_TE[1] || fabs(y) > m_max_FL_TE[1] ) return;
02152 int i = int( (fabs(x)-m_min_FL_TE[0])/m_Dxyz_TE[0]);
02153 int j = int( (fabs(y)-m_min_FL_TE[1])/m_Dxyz_TE[1] );
02154 int k = int( (fabs(z)-m_min_FL_TE[2])/m_Dxyz_TE[2] );
02155
02156 int ijk000 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*k + j ) + i );
02157 int ijk001 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*(k+1) + j ) + i );
02158 int ijk010 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*k + j + 1 ) + i );
02159 int ijk011 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*(k+1) + j + 1) + i );
02160 int ijk100 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*k + j) + i + 1 );
02161 int ijk101 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*(k+1) + j) + i + 1 );
02162 int ijk110 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*k + j + 1) + i + 1 );
02163 int ijk111 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*(k+1) + j + 1 ) + i + 1 );
02164
02165
02166
02167
02168
02169
02170
02171
02172
02173
02174
02175
02176
02177 double cx000 = m_Q_TE[ ijk000 ];
02178 double cx001 = m_Q_TE[ ijk001 ];
02179 double cx010 = m_Q_TE[ ijk010 ];
02180 double cx011 = m_Q_TE[ ijk011 ];
02181 double cx100 = m_Q_TE[ ijk100 ];
02182 double cx101 = m_Q_TE[ ijk101 ];
02183 double cx110 = m_Q_TE[ ijk110 ];
02184 double cx111 = m_Q_TE[ ijk111 ];
02185 double cy000 = m_Q_TE[ ijk000+1 ];
02186 double cy001 = m_Q_TE[ ijk001+1 ];
02187 double cy010 = m_Q_TE[ ijk010+1 ];
02188 double cy011 = m_Q_TE[ ijk011+1 ];
02189 double cy100 = m_Q_TE[ ijk100+1 ];
02190 double cy101 = m_Q_TE[ ijk101+1 ];
02191 double cy110 = m_Q_TE[ ijk110+1 ];
02192 double cy111 = m_Q_TE[ ijk111+1 ];
02193 double cz000 = m_Q_TE[ ijk000+2 ];
02194 double cz001 = m_Q_TE[ ijk001+2 ];
02195 double cz010 = m_Q_TE[ ijk010+2 ];
02196 double cz011 = m_Q_TE[ ijk011+2 ];
02197 double cz100 = m_Q_TE[ ijk100+2 ];
02198 double cz101 = m_Q_TE[ ijk101+2 ];
02199 double cz110 = m_Q_TE[ ijk110+2 ];
02200 double cz111 = m_Q_TE[ ijk111+2 ];
02201 double hx1 = ( fabs(x)-m_min_FL_TE[0]-i*m_Dxyz_TE[0] )/m_Dxyz_TE[0];
02202 double hy1 = ( fabs(y)-m_min_FL_TE[1]-j*m_Dxyz_TE[1] )/m_Dxyz_TE[1];
02203 double hz1 = ( fabs(z)-m_min_FL_TE[2]-k*m_Dxyz_TE[2] )/m_Dxyz_TE[2];
02204 double hx0 = 1.0-hx1;
02205 double hy0 = 1.0-hy1;
02206 double hz0 = 1.0-hz1;
02207 double h000 = hx0*hy0*hz0;
02208 if( fabs(h000) > 0.0 &&
02209 cx000 > 9.0e5 && cy000 > 9.0e5 && cz000 > 9.0e5) return;
02210
02211 double h001 = hx0*hy0*hz1;
02212 if( fabs(h001) > 0.0 &&
02213 cx001 > 9.0e5 && cy001 > 9.0e5 && cz001 > 9.0e5) return;
02214
02215 double h010 = hx0*hy1*hz0;
02216 if( fabs(h010) > 0.0 &&
02217 cx010 > 9.0e5 && cy010 > 9.0e5 && cz010 > 9.0e5) return;
02218
02219 double h011 = hx0*hy1*hz1;
02220 if( fabs(h011) > 0.0 &&
02221 cx011 > 9.0e5 && cy011 > 9.0e5 && cz011 > 9.0e5) return;
02222
02223 double h100 = hx1*hy0*hz0;
02224 if( fabs(h100) > 0.0 &&
02225 cx100 > 9.0e5 && cy100 > 9.0e5 && cz100 > 9.0e5) return;
02226
02227 double h101 = hx1*hy0*hz1;
02228 if( fabs(h101) > 0.0 &&
02229 cx101 > 9.0e5 && cy101 > 9.0e5 && cz101 > 9.0e5) return;
02230
02231 double h110 = hx1*hy1*hz0;
02232 if( fabs(h110) > 0.0 &&
02233 cx110 > 9.0e5 && cy110 > 9.0e5 && cz110 > 9.0e5) return;
02234
02235 double h111 = hx1*hy1*hz1;
02236 if( fabs(h111) > 0.0 &&
02237 cx111 > 9.0e5 && cy111 > 9.0e5 && cz111 > 9.0e5) return;
02238
02239 bf(0) = ( cx000*h000 + cx001*h001 + cx010*h010 + cx011*h011 +
02240 cx100*h100 + cx101*h101 + cx110*h110 + cx111*h111);
02241 bf(1) = ( cy000*h000 + cy001*h001 + cy010*h010 + cy011*h011 +
02242 cy100*h100 + cy101*h101 + cy110*h110 + cy111*h111 );
02243 bf(2) = ( cz000*h000 + cz001*h001 + cz010*h010 + cz011*h011 +
02244 cz100*h100 + cz101*h101 + cz110*h110 + cz111*h111 );
02245
02246
02247 if( r.x() < 0. && r.y() >= 0. && r.z() > 0. ){
02248 bf(0) = -bf(0);
02249 }
02250 else if( r.x() <= 0. && r.y() < 0. && r.z() > 0. ){
02251 bf(0) = -bf(0);
02252 bf(1) = -bf(1);
02253 }
02254 else if( r.x() > 0. && r.y() < 0. && r.z() > 0. ){
02255 bf(1) = -bf(1);
02256 }
02257
02258 else if( r.x() >= 0. && r.y() > 0. && r.z() < 0. ){
02259 bf(0) = -bf(0);
02260 bf(1) = -bf(1);
02261 }
02262 else if( r.x() < 0. && r.y() >= 0. && r.z() < 0. ){
02263 bf(1) = -bf(1);
02264 }
02265 else if( r.x() <= 0. && r.y() < 0. && r.z() < 0. ){
02266 bf(0) = bf(0);
02267 bf(1) = bf(1);
02268 }
02269 else if( r.x() > 0. && r.y() <= 0. && r.z() < 0. ){
02270 bf(0) = -bf(0);
02271 }
02272
02273 return;
02274 }