BOSS 7.1.3
BESIII Offline Software System
Loading...
Searching...
No Matches
MagneticFieldSvc.cxx
Go to the documentation of this file.
1// Include files
2#ifndef BEAN
3#include "GaudiKernel/AlgFactory.h"
4#include "GaudiKernel/SvcFactory.h"
5#include "GaudiKernel/ISvcLocator.h"
6#include "GaudiKernel/MsgStream.h"
7#include "GaudiKernel/Incident.h"
8#include "GaudiKernel/IIncidentSvc.h"
9#include "GaudiKernel/IDataProviderSvc.h"
10#include "GaudiKernel/DataObject.h"
11#include "GaudiKernel/SmartDataPtr.h"
12
13#include "MagneticField/IMagneticFieldSvc.h"
14#endif
15
16#include "MagneticField/MagneticFieldSvc.h"
17#include "MagneticField/MucMagneticField.h"
18
19#include "CLHEP/Units/PhysicalConstants.h"
20
21#ifndef BEAN
22#include "EventModel/EventModel.h"
23#include "EventModel/EventHeader.h"
24#endif
25
26#include <cstring>
27#include <cstdlib>
28#include <fstream>
29#include <iostream>
30using namespace std;
31using namespace CLHEP;
32
33/** @file MagneticFieldSvc.cxx
34 * Implementation of MagneticFieldSvc class
35 *
36 */
37
38#ifndef BEAN
39// Instantiation of a static factory class used by clients to create
40// instances of this service
41//static SvcFactory<MagneticFieldSvc> s_factory;
42//const ISvcFactory& MagneticFieldSvcFactory = s_factory;
43
44//=============================================================================
45// Standard constructor, initializes variables
46//=============================================================================
47DECLARE_COMPONENT(MagneticFieldSvc)
48MagneticFieldSvc::MagneticFieldSvc( const std::string& name,
49 ISvcLocator* svc ) : base_class( name, svc )
50{
51 declareProperty( "TurnOffField", m_turnOffField = false);
52 declareProperty( "FieldMapFile", m_filename );
53 declareProperty( "GridDistance", m_gridDistance = 5);
54 declareProperty( "RunMode", m_runmode = 2);
55 declareProperty( "IfRealField", m_ifRealField = true);
56 declareProperty( "OutLevel", m_outlevel = 1);
57 declareProperty( "Scale", m_scale = 1.0);
58 declareProperty( "UniField", m_uniField = false);
59
60 declareProperty( "Cur_SCQ1_55", m_Cur_SCQ1_55 = 349.4);
61 declareProperty( "Cur_SCQ1_89", m_Cur_SCQ1_89 = 426.2);
62 declareProperty( "Cur_SCQ2_10", m_Cur_SCQ2_10 = 474.2);
63
64 declareProperty( "UseDBFlag", m_useDB = true);
65 declareProperty( "ReadOneTime", m_readOneTime=false);
66 declareProperty( "RunFrom", m_runFrom=8093);
67 declareProperty("RunTo",m_runTo=9025);
68
69 m_Mucfield = new MucMagneticField();
70 if(!m_Mucfield) cout<<"error: can not get MucMagneticField pointer"<<endl;
71
72 m_zfield = -1.0*tesla;
73
74 if(getenv("MAGNETICFIELDROOT") != NULL) {
75 path = std::string(getenv( "MAGNETICFIELDROOT" ));
76 } else {
77 std::cerr<<"Couldn't find MAGNETICFIELDROOT"<<std::endl;
78 }
79
80}
81
82#else
83// -------------------------- BEAN ------------------------------------
84MagneticFieldSvc* MagneticFieldSvc::m_field = 0;
85
86//=============================================================================
87// Standard constructor, initializes variables
88//=============================================================================
89MagneticFieldSvc::MagneticFieldSvc()
90{
91 // use functions instead of "declareProperty"
92 m_gridDistance = 5;
93 m_runmode = 2;
94 m_ifRealField = true;
95 m_outlevel = 1;
96
97 m_Cur_SCQ1_55 = 349.4;
98 m_Cur_SCQ1_89 = 426.2;
99 m_Cur_SCQ2_10 = 474.2;
100
101 m_useDB = true;
102
103 path = "";
104 m_Mucfield = 0;
105
106 m_zfield = -1.0*tesla;
107}
108#endif
109
110//=============================================================================
111// Standard destructor
112//=============================================================================
113MagneticFieldSvc::~MagneticFieldSvc()
114{
115 if(m_Mucfield) delete m_Mucfield;
116}
117
118//=============================================================================
119// Initialize
120//=============================================================================
121#ifdef BEAN
122bool init_mucMagneticField()
123{
124 // reinstall MucMagneticField probably with a new path
125 if( m_Mucfield ) {
126 if( m_Mucfield->getPath() == path ) return true;
127 delete m_Mucfield;
128 }
129
130 m_Mucfield = new(std::nothrow) MucMagneticField(path);
131 if( !m_Mucfield ) {
132 cout<<"error: can not get MucMagneticField pointer"<<endl;
133 return false;
134 }
135
136 return true;
137}
138#endif
139
140StatusCode MagneticFieldSvc::initialize()
141{
142#ifndef BEAN
143 MsgStream log(msgSvc(), name());
144 StatusCode status = Service::initialize();
145 former_m_filename_TE = "First Run";
146 former_m_filename = "First Run";
147 //cout<<"former_m_filename_TE is:"<<former_m_filename_TE<<":former_m_filename is:"<<former_m_filename<<endl;
148 // Get the references to the services that are needed by the ApplicationMgr itself
149 IIncidentSvc* incsvc;
150 status = service("IncidentSvc", incsvc);
151 int priority = 100;
152 if( status.isSuccess() ){
153 incsvc -> addListener(this, "NewRun", priority);
154 }
155
156 status = serviceLocator()->service("EventDataSvc", m_eventSvc, true);
157 if (status.isFailure() ) {
158 log << MSG::ERROR << "Unable to find EventDataSvc " << endreq;
159 return status;
160 }
161#else
162 if( !init_mucMagneticField() ) return false;
163 StatusCode status = true;
164#endif
165 if(m_useDB == false) {
166 if(m_gridDistance == 5) {
167 m_Q.reserve(201250);
168 m_P.reserve(201250);
169 m_Q_TE.reserve(176608);
170 m_P_TE.reserve(176608);
171 }
172 if(m_gridDistance == 10) {
173 m_Q.reserve(27082);
174 m_P.reserve(27082);
175 m_Q_TE.reserve(23833);
176 m_P_TE.reserve(23833);
177 }
178
179 m_filename = path;
180 m_filename_TE = path;
181 if(m_gridDistance == 10) {
182 m_filename_TE += std::string("/dat/TEMap10cm.dat");
183 if(m_runmode == 2)
184 m_filename += std::string( "/dat/2008-04-03BESIII-field-Mode2.dat");
185 else if(m_runmode == 4)
186 m_filename += std::string( "/dat/2008-04-03BESIII-field-Mode4.dat");
187 else if(m_runmode == 6)
188 m_filename += std::string( "/dat/2008-04-03BESIII-field-Mode6.dat");
189 else if(m_runmode == 7)
190 m_filename += std::string( "/dat/2008-04-03BESIII-field-Mode7.dat");
191 else {
192 m_filename += std::string( "/dat/2007-08-28BESIII-field.dat");
193 std::cout<<"WARNING: YOU ARE USING OLD FIELD MAP!!!!"<<std::endl;
194 }
195 }
196 if(m_gridDistance == 5) {
197 m_filename_TE += std::string("/dat/TEMap5cm.dat");
198 if(m_runmode == 1)
199 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode1.dat");
200 else if(m_runmode == 2)
201 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode2.dat");
202 else if(m_runmode == 3)
203 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode3.dat");
204 else if(m_runmode == 4)
205 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode4.dat");
206 else if(m_runmode == 5)
207 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode5.dat");
208 else if(m_runmode == 6)
209 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode6.dat");
210 else if(m_runmode == 7)
211 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode7.dat");
212 else if(m_runmode == 8)
213 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode8.dat");
214 else {
215 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode2.dat");
216 std::cout<<"WARNING: NO RUN MODE, YOU ARE USING DEFAULT FIELD MAP!!!!"<<std::endl;
217 }
218 }
219 cout<<"*** Read field map: "<<m_filename<<endl;
220 cout<<"*** Read field map "<<m_filename_TE<<endl;
221
222 status = parseFile();
223 status = parseFile_TE();
224#ifndef BEAN
225 if ( status.isSuccess() ) {
226 log << MSG::DEBUG << "Magnetic field parsed successfully" << endreq;
227#else
228 if ( status ) {
229 cout << "Magnetic field parsed successfully" << endl;
230#endif
231
232 for (int iC = 0; iC< 2; ++iC ){
233 m_min_FL[iC] = -90.0*cm;
234 m_max_FL[iC] = m_min_FL[iC]+( m_Nxyz[iC]-1 )* m_Dxyz[iC];
235 //for tof and emc
236 m_min_FL_TE[iC] = 0.0*cm;
237 m_max_FL_TE[iC] = m_min_FL_TE[iC]+( m_Nxyz_TE[iC]-1 )* m_Dxyz_TE[iC];
238 } // iC
239 m_min_FL[2] = -120.0*cm;
240 m_max_FL[2] = m_min_FL[2]+( m_Nxyz[2]-1 )* m_Dxyz[2];
241 //for tof and emc
242 m_min_FL_TE[2] = 0.0*cm;
243 m_max_FL_TE[2] = m_min_FL_TE[2]+( m_Nxyz_TE[2]-1 )* m_Dxyz_TE[2];
244 }
245 else {
246#ifndef BEAN
247 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
248#else
249 cout << "Magnetic field parse failled" << endl;
250#endif
251 }
252 }
253 else {
254 m_connect_run = new FieldDBUtil::ConnectionDB();
255 if (m_readOneTime == true){
256 bool e = m_connect_run->getReadSC_MagnetInfo(m_mapMagnetInfo,m_runFrom,m_runTo);
257 e = m_connect_run->getBeamEnergy(m_mapBeamEnergy, m_runFrom,m_runTo);
258 }
259#ifndef BEAN
260 if (!m_connect_run) {
261 log << MSG::ERROR << "Could not open connection to database" << endreq;
262 }
263#endif
264 }
265 return status;
266}
267
268
269#ifndef BEAN
270void MagneticFieldSvc::init_params(std::vector<double> current, std::vector<double> beamEnergy, int runNo)
271{
272 MsgStream log(msgSvc(), name());
273#else
274void init_params(int runNo)
275{
276 if( !init_mucMagneticField() ) {
277 cerr << " STOP! " << endl;
278 exit(1);
279 }
280#endif
281
282 m_Q.clear();
283 m_P.clear();
284 m_Q_1.clear();
285 m_P_1.clear();
286 m_Q_2.clear();
287 m_P_2.clear();
288 m_P_TE.clear();
289 m_Q_TE.clear();
290
291 if(m_gridDistance == 5) {
292 m_Q.reserve(201250);
293 m_P.reserve(201250);
294 m_Q_1.reserve(201250);
295 m_P_1.reserve(201250);
296 m_Q_2.reserve(201250);
297 m_P_2.reserve(201250);
298 m_Q_TE.reserve(176608);
299 m_P_TE.reserve(176608);
300 }
301 if(m_gridDistance == 10) {
302 m_Q.reserve(27082);
303 m_P.reserve(27082);
304 m_Q_1.reserve(27082);
305 m_P_1.reserve(27082);
306 m_Q_2.reserve(27082);
307 m_P_2.reserve(27082);
308 m_Q_TE.reserve(23833);
309 m_P_TE.reserve(23833);
310 }
311
312 char BPR_PRB[200];
313 char BER_PRB[200];
314 sprintf(BPR_PRB,"%f ",beamEnergy[0]);
315 sprintf(BER_PRB,"%f ",beamEnergy[1]);
316 setenv("BEPCII_INFO.BPR_PRB", BPR_PRB, 1);
317 setenv("BEPCII_INFO.BER_PRB", BER_PRB, 1);
318 int ssm_curr = int (current[0]);
319 double scql_curr = current[1];
320 double scqr_curr = current[2];
321// cout<<"curent is:"<<ssm_curr<<"scql_curr is:"<<scql_curr<<"scqr_curr is:"<<scqr_currendl;
322
323
324
325#ifndef BEAN
326 log << MSG::INFO << "SSM current: " << current[0] << " SCQL current: " << current[1] << " SCQR current: " << current[2] << " in Run " << runNo << endreq;
327 log << MSG::INFO << "beamEnergy is:"<< beamEnergy[0] <<" BER_PRB:"<< beamEnergy[1] << " in Run " << runNo << endreq;
328#else
329
330 cout << "SSM current: " << current[0] << " SCQL current: " << current[1]
331 << " SCQR current: " << current[2] << " in Run " << runNo << endl;
332 log << MSG::INFO << "beamEnergy is:"<< beamEnergy[0] <<" BER_PRB:"<< beamEnergy[1] << " in Run " << runNo << endreq;
333#endif
334
335 //int ssm_curr = 3369;
336 //double scql_curr = 426.2;
337 //double scqr_curr = 426.2;
338 //for the energy less than 1.89GeV
339 if(((ssm_curr >= 3367) && (ssm_curr <= 3370) && ((scql_curr+scqr_curr)/2 < m_Cur_SCQ1_89))) {
340 m_zfield = -1.0*tesla;
341 m_filename = path;
342 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode2.dat");
343#ifndef BEAN
344 log << MSG::INFO << "Select field map: " << m_filename << endreq;
345#else
346 cout << "Select field map: " << m_filename << endl;
347#endif
348
349 if(((former_m_filename == "First Run") || (former_m_filename != m_filename))&&(m_readOneTime == true))
350 {
351 former_m_filename = m_filename;
352 }
353 else if(m_readOneTime == false) {
354 }
355 StatusCode status = parseFile();
356#ifndef BEAN
357 if ( !status.isSuccess() ) {
358 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
359 }
360#else
361 if ( !status ) {
362 cout << "Magnetic field parse failled" << endl;
363 }
364#endif
365
366 m_Q_1 = m_Q;
367 m_P_1 = m_P;
368
369 m_Q.clear();
370 m_P.clear();
371 if(m_gridDistance == 5) {
372 m_Q.reserve(201250);
373 m_P.reserve(201250);
374 }
375 if(m_gridDistance == 10) {
376 m_Q.reserve(27082);
377 m_P.reserve(27082);
378 }
379
380 m_filename = path;
381 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode3.dat");
382#ifndef BEAN
383 log << MSG::INFO << "Select field map: " << m_filename << endreq;
384#else
385 cout << "Select field map: " << m_filename << endl;
386#endif
387
388 if(((former_m_filename == "First Run") || (former_m_filename != m_filename))&&(m_readOneTime == true))
389 {
390 former_m_filename = m_filename;
391 }
392 else if(m_readOneTime == false){
393 }
394 status = parseFile();
395#ifndef BEAN
396 if ( !status.isSuccess() ) {
397 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
398 }
399#else
400 if ( !status ) {
401 cout << "Magnetic field parse failled" << endl;
402 }
403#endif
404
405 m_Q_2 = m_Q;
406 m_P_2 = m_P;
407
408 m_Q.clear();
409 m_P.clear();
410 if(m_gridDistance == 5) {
411 m_Q.reserve(201250);
412 m_P.reserve(201250);
413 }
414 if(m_gridDistance == 10) {
415 m_Q.reserve(27082);
416 m_P.reserve(27082);
417 }
418 //check
419 if(m_Q_2.size() != m_Q_1.size()) {
420#ifndef BEAN
421 log << MSG::FATAL << "The two field maps used in this run are wrong!!!" << endreq;
422#else
423 cout << "The two field maps used in this run are wrong!!!" << endl;
424#endif
425 exit(1);
426 }
427
428 for(unsigned int iQ = 0; iQ < m_Q_1.size(); iQ++) {
429 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];
430 m_Q.push_back(fieldvalue);
431 m_P.push_back(m_P_2[iQ]);
432 }
433 }
434 //for the energy greater than 1.89GeV
435 if(((ssm_curr >= 3367) && (ssm_curr <= 3370) && ((scql_curr+scqr_curr)/2 >= m_Cur_SCQ1_89))){
436 m_zfield = -1.0*tesla;
437 m_filename = path;
438 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode3.dat");
439#ifndef BEAN
440 log << MSG::INFO << "Select field map: " << m_filename << endreq;
441#else
442 cout << "Select field map: " << m_filename << endl;
443#endif
444 if(((former_m_filename == "First Run") || (former_m_filename != m_filename))&&(m_readOneTime == true))
445 {
446 former_m_filename = m_filename;
447 }
448 else if(m_readOneTime == false){
449 }
450 StatusCode status = parseFile();
451#ifndef BEAN
452 if ( !status.isSuccess() ) {
453 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
454 }
455#else
456 if ( !status ) {
457 cout << "Magnetic field parse failled" << endl;
458 }
459#endif
460
461 m_Q_1 = m_Q;
462 m_P_1 = m_P;
463
464 m_Q.clear();
465 m_P.clear();
466 if(m_gridDistance == 5) {
467 m_Q.reserve(201250);
468 m_P.reserve(201250);
469 }
470 if(m_gridDistance == 10) {
471 m_Q.reserve(27082);
472 m_P.reserve(27082);
473 }
474
475 m_filename = path;
476 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode4.dat");
477#ifndef BEAN
478 log << MSG::INFO << "Select field map: " << m_filename << endreq;
479#else
480 cout << "Select field map: " << m_filename << endl;
481#endif
482
483 if(((former_m_filename == "First Run") || (former_m_filename != m_filename))&&(m_readOneTime == true))
484 {
485 former_m_filename = m_filename;
486 }
487 else if(m_readOneTime == false){
488 }
489 status = parseFile();
490#ifndef BEAN
491 if ( !status.isSuccess() ) {
492 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
493 }
494#else
495 if ( !status ) {
496 cout << "Magnetic field parse failled" << endl;
497 }
498#endif
499
500 m_Q_2 = m_Q;
501 m_P_2 = m_P;
502
503 m_Q.clear();
504 m_P.clear();
505 if(m_gridDistance == 5) {
506 m_Q.reserve(201250);
507 m_P.reserve(201250);
508 }
509 if(m_gridDistance == 10) {
510 m_Q.reserve(27082);
511 m_P.reserve(27082);
512 }
513
514 //check
515 if(m_Q_2.size() != m_Q_1.size()) {
516#ifndef BEAN
517
518 log << MSG::FATAL << "The two field maps used in this run are wrong!!!" << endreq;
519
520#else
521
522 cout << "The two field maps used in this run are wrong!!!" << endl;
523
524#endif
525 exit(1);
526 }
527
528 for(unsigned int iQ = 0; iQ < m_Q_1.size(); iQ++) {
529 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];
530 m_Q.push_back(fieldvalue);
531 m_P.push_back(m_P_2[iQ]);
532 }
533 }
534 if((ssm_curr >= 3033) && (ssm_curr <= 3035)) {
535 m_zfield = -0.9*tesla;
536 m_filename = path;
537 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode7.dat");
538#ifndef BEAN
539 log << MSG::INFO << "Select field map: " << m_filename << endreq;
540#else
541 cout << "Select field map: " << m_filename << endl;
542#endif
543
544 if(((former_m_filename == "First Run") || (former_m_filename != m_filename))&&(m_readOneTime == true))
545 {
546 former_m_filename = m_filename;
547 }
548 else if(m_readOneTime == false){
549 }
550 StatusCode status = parseFile();
551#ifndef BEAN
552 if ( !status.isSuccess() ) {
553 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
554 }
555#else
556 if ( !status ) {
557 cout << "Magnetic field parse failled" << endl;
558 }
559#endif
560
561 m_Q_1 = m_Q;
562 m_P_1 = m_P;
563
564 m_Q.clear();
565 m_P.clear();
566 if(m_gridDistance == 5) {
567 m_Q.reserve(201250);
568 m_P.reserve(201250);
569 }
570 if(m_gridDistance == 10) {
571 m_Q.reserve(27082);
572 m_P.reserve(27082);
573 }
574
575 m_filename = path;
576 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode8.dat");
577#ifndef BEAN
578 log << MSG::INFO << "Select field map: " << m_filename << endreq;
579#else
580 cout << "Select field map: " << m_filename << endl;
581#endif
582
583 if(((former_m_filename == "First Run") || (former_m_filename != m_filename))&&(m_readOneTime == true)){
584 former_m_filename = m_filename;
585 }
586 else if(m_readOneTime == false){
587 }
588 status = parseFile();
589#ifndef BEAN
590 if ( !status.isSuccess() ) {
591 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
592 }
593#else
594 if ( !status ) {
595 cout << "Magnetic field parse failled" << endl;
596 }
597#endif
598
599 m_Q_2 = m_Q;
600 m_P_2 = m_P;
601
602 m_Q.clear();
603 m_P.clear();
604 if(m_gridDistance == 5) {
605 m_Q.reserve(201250);
606 m_P.reserve(201250);
607 }
608 if(m_gridDistance == 10) {
609 m_Q.reserve(27082);
610 m_P.reserve(27082);
611 }
612 //check
613 if(m_Q_2.size() != m_Q_1.size()) {
614#ifndef BEAN
615 log << MSG::FATAL << "The two field maps used in this run are wrong!!!" << endreq;
616#else
617 cout << "The two field maps used in this run are wrong!!!" << endl;
618#endif
619 exit(1);
620 }
621
622 for(unsigned int iQ = 0; iQ < m_Q_1.size(); iQ++) {
623 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];
624 m_Q.push_back(fieldvalue);
625 m_P.push_back(m_P_2[iQ]);
626 }
627 }
628 if((!((ssm_curr >= 3367) && (ssm_curr <= 3370)) && !((ssm_curr >= 3033) && (ssm_curr <= 3035))) || scql_curr == 0 || scqr_curr == 0) {
629#ifndef BEAN
630 log << MSG::FATAL << "The current of run " << runNo << " is abnormal in database, please check it! " << endreq;
631 log << MSG::FATAL << "The current of SSM is " << ssm_curr << ", SCQR is " << scqr_curr << ", SCQL is " << scql_curr << endreq;
632#else
633 cout << "The current of run " << runNo
634 << " is abnormal in database, please check it! " << endl;
635 cout << "The current of SSM is " << ssm_curr
636 << ", SCQR is " << scqr_curr << ", SCQL is " << scql_curr << endl;
637#endif
638 exit(1);
639
640 }
641 if(m_Q.size() == 0) {
642#ifndef BEAN
643 log << MSG::FATAL << "This size of field map vector is ZERO, please check the current of run " << runNo << " in database!" << endreq;
644#else
645 cout << "This size of field map vector is ZERO,"
646 << " please check the current of run " << runNo
647 << " in database!" << endl;
648#endif
649 exit(1);
650 }
651
652 m_filename_TE = path;
653 if(m_gridDistance == 10) m_filename_TE += std::string( "/dat/TEMap10cm.dat");
654 if(m_gridDistance == 5) m_filename_TE += std::string( "/dat/TEMap5cm.dat");
655#ifndef BEAN
656 log << MSG::INFO << "Select field map: " << m_filename_TE << endreq;
657#else
658 cout << "Select field map: " << m_filename_TE << endl;
659#endif
660 if(((former_m_filename_TE == "First Run") || (former_m_filename_TE != m_filename_TE))&&(m_readOneTime == true))
661 {
662 former_m_filename_TE = m_filename_TE;
663 }
664 else if (m_readOneTime == false)
665 {}
666 StatusCode status = parseFile_TE();
667#ifndef BEAN
668 if ( !status.isSuccess() ) {
669 log << MSG::DEBUG << "Magnetic field parse failled" << endreq;
670 }
671#else
672 if ( !status ) {
673 cout << "Magnetic field parse failled" << endl;
674 }
675#endif
676
677 for (int iC = 0; iC< 2; ++iC ){
678 m_min_FL[iC] = -90.0*cm;
679 m_max_FL[iC] = m_min_FL[iC]+( m_Nxyz[iC]-1 )* m_Dxyz[iC];
680 //for tof and emc
681 m_min_FL_TE[iC] = 0.0*cm;
682 m_max_FL_TE[iC] = m_min_FL_TE[iC]+( m_Nxyz_TE[iC]-1 )* m_Dxyz_TE[iC];
683 } // iC
684 m_min_FL[2] = -120.0*cm;
685 m_max_FL[2] = m_min_FL[2]+( m_Nxyz[2]-1 )* m_Dxyz[2];
686 //for tof and emc
687 m_min_FL_TE[2] = 0.0*cm;
688 m_max_FL_TE[2] = m_min_FL_TE[2]+( m_Nxyz_TE[2]-1 )* m_Dxyz_TE[2];
689}
690
691#ifndef BEAN
692//=============================================================================
693// Finalize
694//=============================================================================
695StatusCode MagneticFieldSvc::finalize()
696{
697 MsgStream log( msgSvc(), name() );
698 //if(m_connect_run) delete m_connect_run;
699 StatusCode status = Service::finalize();
700
701 if ( status.isSuccess() )
702 log << MSG::INFO << "Service finalized successfully" << endreq;
703 return status;
704}
705
706//=============================================================================
707// QueryInterface
708//=============================================================================
709/*StatusCode MagneticFieldSvc::queryInterface( const InterfaceID& riid,
710 void** ppvInterface )
711{
712 StatusCode sc = StatusCode::FAILURE;
713 if ( ppvInterface ) {
714 *ppvInterface = 0;
715
716 if ( riid == IID_IMagneticFieldSvc ) {
717 *ppvInterface = static_cast<IMagneticFieldSvc*>(this);
718 sc = StatusCode::SUCCESS;
719 addRef();
720 }
721 else
722 sc = Service::queryInterface( riid, ppvInterface );
723 }
724 return sc;
725}
726*/
727void MagneticFieldSvc::handle(const Incident& inc) {
728 MsgStream log( messageService(), name() );
729 log << MSG::DEBUG << "handle: " << inc.type() << endreq;
730 if ( inc.type() != "NewRun" ){
731 return;
732 }
733 log << MSG::DEBUG << "Begin New Runcc" << endreq;
734
735 SmartDataPtr<Event::EventHeader> eventHeader(m_eventSvc,"/Event/EventHeader");
736 int new_run = eventHeader->runNumber();
737 if(new_run<0) new_run=-new_run;
738
739
740 if((m_useDB == true)&&(m_readOneTime == false)) {
741 using FieldDBUtil::ConnectionDB;
742 ConnectionDB::eRet e = m_connect_run->getReadSC_MagnetInfo(current,std::abs(new_run));
743 e = m_connect_run->getBeamEnergy(beamEnergy, std::abs(new_run));
744 init_params(current,beamEnergy,new_run);
745 }
746 else if ((m_useDB == true)&&(m_readOneTime == true))
747 {
748 std::vector<double> beamEnergy;
749 std::vector<double> current(3,0.);
750 //if (m_mapBeamEnergy == NULL ||m_mapBeamEnergy.size() == 0))
751 if (m_mapBeamEnergy.size() == 0)
752 {
753 cout<<"Run:"<<new_run<<" BeamEnergy is NULL, please check!!"<<endl;
754 exit(1);
755 }
756 beamEnergy.push_back((m_mapBeamEnergy[new_run])[0]);
757 beamEnergy.push_back((m_mapBeamEnergy[new_run])[1]);
758 //if (m_mapMagnetInfo.size()<1||m_mapMagnetInfo.isEmpty())
759 if(m_mapMagnetInfo.size() == 0)
760 {
761 cout<<"Run:"<<new_run<<" MagnetInfo is NULL, please check!!"<<endl;
762 exit(1);
763 }
764 current[0]=(m_mapMagnetInfo[new_run])[0];
765 current[1]=(m_mapMagnetInfo[new_run])[1];
766 current[2]=(m_mapMagnetInfo[new_run])[2];
767 init_params(current,beamEnergy,new_run);
768 }
769}
770
771#else
772void MagneticFieldSvc::handle(int new_run) {
773 static int save_run = 0;
774 if( new_run == save_run ) return;
775
776 cout << "Begin New Run " << new_run << endl;
777 save_run = new_run;
778 if(m_useDB == true) {
779 init_params(current,beamEnergy,new_run);
780 }
781}
782#endif
783
784// ---------------------------------------------------------------------------
785// Routine: parseFile
786// Purpose: Parses the file and fill a magnetic field vector
787// ---------------------------------------------------------------------------
788StatusCode MagneticFieldSvc::parseFile() {
789#ifndef BEAN
790 StatusCode sc = StatusCode::FAILURE;
791
792 MsgStream log( msgSvc(), name() );
793#else
794 StatusCode sc = false;
795#endif
796
797 char line[ 255 ];
798 std::ifstream infile( m_filename.c_str() );
799 if ( infile ) {
800#ifndef BEAN
801 sc = StatusCode::SUCCESS;
802#else
803 sc = true;
804#endif
805
806 // Skip the header till PARAMETER
807 do{
808 infile.getline( line, 255 );
809 } while( line[0] != 'P' );
810
811 // Get the PARAMETER
812 std::string sPar[2];
813 char* token = strtok( line, " " );
814 int ip = 0;
815 do{
816 if ( token ) { sPar[ip] = token; token = strtok( NULL, " " );}
817 else continue;
818 ip++;
819 } while ( token != NULL );
820 long int npar = atoi( sPar[1].c_str() );
821
822 // Skip the header till GEOMETRY
823 do{
824 infile.getline( line, 255 );
825 } while( line[0] != 'G' );
826
827 // Skip any comment before GEOMETRY
828 do{
829 infile.getline( line, 255 );
830 } while( line[0] != '#' );
831
832 // Get the GEOMETRY
833 infile.getline( line, 255 );
834 std::string sGeom[7];
835 token = strtok( line, " " );
836 int ig = 0;
837 do{
838 if ( token ) { sGeom[ig] = token; token = strtok( NULL, " " );}
839 else continue;
840 ig++;
841 } while (token != NULL);
842
843 // Grid dimensions are given in cm in CDF file. Convert to CLHEP units
844 m_Dxyz[0] = atof( sGeom[0].c_str() ) * cm;
845 m_Dxyz[1] = atof( sGeom[1].c_str() ) * cm;
846 m_Dxyz[2] = atof( sGeom[2].c_str() ) * cm;
847 m_Nxyz[0] = atoi( sGeom[3].c_str() );
848 m_Nxyz[1] = atoi( sGeom[4].c_str() );
849 m_Nxyz[2] = atoi( sGeom[5].c_str() );
850 m_zOffSet = atof( sGeom[6].c_str() ) * cm;
851 // Number of lines with data to be read
852 long int nlines = ( npar - 7 ) / 3;
853
854 // Check number of lines with data read in the loop
855 long int ncheck = 0;
856
857 // Skip comments and fill a vector of magnetic components for the
858 // x, y and z positions given in GEOMETRY
859
860 while( infile ) {
861 // parse each line of the file,
862 // comment lines begin with '#' in the cdf file
863 infile.getline( line, 255 );
864 if ( line[0] == '#' ) continue;
865 std::string gridx, gridy, gridz, sFx, sFy, sFz;
866 char* token = strtok( line, " " );
867 if ( token ) { gridx = token; token = strtok( NULL, " " );} else continue;
868 if ( token ) { gridy = token; token = strtok( NULL, " " );} else continue;
869 if ( token ) { gridz = token; token = strtok( NULL, " " );} else continue;
870 if ( token ) { sFx = token; token = strtok( NULL, " " );} else continue;
871 if ( token ) { sFy = token; token = strtok( NULL, " " );} else continue;
872 if ( token ) { sFz = token; token = strtok( NULL, " " );} else continue;
873 if ( token != NULL ) continue;
874 //Grid position
875 double gx = atof( gridx.c_str() ) * m;
876 double gy = atof( gridy.c_str() ) * m;
877 double gz = atof( gridz.c_str() ) * m;
878 // Field values are given in tesla in CDF file. Convert to CLHEP units
879 double fx = atof( sFx.c_str() ) * tesla;
880 double fy = atof( sFy.c_str() ) * tesla;
881 double fz = atof( sFz.c_str() ) * tesla;
882 //for debug
883 if(m_outlevel == 0) {
884#ifndef BEAN
885 log << MSG::DEBUG << "grid x, y, z is "<< gx << ", " << gy << ", " << gz << endreq;
886 log << MSG::DEBUG << "field x, y, z is "<< fx << ", " << fy << ", " << fz << endreq;
887#else
888 cout << "grid x, y, z is "<< gx << ", " << gy << ", " << gz << endl;
889 cout << "field x, y, z is "<< fx << ", " << fy << ", " << fz << endl;
890#endif
891 }
892 //Fill the postion to the vector
893 m_P.push_back( gx );
894 m_P.push_back( gy );
895 m_P.push_back( gz );
896 // Add the magnetic field components of each point to
897 // sequentialy in a vector
898 m_Q.push_back( fx );
899 m_Q.push_back( fy );
900 m_Q.push_back( fz );
901 // counts after reading and filling to match the number of lines
902 ncheck++;
903 }
904 infile.close();
905 if ( nlines != ncheck) {
906#ifndef BEAN
907 log << MSG::ERROR << "nline is " << nlines << "; ncheck is " << ncheck << " Number of points in field map does not match!"
908 << endreq;
909 return StatusCode::FAILURE;
910#else
911 cout << "nline is " << nlines << "; ncheck is " << ncheck << " Number of points in field map does not match!"
912 << endl;
913 return false;
914#endif
915 }
916 }
917 else {
918#ifndef BEAN
919 log << MSG::ERROR << "Unable to open magnetic field file : "
920 << m_filename << endreq;
921#else
922 cout << "Unable to open magnetic field file : "
923 << m_filename << endl;
924#endif
925 }
926
927 return sc;
928}
929
930
931// ---------------------------------------------------------------------------
932// Routine: parseFile_TE
933// Purpose: Parses the file and fill a magnetic field vector for tof and emc
934// ---------------------------------------------------------------------------
935StatusCode MagneticFieldSvc::parseFile_TE() {
936#ifndef BEAN
937 StatusCode sc = StatusCode::FAILURE;
938
939 MsgStream log( msgSvc(), name() );
940#else
941 StatusCode sc = false;
942#endif
943
944 char line[ 255 ];
945 std::ifstream infile( m_filename_TE.c_str() );
946
947 if ( infile ) {
948#ifndef BEAN
949 sc = StatusCode::SUCCESS;
950#else
951 sc = true;
952#endif
953
954 // Skip the header till PARAMETER
955 do{
956 infile.getline( line, 255 );
957 } while( line[0] != 'P' );
958
959 // Get the PARAMETER
960 std::string sPar[2];
961 char* token = strtok( line, " " );
962 int ip = 0;
963 do{
964 if ( token ) { sPar[ip] = token; token = strtok( NULL, " " );}
965 else continue;
966 ip++;
967 } while ( token != NULL );
968 long int npar = atoi( sPar[1].c_str() );
969
970 // Skip the header till GEOMETRY
971 do{
972 infile.getline( line, 255 );
973 } while( line[0] != 'G' );
974
975 // Skip any comment before GEOMETRY
976 do{
977 infile.getline( line, 255 );
978 } while( line[0] != '#' );
979
980 // Get the GEOMETRY
981 infile.getline( line, 255 );
982 std::string sGeom[7];
983 token = strtok( line, " " );
984 int ig = 0;
985 do{
986 if ( token ) { sGeom[ig] = token; token = strtok( NULL, " " );}
987 else continue;
988 ig++;
989 } while (token != NULL);
990
991 // Grid dimensions are given in cm in CDF file. Convert to CLHEP units
992 m_Dxyz_TE[0] = atof( sGeom[0].c_str() ) * cm;
993 m_Dxyz_TE[1] = atof( sGeom[1].c_str() ) * cm;
994 m_Dxyz_TE[2] = atof( sGeom[2].c_str() ) * cm;
995 m_Nxyz_TE[0] = atoi( sGeom[3].c_str() );
996 m_Nxyz_TE[1] = atoi( sGeom[4].c_str() );
997 m_Nxyz_TE[2] = atoi( sGeom[5].c_str() );
998 m_zOffSet_TE = atof( sGeom[6].c_str() ) * cm;
999 // Number of lines with data to be read
1000 long int nlines = ( npar - 7 ) / 3;
1001
1002 // Check number of lines with data read in the loop
1003 long int ncheck = 0;
1004
1005 // Skip comments and fill a vector of magnetic components for the
1006 // x, y and z positions given in GEOMETRY
1007
1008 while( infile ) {
1009 // parse each line of the file,
1010 // comment lines begin with '#' in the cdf file
1011 infile.getline( line, 255 );
1012 if ( line[0] == '#' ) continue;
1013 std::string gridx, gridy, gridz, sFx, sFy, sFz;
1014 char* token = strtok( line, " " );
1015 if ( token ) { gridx = token; token = strtok( NULL, " " );} else continue;
1016 if ( token ) { gridy = token; token = strtok( NULL, " " );} else continue;
1017 if ( token ) { gridz = token; token = strtok( NULL, " " );} else continue;
1018 if ( token ) { sFx = token; token = strtok( NULL, " " );} else continue;
1019 if ( token ) { sFy = token; token = strtok( NULL, " " );} else continue;
1020 if ( token ) { sFz = token; token = strtok( NULL, " " );} else continue;
1021 if ( token != NULL ) continue;
1022 //Grid position
1023 double gx = atof( gridx.c_str() ) * m;
1024 double gy = atof( gridy.c_str() ) * m;
1025 double gz = atof( gridz.c_str() ) * m;
1026 // Field values are given in tesla in CDF file. Convert to CLHEP units
1027 double fx = atof( sFx.c_str() ) * tesla;
1028 double fy = atof( sFy.c_str() ) * tesla;
1029 double fz = atof( sFz.c_str() ) * tesla;
1030 //for debug
1031 if(m_outlevel == 0) {
1032#ifndef BEAN
1033 log << MSG::DEBUG << "grid x, y, z is "<< gx << ", " << gy << ", " << gz << endreq;
1034 log << MSG::DEBUG << "field x, y, z is "<< fx << ", " << fy << ", " << fz << endreq;
1035#else
1036 cout << "grid x, y, z is "<< gx << ", " << gy << ", " << gz << endl;
1037 cout << "field x, y, z is "<< fx << ", " << fy << ", " << fz << endl;
1038#endif
1039 }
1040 //Fill the postion to the vector
1041 m_P_TE.push_back( gx );
1042 m_P_TE.push_back( gy );
1043 m_P_TE.push_back( gz );
1044 // Add the magnetic field components of each point to
1045 // sequentialy in a vector
1046 m_Q_TE.push_back( fx );
1047 m_Q_TE.push_back( fy );
1048 m_Q_TE.push_back( fz );
1049 // counts after reading and filling to match the number of lines
1050 ncheck++;
1051 }
1052 infile.close();
1053 if ( nlines != ncheck) {
1054#ifndef BEAN
1055 log << MSG::ERROR << "nline is " << nlines << "; ncheck is " << ncheck << " Number of points in field map does not match!"
1056 << endreq;
1057 return StatusCode::FAILURE;
1058#else
1059 cout << "nline is " << nlines << "; ncheck is " << ncheck << " Number of points in field map does not match!"
1060 << endl;
1061 return false;
1062#endif
1063 }
1064 }
1065 else {
1066#ifndef BEAN
1067 log << MSG::ERROR << "Unable to open magnetic field file : "
1068 << m_filename_TE << endreq;
1069#else
1070 cout << "Unable to open magnetic field file : "
1071 << m_filename_TE << endl;
1072#endif
1073 }
1074
1075 return sc;
1076}
1077//=============================================================================
1078// FieldVector: find the magnetic field value at a given point in space
1079//=============================================================================
1080StatusCode MagneticFieldSvc::fieldVector(const HepPoint3D& newr,
1081 HepVector3D& newb) const
1082{
1083
1084 if(m_turnOffField == true) {
1085 newb[0] = 0.;
1086 newb[1] = 0.;
1087 newb[2] = 0.;
1088#ifndef BEAN
1089 return StatusCode::SUCCESS;
1090#else
1091 return true;
1092#endif
1093 }
1094
1095 // wll added 2012-08-27
1096 if(m_uniField) {
1097 uniFieldVector(newr,newb);
1098#ifndef BEAN
1099 return StatusCode::SUCCESS;
1100#else
1101 return true;
1102#endif
1103 }
1104
1105
1106 //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
1107 double new_x = -newr.x();
1108 double new_y = newr.y();
1109 double new_z = -newr.z();
1110 HepPoint3D r(new_x,new_y,new_z);
1111 HepVector3D b;
1112 b[0]=0.0*tesla;
1113 b[1]=0.0*tesla;
1114 b[2]=0.0*tesla;
1115 // This routine is now dummy. Was previously converting to/from CLHEP units
1116 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)
1117 {
1118 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)
1119 //if(std::abs(r.z())<1.2*m&&std::abs(r.x())<=0.9*m&&std::abs(r.y())<=0.9*m)
1120 {
1121 this->fieldGrid( r, b );
1122 }
1123 else
1124 {
1125 this->fieldGrid_TE( r, b );
1126 }
1127 }
1128 if((fabs(r.z())<=1970*mm && sqrt(r.x()*r.x()+r.y()*r.y()) >= 1740*mm) || (fabs(r.z())>=2050*mm))
1129 {
1130 HepPoint3D mr;
1131 HepVector3D tb;
1132 int part = 0, layer = 0, mat = 0;
1133 double theta;
1134 bool ifbar = false;
1135 bool ifend = false;
1136 if(r.x()!=0.){
1137 theta = atan2(fabs(r.y()),fabs(r.x()));
1138 if(0<=theta&&theta<pi/8) {
1139 mr[0] = fabs(r.x()); mr[1] = -fabs(r.y()); mr[2] = fabs(r.z());
1140 if(mr[2]<=1970*mm&&1740*mm<=mr[0]&&mr[0]<=2620*mm){
1141 part = 1;
1142 if(1740*mm<=mr[0]&&mr[0]<1770*mm) { layer = 0; mat = 0; }
1143 if(1770*mm<=mr[0]&&mr[0]<1810*mm) { layer = 0; mat = 1; }
1144 if(1810*mm<=mr[0]&&mr[0]<1840*mm) { layer = 1; mat = 0; }
1145 if(1840*mm<=mr[0]&&mr[0]<1880*mm) { layer = 1; mat = 1; }
1146 if(1880*mm<=mr[0]&&mr[0]<1910*mm) { layer = 2; mat = 0; }
1147 if(1910*mm<=mr[0]&&mr[0]<1950*mm) { layer = 2; mat = 1; }
1148 if(1950*mm<=mr[0]&&mr[0]<1990*mm) { layer = 3; mat = 0; }
1149 if(1990*mm<=mr[0]&&mr[0]<2030*mm) { layer = 3; mat = 1; }
1150 if(2030*mm<=mr[0]&&mr[0]<2070*mm) { layer = 4; mat = 0; }
1151 if(2070*mm<=mr[0]&&mr[0]<2110*mm) { layer = 4; mat = 1; }
1152 if(2110*mm<=mr[0]&&mr[0]<2190*mm) { layer = 5; mat = 0; }
1153 if(2190*mm<=mr[0]&&mr[0]<2230*mm) { layer = 5; mat = 1; }
1154 if(2230*mm<=mr[0]&&mr[0]<2310*mm) { layer = 6; mat = 0; }
1155 if(2310*mm<=mr[0]&&mr[0]<2350*mm) { layer = 6; mat = 1; }
1156 if(2350*mm<=mr[0]&&mr[0]<2430*mm) { layer = 7; mat = 0; }
1157 if(2430*mm<=mr[0]&&mr[0]<2470*mm) { layer = 7; mat = 1; }
1158 if(2470*mm<=mr[0]&&mr[0]<=2620*mm) { layer = 8; mat = 0; }
1159 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1160 b[0] = tb[0];
1161 b[1] = -tb[1];
1162 b[2] = tb[2];
1163 ifbar = true;
1164 }
1165 if(2050*mm<=mr[2]&&mr[2]<2090*mm&&1034*mm<=mr[0]&&mr[0]<=2500*mm){
1166 part = 0; layer = 0; mat = 0;
1167 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1168 b[0] = tb[0];
1169 b[1] = -tb[1];
1170 b[2] = tb[2];
1171 ifend = true;
1172 }
1173 if(2090*mm<=mr[2]&&mr[2]<2130*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
1174 part = 0; layer = 0; mat = 1;
1175 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1176 b[0] = tb[0];
1177 b[1] = -tb[1];
1178 b[2] = tb[2];
1179 ifend = true;
1180 }
1181 if(2130*mm<=mr[2]&&mr[2]<2170*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
1182 part = 0; layer = 1; mat = 0;
1183 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1184 b[0] = tb[0];
1185 b[1] = -tb[1];
1186 b[2] = tb[2];
1187 ifend = true;
1188 }
1189 if(2170*mm<=mr[2]&&mr[2]<2210*mm&&1100*mm<=mr[0]&&mr[0]<=2500*mm) {
1190 part = 0; layer = 1; mat = 1;
1191 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1192 b[0] = tb[0];
1193 b[1] = -tb[1];
1194 b[2] = tb[2];
1195 ifend = true;
1196 }
1197 if(2210*mm<=mr[2]&&mr[2]<2240*mm&&1100*mm<mr[0]&&mr[0]<=2500*mm) {
1198 part = 0; layer = 2; mat = 0;
1199 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1200 b[0] = tb[0];
1201 b[1] = -tb[1];
1202 b[2] = tb[2];
1203 ifend = true;
1204 }
1205 if(2240*mm<=mr[2]&&mr[2]<2280*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
1206 part = 0; layer = 2; mat = 1;
1207 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1208 b[0] = tb[0];
1209 b[1] = -tb[1];
1210 b[2] = tb[2];
1211 ifend = true;
1212 }
1213 if(2280*mm<=mr[2]&&mr[2]<2310*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
1214 part = 0; layer = 3; mat = 0;
1215 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1216 b[0] = tb[0];
1217 b[1] = -tb[1];
1218 b[2] = tb[2];
1219 ifend = true;
1220 }
1221 if(2310*mm<=mr[2]&&mr[2]<2350*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
1222 part = 0; layer = 3; mat = 1;
1223 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1224 b[0] = tb[0];
1225 b[1] = -tb[1];
1226 b[2] = tb[2];
1227 ifend = true;
1228 }
1229 if(2350*mm<=mr[2]&&mr[2]<2380*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
1230 part = 0; layer = 4; mat = 0;
1231 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1232 b[0] = tb[0];
1233 b[1] = -tb[1];
1234 b[2] = tb[2];
1235 ifend = true;
1236 }
1237 if(2380*mm<=mr[2]&&mr[2]<2420*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
1238 part = 0; layer = 4; mat = 1;
1239 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1240 b[0] = tb[0];
1241 b[1] = -tb[1];
1242 b[2] = tb[2];
1243 ifend = true;
1244 }
1245 if(2420*mm<=mr[2]&&mr[2]<2470*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
1246 part = 0; layer = 5; mat = 0;
1247 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1248 b[0] = tb[0];
1249 b[1] = -tb[1];
1250 b[2] = tb[2];
1251 ifend = true;
1252 }
1253 if(2470*mm<=mr[2]&&mr[2]<2510*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
1254 part = 0; layer = 5; mat =1;
1255 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1256 b[0] = tb[0];
1257 b[1] = -tb[1];
1258 b[2] = tb[2];
1259 ifend = true;
1260 }
1261 if(2510*mm<=mr[2]&&mr[2]<2590*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
1262 part = 0; layer = 6; mat = 0;
1263 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1264 b[0] = tb[0];
1265 b[1] = -tb[1];
1266 b[2] = tb[2];
1267 ifend = true;
1268 }
1269 if(2590*mm<=mr[2]&&mr[2]<2630*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
1270 part = 0; layer = 6; mat = 1;
1271 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1272 b[0] = tb[0];
1273 b[1] = -tb[1];
1274 b[2] = tb[2];
1275 ifend = true;
1276 }
1277 if(2630*mm<=mr[2]&&mr[2]<2710*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
1278 part = 0; layer = 7; mat = 0;
1279 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1280 b[0] = tb[0];
1281 b[1] = -tb[1];
1282 b[2] = tb[2];
1283 ifend = true;
1284 }
1285 if(2710*mm<=mr[2]&&mr[2]<2750*mm&&1362*mm<=mr[0]&&mr[0]<=2500*mm) {
1286 part = 0; layer = 7; mat = 1;
1287 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1288 b[0] = tb[0];
1289 b[1] = -tb[1];
1290 b[2] = tb[2];
1291 ifend = true;
1292 }
1293 if(2750*mm<=mr[2]&&mr[2]<=2800*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
1294 part = 0; layer = 8; mat = 0;
1295 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1296 b[0] = tb[0];
1297 b[1] = -tb[1];
1298 b[2] = tb[2];
1299 ifend = true;
1300 }
1301 }
1302 if(pi/8<=theta&&theta<pi/4) {
1303 mr[0] = fabs(r.x())*cos(pi/4)+fabs(r.y())*sin(pi/4);
1304 mr[1] = -fabs(r.x())*sin(pi/4)+fabs(r.y())*cos(pi/4);
1305 mr[2] = fabs(r.z());
1306 if(mr[2]<=1970*mm&&1740*mm<=mr[0]&&mr[0]<=2620*mm) {
1307 part = 1;
1308 if(1740*mm<=mr[0]&&mr[0]<1770*mm) { layer = 0; mat = 0; }
1309 if(1770*mm<=mr[0]&&mr[0]<1810*mm) { layer = 0; mat = 1; }
1310 if(1810*mm<=mr[0]&&mr[0]<1840*mm) { layer = 1; mat = 0; }
1311 if(1840*mm<=mr[0]&&mr[0]<1880*mm) { layer = 1; mat = 1; }
1312 if(1880*mm<=mr[0]&&mr[0]<1910*mm) { layer = 2; mat = 0; }
1313 if(1910*mm<=mr[0]&&mr[0]<1950*mm) { layer = 2; mat = 1; }
1314 if(1950*mm<=mr[0]&&mr[0]<1990*mm) { layer = 3; mat = 0; }
1315 if(1990*mm<=mr[0]&&mr[0]<2030*mm) { layer = 3; mat = 1; }
1316 if(2030*mm<=mr[0]&&mr[0]<2070*mm) { layer = 4; mat = 0; }
1317 if(2070*mm<=mr[0]&&mr[0]<2110*mm) { layer = 4; mat = 1; }
1318 if(2110*mm<=mr[0]&&mr[0]<2190*mm) { layer = 5; mat = 0; }
1319 if(2190*mm<=mr[0]&&mr[0]<2230*mm) { layer = 5; mat = 1; }
1320 if(2230*mm<=mr[0]&&mr[0]<2310*mm) { layer = 6; mat = 0; }
1321 if(2310*mm<=mr[0]&&mr[0]<2350*mm) { layer = 6; mat = 1; }
1322 if(2350*mm<=mr[0]&&mr[0]<2430*mm) { layer = 7; mat = 0; }
1323 if(2430*mm<=mr[0]&&mr[0]<2470*mm) { layer = 7; mat = 1; }
1324 if(2470*mm<=mr[0]&&mr[0]<=2620*mm) { layer = 8; mat = 0; }
1325 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1326 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
1327 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
1328 b[2] = tb[2];
1329 ifbar = true;
1330 }
1331 if(2050*mm<=mr[2]&&mr[2]<2090*mm&&1034*mm<=mr[0]&&mr[0]<=2500*mm){
1332 part = 0; layer = 0; mat = 0;
1333 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1334 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
1335 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
1336 b[2] = tb[2];
1337 ifend = true;
1338 }
1339 if(2090*mm<=mr[2]&&mr[2]<2130*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
1340 part = 0; layer = 0; mat = 1;
1341 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1342 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
1343 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
1344 b[2] = tb[2];
1345 ifend = true;
1346 }
1347 if(2130*mm<=mr[2]&&mr[2]<2170*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
1348 part = 0; layer = 1; mat = 0;
1349 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1350 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
1351 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
1352 b[2] = tb[2];
1353 ifend = true;
1354 }
1355 if(2170*mm<=mr[2]&&mr[2]<2210*mm&&1100*mm<=mr[0]&&mr[0]<=2500*mm) {
1356 part = 0; layer = 1; mat = 1;
1357 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1358 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
1359 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
1360 b[2] = tb[2];
1361 ifend = true;
1362 }
1363 if(2210*mm<=mr[2]&&mr[2]<2240*mm&&1100*mm<mr[0]&&mr[0]<=2500*mm) {
1364 part = 0; layer = 2; mat = 0;
1365 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1366 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
1367 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
1368 b[2] = tb[2];
1369 ifend = true;
1370 }
1371 if(2240*mm<=mr[2]&&mr[2]<2280*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
1372 part = 0; layer = 2; mat = 1;
1373 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1374 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
1375 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
1376 b[2] = tb[2];
1377 ifend = true;
1378 }
1379 if(2280*mm<=mr[2]&&mr[2]<2310*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
1380 part = 0; layer = 3; mat = 0;
1381 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1382 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
1383 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
1384 b[2] = tb[2];
1385 ifend = true;
1386 }
1387 if(2310*mm<=mr[2]&&mr[2]<2350*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
1388 part = 0; layer = 3; mat = 1;
1389 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1390 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
1391 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
1392 b[2] = tb[2];
1393 ifend = true;
1394 }
1395 if(2350*mm<=mr[2]&&mr[2]<2380*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
1396 part = 0; layer = 4; mat = 0;
1397 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1398 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
1399 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
1400 b[2] = tb[2];
1401 ifend = true;
1402 }
1403 if(2380*mm<=mr[2]&&mr[2]<2420*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
1404 part = 0; layer = 4; mat = 1;
1405 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1406 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
1407 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
1408 b[2] = tb[2];
1409 ifend = true;
1410 }
1411 if(2420*mm<=mr[2]&&mr[2]<2470*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
1412 part = 0; layer = 5; mat = 0;
1413 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1414 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
1415 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
1416 b[2] = tb[2];
1417 ifend = true;
1418 }
1419 if(2470*mm<=mr[2]&&mr[2]<2510*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
1420 part = 0; layer = 5; mat =1;
1421 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1422 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
1423 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
1424 b[2] = tb[2];
1425 ifend = true;
1426 }
1427 if(2510*mm<=mr[2]&&mr[2]<2590*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
1428 part = 0; layer = 6; mat = 0;
1429 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1430 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
1431 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
1432 b[2] = tb[2];
1433 ifend = true;
1434 }
1435 if(2590*mm<=mr[2]&&mr[2]<2630*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
1436 part = 0; layer = 6; mat = 1;
1437 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1438 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
1439 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
1440 b[2] = tb[2];
1441 ifend = true;
1442 }
1443 if(2630*mm<=mr[2]&&mr[2]<2710*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
1444 part = 0; layer = 7; mat = 0;
1445 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1446 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
1447 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
1448 b[2] = tb[2];
1449 ifend = true;
1450 }
1451 if(2710*mm<=mr[2]&&mr[2]<2750*mm&&1362*mm<=mr[0]&&mr[0]<=2500*mm) {
1452 part = 0; layer = 7; mat = 1;
1453 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1454 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
1455 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
1456 b[2] = tb[2];
1457 ifend = true;
1458 }
1459 if(2750*mm<=mr[2]&&mr[2]<=2800*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
1460 part = 0; layer = 8; mat = 0;
1461 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1462 b[0] = tb[0]*cos(pi/4)-tb[1]*sin(pi/4);
1463 b[1] = tb[0]*sin(pi/4)+tb[1]*cos(pi/4);
1464 b[2] = tb[2];
1465 ifend = true;
1466 }
1467 }
1468 if(pi/4<=theta&&theta<3*pi/8) {
1469 mr[0] = fabs(r.x())*cos(pi/4)+fabs(r.y())*sin(pi/4);
1470 mr[1] = fabs(r.x())*sin(pi/4)-fabs(r.y())*cos(pi/4);
1471 mr[2] = fabs(r.z());
1472 if(mr[2]<=1970*mm&&1740*mm<=mr[0]&&mr[0]<=2620*mm) {
1473 part = 1;
1474 if(1740*mm<=mr[0]&&mr[0]<1770*mm) { layer = 0; mat = 0; }
1475 if(1770*mm<=mr[0]&&mr[0]<1810*mm) { layer = 0; mat = 1; }
1476 if(1810*mm<=mr[0]&&mr[0]<1840*mm) { layer = 1; mat = 0; }
1477 if(1840*mm<=mr[0]&&mr[0]<1880*mm) { layer = 1; mat = 1; }
1478 if(1880*mm<=mr[0]&&mr[0]<1910*mm) { layer = 2; mat = 0; }
1479 if(1910*mm<=mr[0]&&mr[0]<1950*mm) { layer = 2; mat = 1; }
1480 if(1950*mm<=mr[0]&&mr[0]<1990*mm) { layer = 3; mat = 0; }
1481 if(1990*mm<=mr[0]&&mr[0]<2030*mm) { layer = 3; mat = 1; }
1482 if(2030*mm<=mr[0]&&mr[0]<2070*mm) { layer = 4; mat = 0; }
1483 if(2070*mm<=mr[0]&&mr[0]<2110*mm) { layer = 4; mat = 1; }
1484 if(2110*mm<=mr[0]&&mr[0]<2190*mm) { layer = 5; mat = 0; }
1485 if(2190*mm<=mr[0]&&mr[0]<2230*mm) { layer = 5; mat = 1; }
1486 if(2230*mm<=mr[0]&&mr[0]<2310*mm) { layer = 6; mat = 0; }
1487 if(2310*mm<=mr[0]&&mr[0]<2350*mm) { layer = 6; mat = 1; }
1488 if(2350*mm<=mr[0]&&mr[0]<2430*mm) { layer = 7; mat = 0; }
1489 if(2430*mm<=mr[0]&&mr[0]<2470*mm) { layer = 7; mat = 1; }
1490 if(2470*mm<=mr[0]&&mr[0]<=2620*mm) { layer = 8; mat = 0; }
1491 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1492 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
1493 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
1494 b[2] = tb[2];
1495 ifbar = true;
1496 }
1497 if(2050*mm<=mr[2]&&mr[2]<2090*mm&&1034*mm<=mr[0]&&mr[0]<=2500*mm){
1498 part = 0; layer = 0; mat = 0;
1499 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1500 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
1501 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
1502 b[2] = tb[2];
1503 ifend = true;
1504 }
1505 if(2090*mm<=mr[2]&&mr[2]<2130*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
1506 part = 0; layer = 0; mat = 1;
1507 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1508 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
1509 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
1510 b[2] = tb[2];
1511 ifend = true;
1512 }
1513 if(2130*mm<=mr[2]&&mr[2]<2170*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
1514 part = 0; layer = 1; mat = 0;
1515 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1516 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
1517 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
1518 b[2] = tb[2];
1519 ifend = true;
1520 }
1521 if(2170*mm<=mr[2]&&mr[2]<2210*mm&&1100*mm<=mr[0]&&mr[0]<=2500*mm) {
1522 part = 0; layer = 1; mat = 1;
1523 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1524 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
1525 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
1526 b[2] = tb[2];
1527 ifend = true;
1528 }
1529 if(2210*mm<=mr[2]&&mr[2]<2240*mm&&1100*mm<mr[0]&&mr[0]<=2500*mm) {
1530 part = 0; layer = 2; mat = 0;
1531 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1532 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
1533 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
1534 b[2] = tb[2];
1535 ifend = true;
1536 }
1537 if(2240*mm<=mr[2]&&mr[2]<2280*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
1538 part = 0; layer = 2; mat = 1;
1539 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1540 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
1541 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
1542 b[2] = tb[2];
1543 ifend = true;
1544 }
1545 if(2280*mm<=mr[2]&&mr[2]<2310*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
1546 part = 0; layer = 3; mat = 0;
1547 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1548 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
1549 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
1550 b[2] = tb[2];
1551 ifend = true;
1552 }
1553 if(2310*mm<=mr[2]&&mr[2]<2350*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
1554 part = 0; layer = 3; mat = 1;
1555 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1556 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
1557 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
1558 b[2] = tb[2];
1559 ifend = true;
1560 }
1561 if(2350*mm<=mr[2]&&mr[2]<2380*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
1562 part = 0; layer = 4; mat = 0;
1563 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1564 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
1565 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
1566 b[2] = tb[2];
1567 ifend = true;
1568 }
1569 if(2380*mm<=mr[2]&&mr[2]<2420*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
1570 part = 0; layer = 4; mat = 1;
1571 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1572 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
1573 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
1574 b[2] = tb[2];
1575 ifend = true;
1576 }
1577 if(2420*mm<=mr[2]&&mr[2]<2470*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
1578 part = 0; layer = 5; mat = 0;
1579 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1580 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
1581 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
1582 b[2] = tb[2];
1583 ifend = true;
1584 }
1585 if(2470*mm<=mr[2]&&mr[2]<2510*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
1586 part = 0; layer = 5; mat =1;
1587 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1588 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
1589 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
1590 b[2] = tb[2];
1591 ifend = true;
1592 }
1593 if(2510*mm<=mr[2]&&mr[2]<2590*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
1594 part = 0; layer = 6; mat = 0;
1595 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1596 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
1597 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
1598 b[2] = tb[2];
1599 ifend = true;
1600 }
1601 if(2590*mm<=mr[2]&&mr[2]<2630*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
1602 part = 0; layer = 6; mat = 1;
1603 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1604 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
1605 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
1606 b[2] = tb[2];
1607 ifend = true;
1608 }
1609 if(2630*mm<=mr[2]&&mr[2]<2710*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
1610 part = 0; layer = 7; mat = 0;
1611 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1612 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
1613 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
1614 b[2] = tb[2];
1615 ifend = true;
1616 }
1617 if(2710*mm<=mr[2]&&mr[2]<2750*mm&&1362*mm<=mr[0]&&mr[0]<=2500*mm) {
1618 part = 0; layer = 7; mat = 1;
1619 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1620 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
1621 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
1622 b[2] = tb[2];
1623 ifend = true;
1624 }
1625 if(2750*mm<=mr[2]&&mr[2]<=2800*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
1626 part = 0; layer = 8; mat = 0;
1627 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1628 b[0] = tb[0]*cos(pi/4)+tb[1]*sin(pi/4);
1629 b[1] = tb[0]*sin(pi/4)-tb[1]*cos(pi/4);
1630 b[2] = tb[2];
1631 ifend = true;
1632 }
1633 }
1634 if(3*pi/8<=theta&&theta<pi/2) {
1635 mr[0] = fabs(r.y()); mr[1] = -fabs(r.x()); mr[2] = fabs(r.z());
1636 if(mr[2]<=1970*mm&&1740*mm<=mr[0]&&mr[0]<=2620*mm) {
1637 part = 1;
1638 if(1740*mm<=mr[0]&&mr[0]<1770*mm) { layer = 0; mat = 0; }
1639 if(1770*mm<=mr[0]&&mr[0]<1810*mm) { layer = 0; mat = 1; }
1640 if(1810*mm<=mr[0]&&mr[0]<1840*mm) { layer = 1; mat = 0; }
1641 if(1840*mm<=mr[0]&&mr[0]<1880*mm) { layer = 1; mat = 1; }
1642 if(1880*mm<=mr[0]&&mr[0]<1910*mm) { layer = 2; mat = 0; }
1643 if(1910*mm<=mr[0]&&mr[0]<1950*mm) { layer = 2; mat = 1; }
1644 if(1950*mm<=mr[0]&&mr[0]<1990*mm) { layer = 3; mat = 0; }
1645 if(1990*mm<=mr[0]&&mr[0]<2030*mm) { layer = 3; mat = 1; }
1646 if(2030*mm<=mr[0]&&mr[0]<2070*mm) { layer = 4; mat = 0; }
1647 if(2070*mm<=mr[0]&&mr[0]<2110*mm) { layer = 4; mat = 1; }
1648 if(2110*mm<=mr[0]&&mr[0]<2190*mm) { layer = 5; mat = 0; }
1649 if(2190*mm<=mr[0]&&mr[0]<2230*mm) { layer = 5; mat = 1; }
1650 if(2230*mm<=mr[0]&&mr[0]<2310*mm) { layer = 6; mat = 0; }
1651 if(2310*mm<=mr[0]&&mr[0]<2350*mm) { layer = 6; mat = 1; }
1652 if(2350*mm<=mr[0]&&mr[0]<2430*mm) { layer = 7; mat = 0; }
1653 if(2430*mm<=mr[0]&&mr[0]<2470*mm) { layer = 7; mat = 1; }
1654 if(2470*mm<=mr[0]&&mr[0]<=2620*mm) { layer = 8; mat = 0; }
1655 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1656 b[0] = -tb[1];
1657 b[1] = tb[0];
1658 b[2] = tb[2];
1659 ifbar = true;
1660 }
1661 if(2050*mm<=mr[2]&&mr[2]<2090*mm&&1034*mm<=mr[0]&&mr[0]<=2500*mm){
1662 part = 0; layer = 0; mat = 0;
1663 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1664 b[0] = -tb[1];
1665 b[1] = tb[0];
1666 b[2] = tb[2];
1667 ifend = true;
1668 }
1669 if(2090*mm<=mr[2]&&mr[2]<2130*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
1670 part = 0; layer = 0; mat = 1;
1671 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1672 b[0] = -tb[1];
1673 b[1] = tb[0];
1674 b[2] = tb[2];
1675 ifend = true;
1676 }
1677 if(2130*mm<=mr[2]&&mr[2]<2170*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
1678 part = 0; layer = 1; mat = 0;
1679 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1680 b[0] = -tb[1];
1681 b[1] = tb[0];
1682 b[2] = tb[2];
1683 ifend = true;
1684 }
1685 if(2170*mm<=mr[2]&&mr[2]<2210*mm&&1100*mm<=mr[0]&&mr[0]<=2500*mm) {
1686 part = 0; layer = 1; mat = 1;
1687 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1688 b[0] = -tb[1];
1689 b[1] = tb[0];
1690 b[2] = tb[2];
1691 ifend = true;
1692 }
1693 if(2210*mm<=mr[2]&&mr[2]<2240*mm&&1100*mm<mr[0]&&mr[0]<=2500*mm) {
1694 part = 0; layer = 2; mat = 0;
1695 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1696 b[0] = -tb[1];
1697 b[1] = tb[0];
1698 b[2] = tb[2];
1699 ifend = true;
1700 }
1701 if(2240*mm<=mr[2]&&mr[2]<2280*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
1702 part = 0; layer = 2; mat = 1;
1703 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1704 b[0] = -tb[1];
1705 b[1] = tb[0];
1706 b[2] = tb[2];
1707 ifend = true;
1708 }
1709 if(2280*mm<=mr[2]&&mr[2]<2310*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
1710 part = 0; layer = 3; mat = 0;
1711 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1712 b[0] = -tb[1];
1713 b[1] = tb[0];
1714 b[2] = tb[2];
1715 ifend = true;
1716 }
1717 if(2310*mm<=mr[2]&&mr[2]<2350*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
1718 part = 0; layer = 3; mat = 1;
1719 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1720 b[0] = -tb[1];
1721 b[1] = tb[0];
1722 b[2] = tb[2];
1723 ifend = true;
1724 }
1725 if(2350*mm<=mr[2]&&mr[2]<2380*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
1726 part = 0; layer = 4; mat = 0;
1727 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1728 b[0] = -tb[1];
1729 b[1] = tb[0];
1730 b[2] = tb[2];
1731 ifend = true;
1732 }
1733 if(2380*mm<=mr[2]&&mr[2]<2420*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
1734 part = 0; layer = 4; mat = 1;
1735 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1736 b[0] = -tb[1];
1737 b[1] = tb[0];
1738 b[2] = tb[2];
1739 ifend = true;
1740 }
1741 if(2420*mm<=mr[2]&&mr[2]<2470*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
1742 part = 0; layer = 5; mat = 0;
1743 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1744 b[0] = -tb[1];
1745 b[1] = tb[0];
1746 b[2] = tb[2];
1747 ifend = true;
1748 }
1749 if(2470*mm<=mr[2]&&mr[2]<2510*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
1750 part = 0; layer = 5; mat =1;
1751 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1752 b[0] = -tb[1];
1753 b[1] = tb[0];
1754 b[2] = tb[2];
1755 ifend = true;
1756 }
1757 if(2510*mm<=mr[2]&&mr[2]<2590*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
1758 part = 0; layer = 6; mat = 0;
1759 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1760 b[0] = -tb[1];
1761 b[1] = tb[0];
1762 b[2] = tb[2];
1763 ifend = true;
1764 }
1765 if(2590*mm<=mr[2]&&mr[2]<2630*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
1766 part = 0; layer = 6; mat = 1;
1767 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1768 b[0] = -tb[1];
1769 b[1] = tb[0];
1770 b[2] = tb[2];
1771 ifend = true;
1772 }
1773 if(2630*mm<=mr[2]&&mr[2]<2710*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
1774 part = 0; layer = 7; mat = 0;
1775 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1776 b[0] = -tb[1];
1777 b[1] = tb[0];
1778 b[2] = tb[2];
1779 ifend = true;
1780 }
1781 if(2710*mm<=mr[2]&&mr[2]<2750*mm&&1362*mm<=mr[0]&&mr[0]<=2500*mm) {
1782 part = 0; layer = 7; mat = 1;
1783 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1784 b[0] = -tb[1];
1785 b[1] = tb[0];
1786 b[2] = tb[2];
1787 ifend = true;
1788 }
1789 if(2750*mm<=mr[2]&&mr[2]<=2800*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
1790 part = 0; layer = 8; mat = 0;
1791 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1792 b[0] = -tb[1];
1793 b[1] = tb[0];
1794 b[2] = tb[2];
1795 ifend = true;
1796 }
1797 }
1798 }
1799 if(r.x()==0.) {
1800 mr[0] = fabs(r.y()); mr[1] = -fabs(r.x()); mr[2] = fabs(r.z());
1801 if(mr[2]<=1970*mm&&1740*mm<=mr[0]&&mr[0]<=2620*mm) {
1802 part = 1;
1803 if(1740*mm<=mr[0]&&mr[0]<1770*mm) { layer = 0; mat = 0; }
1804 if(1770*mm<=mr[0]&&mr[0]<1810*mm) { layer = 0; mat = 1; }
1805 if(1810*mm<=mr[0]&&mr[0]<1840*mm) { layer = 1; mat = 0; }
1806 if(1840*mm<=mr[0]&&mr[0]<1880*mm) { layer = 1; mat = 1; }
1807 if(1880*mm<=mr[0]&&mr[0]<1910*mm) { layer = 2; mat = 0; }
1808 if(1910*mm<=mr[0]&&mr[0]<1950*mm) { layer = 2; mat = 1; }
1809 if(1950*mm<=mr[0]&&mr[0]<1990*mm) { layer = 3; mat = 0; }
1810 if(1990*mm<=mr[0]&&mr[0]<2030*mm) { layer = 3; mat = 1; }
1811 if(2030*mm<=mr[0]&&mr[0]<2070*mm) { layer = 4; mat = 0; }
1812 if(2070*mm<=mr[0]&&mr[0]<2110*mm) { layer = 4; mat = 1; }
1813 if(2110*mm<=mr[0]&&mr[0]<2190*mm) { layer = 5; mat = 0; }
1814 if(2190*mm<=mr[0]&&mr[0]<2230*mm) { layer = 5; mat = 1; }
1815 if(2230*mm<=mr[0]&&mr[0]<2310*mm) { layer = 6; mat = 0; }
1816 if(2310*mm<=mr[0]&&mr[0]<2350*mm) { layer = 6; mat = 1; }
1817 if(2350*mm<=mr[0]&&mr[0]<2430*mm) { layer = 7; mat = 0; }
1818 if(2430*mm<=mr[0]&&mr[0]<2470*mm) { layer = 7; mat = 1; }
1819 if(2470*mm<=mr[0]&&mr[0]<=2620*mm) { layer = 8; mat = 0; }
1820 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1821 b[0] = -tb[1];
1822 b[1] = tb[0];
1823 b[2] = tb[2];
1824 ifbar = true;
1825 }
1826 if(2050*mm<=mr[2]&&mr[2]<2090*mm&&1034*mm<=mr[0]&&mr[0]<=2500*mm){
1827 part = 0; layer = 0; mat = 0;
1828 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1829 b[0] = -tb[1];
1830 b[1] = tb[0];
1831 b[2] = tb[2];
1832 ifend = true;
1833 }
1834 if(2090*mm<=mr[2]&&mr[2]<2130*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
1835 part = 0; layer = 0; mat = 1;
1836 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1837 b[0] = -tb[1];
1838 b[1] = tb[0];
1839 b[2] = tb[2];
1840 ifend = true;
1841 }
1842 if(2130*mm<=mr[2]&&mr[2]<2170*mm&&1067*mm<=mr[0]&&mr[0]<=2500*mm) {
1843 part = 0; layer = 1; mat = 0;
1844 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1845 b[0] = -tb[1];
1846 b[1] = tb[0];
1847 b[2] = tb[2];
1848 ifend = true;
1849 }
1850 if(2170*mm<=mr[2]&&mr[2]<2210*mm&&1100*mm<=mr[0]&&mr[0]<=2500*mm) {
1851 part = 0; layer = 1; mat = 1;
1852 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1853 b[0] = -tb[1];
1854 b[1] = tb[0];
1855 b[2] = tb[2];
1856 ifend = true;
1857 }
1858 if(2210*mm<=mr[2]&&mr[2]<2240*mm&&1100*mm<mr[0]&&mr[0]<=2500*mm) {
1859 part = 0; layer = 2; mat = 0;
1860 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1861 b[0] = -tb[1];
1862 b[1] = tb[0];
1863 b[2] = tb[2];
1864 ifend = true;
1865 }
1866 if(2240*mm<=mr[2]&&mr[2]<2280*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
1867 part = 0; layer = 2; mat = 1;
1868 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1869 b[0] = -tb[1];
1870 b[1] = tb[0];
1871 b[2] = tb[2];
1872 ifend = true;
1873 }
1874 if(2280*mm<=mr[2]&&mr[2]<2310*mm&&1133*mm<=mr[0]&&mr[0]<=2500*mm) {
1875 part = 0; layer = 3; mat = 0;
1876 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1877 b[0] = -tb[1];
1878 b[1] = tb[0];
1879 b[2] = tb[2];
1880 ifend = true;
1881 }
1882 if(2310*mm<=mr[2]&&mr[2]<2350*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
1883 part = 0; layer = 3; mat = 1;
1884 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1885 b[0] = -tb[1];
1886 b[1] = tb[0];
1887 b[2] = tb[2];
1888 ifend = true;
1889 }
1890 if(2350*mm<=mr[2]&&mr[2]<2380*mm&&1167*mm<=mr[0]&&mr[0]<=2500*mm) {
1891 part = 0; layer = 4; mat = 0;
1892 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1893 b[0] = -tb[1];
1894 b[1] = tb[0];
1895 b[2] = tb[2];
1896 ifend = true;
1897 }
1898 if(2380*mm<=mr[2]&&mr[2]<2420*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
1899 part = 0; layer = 4; mat = 1;
1900 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1901 b[0] = -tb[1];
1902 b[1] = tb[0];
1903 b[2] = tb[2];
1904 ifend = true;
1905 }
1906 if(2420*mm<=mr[2]&&mr[2]<2470*mm&&1203*mm<=mr[0]&&mr[0]<=2500*mm) {
1907 part = 0; layer = 5; mat = 0;
1908 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1909 b[0] = -tb[1];
1910 b[1] = tb[0];
1911 b[2] = tb[2];
1912 ifend = true;
1913 }
1914 if(2470*mm<=mr[2]&&mr[2]<2510*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
1915 part = 0; layer = 5; mat =1;
1916 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1917 b[0] = -tb[1];
1918 b[1] = tb[0];
1919 b[2] = tb[2];
1920 ifend = true;
1921 }
1922 if(2510*mm<=mr[2]&&mr[2]<2590*mm&&1241*mm<=mr[0]&&mr[0]<=2500*mm) {
1923 part = 0; layer = 6; mat = 0;
1924 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1925 b[0] = -tb[1];
1926 b[1] = tb[0];
1927 b[2] = tb[2];
1928 ifend = true;
1929 }
1930 if(2590*mm<=mr[2]&&mr[2]<2630*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
1931 part = 0; layer = 6; mat = 1;
1932 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1933 b[0] = -tb[1];
1934 b[1] = tb[0];
1935 b[2] = tb[2];
1936 ifend = true;
1937 }
1938 if(2630*mm<=mr[2]&&mr[2]<2710*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
1939 part = 0; layer = 7; mat = 0;
1940 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1941 b[0] = -tb[1];
1942 b[1] = tb[0];
1943 b[2] = tb[2];
1944 ifend = true;
1945 }
1946 if(2710*mm<=mr[2]&&mr[2]<2750*mm&&1362*mm<=mr[0]&&mr[0]<=2500*mm) {
1947 part = 0; layer = 7; mat = 1;
1948 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1949 b[0] = -tb[1];
1950 b[1] = tb[0];
1951 b[2] = tb[2];
1952 ifend = true;
1953 }
1954 if(2750*mm<=mr[2]&&mr[2]<=2800*mm&&1302*mm<=mr[0]&&mr[0]<=2500*mm) {
1955 part = 0; layer = 8; mat = 0;
1956 m_Mucfield->getMucField(part,layer,mat,mr,tb);
1957 b[0] = -tb[1];
1958 b[1] = tb[0];
1959 b[2] = tb[2];
1960 ifend = true;
1961 }
1962 }
1963 if(ifbar==true||ifend==true) {
1964 if( r.x() < 0. && r.y() >= 0. && r.z() > 0. ){
1965 b[0] = -b[0];
1966 }
1967 else if( r.x() <= 0. && r.y() < 0. && r.z() > 0. ){
1968 b[0] = -b[0];
1969 b[1] = -b[1];
1970 }
1971 else if( r.x() > 0. && r.y() < 0. && r.z() > 0. ){
1972 b[1] = -b[1];
1973 }
1974 else if( r.x() >= 0. && r.y() > 0. && r.z() < 0. ){
1975 b[0] = -b[0];
1976 b[1] = -b[1];
1977 }
1978 else if( r.x() < 0. && r.y() >= 0. && r.z() < 0. ){
1979 b[1] = -b[1];
1980 }
1981 else if( r.x() <= 0. && r.y() < 0. && r.z() < 0. ){
1982 b[0] = b[0];
1983 b[1] = b[1];
1984 }
1985 else if( r.x() > 0. && r.y() <= 0. && r.z() < 0. ){
1986 b[0] = -b[0];
1987 }
1988 }
1989 }
1990
1991 //reference frames defined by simulation and SCM are different.
1992 newb[0] = -b[0] * m_scale;
1993 newb[1] = b[1] * m_scale;
1994 newb[2] = -b[2] * m_scale;
1995/* 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) {
1996 return StatusCode::SUCCESS;
1997 }
1998 else {
1999 newb[0] = newb[0] - 0.10*tesla;
2000 newb[1] = newb[1] + 0.10*tesla;
2001 newb[2] = newb[2] - 0.10*tesla;
2002 }*/
2003
2004 //newb[0] = b[0];
2005 //newb[1] = b[1];
2006 //newb[2] = b[2];
2007
2008#ifndef BEAN
2009 return StatusCode::SUCCESS;
2010#else
2011 return true;
2012#endif
2013}
2014
2015StatusCode MagneticFieldSvc::uniFieldVector(const HepPoint3D& r,
2016 HepVector3D& b) const
2017{
2018 if(m_runmode == 8 || m_runmode == 7) {
2019 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)
2020 {
2021 b[0]=0.0;
2022 b[1]=0.0;
2023 b[2]=-0.9*tesla;
2024 }
2025 else
2026 {
2027 b[0]=0.0;
2028 b[1]=0.0;
2029 b[2]=0.0;
2030 }
2031 }
2032 else {
2033 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)
2034 {
2035 b[0]=0.0;
2036 b[1]=0.0;
2037 b[2]=-1.0*tesla;
2038 }
2039 else
2040 {
2041 b[0]=0.0;
2042 b[1]=0.0;
2043 b[2]=0.0;
2044 }
2045 }
2046
2047 if(m_turnOffField == true) {
2048 b[0] = 0.;
2049 b[1] = 0.;
2050 b[2] = 0.;
2051 }
2052 //yzhang add 2012-04-25
2053 b[0] *= m_scale;
2054 b[1] *= m_scale;
2055 b[2] *= m_scale;
2056#ifndef BEAN
2057 return StatusCode::SUCCESS;
2058#else
2059 return true;
2060#endif
2061}
2062
2063double MagneticFieldSvc::getReferField()
2064{
2065 if(m_useDB == false) {
2066 if(m_runmode == 8 || m_runmode == 7) m_zfield = -0.9*tesla;
2067 else m_zfield = -1.0*tesla;
2068 }
2069
2070 if(m_turnOffField == true) {
2071 m_zfield = 0.;
2072 }
2073 return m_zfield * m_scale;
2074}
2075
2076bool MagneticFieldSvc::ifRealField() const {
2077 return m_ifRealField;
2078}
2079
2080//=============================================================================
2081// routine to fill the field vector
2082//=============================================================================
2083void MagneticFieldSvc::fieldGrid (const HepPoint3D& r,
2084 HepVector3D& bf ) const {
2085
2086 bf[0] = 0.0;
2087 bf[1] = 0.0;
2088 bf[2] = 0.0;
2089
2090 /// Linear interpolated field
2091 double z = r.z() - m_zOffSet;
2092 if( z < m_min_FL[2] || z > m_max_FL[2] ) return;
2093 double x = r.x();
2094 if( x < m_min_FL[0] || x > m_max_FL[0] ) return;
2095 double y = r.y();
2096 if( y < m_min_FL[1] || y > m_max_FL[1] ) return;
2097 int i = int( (x-m_min_FL[0])/m_Dxyz[0]);
2098 int j = int( (y-m_min_FL[1])/m_Dxyz[1] );
2099 int k = int( (z-m_min_FL[2])/m_Dxyz[2] );
2100
2101 int ijk000 = 3*( m_Nxyz[0]*( m_Nxyz[1]*k + j ) + i );
2102 int ijk001 = 3*( m_Nxyz[0]*( m_Nxyz[1]*(k+1) + j ) + i );
2103 int ijk010 = 3*( m_Nxyz[0]*( m_Nxyz[1]*k + j + 1 ) + i );
2104 int ijk011 = 3*( m_Nxyz[0]*( m_Nxyz[1]*(k+1) + j + 1) + i );
2105 int ijk100 = 3*( m_Nxyz[0]*( m_Nxyz[1]*k + j) + i + 1 );
2106 int ijk101 = 3*( m_Nxyz[0]*( m_Nxyz[1]*(k+1) + j) + i + 1 );
2107 int ijk110 = 3*( m_Nxyz[0]*( m_Nxyz[1]*k + j + 1) + i + 1 );
2108 int ijk111 = 3*( m_Nxyz[0]*( m_Nxyz[1]*(k+1) + j + 1 ) + i + 1 );
2109
2110 // auxiliary variables defined at the vertices of the cube that
2111 // contains the (x, y, z) point where the field is interpolated
2112/* std::cout<<"x,y,z: "<<x<<","<<y<<","<<z<<std::endl;
2113 std::cout<<"point1(x,y,z): "<<m_P[ijk000]<<","<<m_P[ijk000+1]<<","<<m_P[ijk000+2]<<std::endl;
2114 std::cout<<"point2(x,y,z): "<<m_P[ijk001]<<","<<m_P[ijk001+1]<<","<<m_P[ijk001+2]<<std::endl;
2115 std::cout<<"point3(x,y,z): "<<m_P[ijk010]<<","<<m_P[ijk010+1]<<","<<m_P[ijk010+2]<<std::endl;
2116 std::cout<<"point4(x,y,z): "<<m_P[ijk011]<<","<<m_P[ijk011+1]<<","<<m_P[ijk011+2]<<std::endl;
2117 std::cout<<"point5(x,y,z): "<<m_P[ijk100]<<","<<m_P[ijk100+1]<<","<<m_P[ijk100+2]<<std::endl;
2118 std::cout<<"point6(x,y,z): "<<m_P[ijk101]<<","<<m_P[ijk101+1]<<","<<m_P[ijk101+2]<<std::endl;
2119 std::cout<<"point7(x,y,z): "<<m_P[ijk110]<<","<<m_P[ijk110+1]<<","<<m_P[ijk110+2]<<std::endl;
2120 std::cout<<"point8(x,y,z): "<<m_P[ijk111]<<","<<m_P[ijk111+1]<<","<<m_P[ijk111+2]<<std::endl;
2121
2122 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])
2123 std::cout<<"FATALERRORX****************************"<<std::endl;
2124 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])
2125 std::cout<<"FATALERRORY***************************"<<std::endl;
2126 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])
2127 std::cout<<"FATALERRORZ****************************"<<std::endl; */
2128 double cx000 = m_Q[ ijk000 ];
2129 double cx001 = m_Q[ ijk001 ];
2130 double cx010 = m_Q[ ijk010 ];
2131 double cx011 = m_Q[ ijk011 ];
2132 double cx100 = m_Q[ ijk100 ];
2133 double cx101 = m_Q[ ijk101 ];
2134 double cx110 = m_Q[ ijk110 ];
2135 double cx111 = m_Q[ ijk111 ];
2136 double cy000 = m_Q[ ijk000+1 ];
2137 double cy001 = m_Q[ ijk001+1 ];
2138 double cy010 = m_Q[ ijk010+1 ];
2139 double cy011 = m_Q[ ijk011+1 ];
2140 double cy100 = m_Q[ ijk100+1 ];
2141 double cy101 = m_Q[ ijk101+1 ];
2142 double cy110 = m_Q[ ijk110+1 ];
2143 double cy111 = m_Q[ ijk111+1 ];
2144 double cz000 = m_Q[ ijk000+2 ];
2145 double cz001 = m_Q[ ijk001+2 ];
2146 double cz010 = m_Q[ ijk010+2 ];
2147 double cz011 = m_Q[ ijk011+2 ];
2148 double cz100 = m_Q[ ijk100+2 ];
2149 double cz101 = m_Q[ ijk101+2 ];
2150 double cz110 = m_Q[ ijk110+2 ];
2151 double cz111 = m_Q[ ijk111+2 ];
2152 double hx1 = ( x-m_min_FL[0]-i*m_Dxyz[0] )/m_Dxyz[0];
2153 double hy1 = ( y-m_min_FL[1]-j*m_Dxyz[1] )/m_Dxyz[1];
2154 double hz1 = ( z-m_min_FL[2]-k*m_Dxyz[2] )/m_Dxyz[2];
2155 double hx0 = 1.0-hx1;
2156 double hy0 = 1.0-hy1;
2157 double hz0 = 1.0-hz1;
2158 double h000 = hx0*hy0*hz0;
2159 if( fabs(h000) > 0.0 &&
2160 cx000 > 9.0e5 && cy000 > 9.0e5 && cz000 > 9.0e5) return;
2161
2162 double h001 = hx0*hy0*hz1;
2163 if( fabs(h001) > 0.0 &&
2164 cx001 > 9.0e5 && cy001 > 9.0e5 && cz001 > 9.0e5) return;
2165
2166 double h010 = hx0*hy1*hz0;
2167 if( fabs(h010) > 0.0 &&
2168 cx010 > 9.0e5 && cy010 > 9.0e5 && cz010 > 9.0e5) return;
2169
2170 double h011 = hx0*hy1*hz1;
2171 if( fabs(h011) > 0.0 &&
2172 cx011 > 9.0e5 && cy011 > 9.0e5 && cz011 > 9.0e5) return;
2173
2174 double h100 = hx1*hy0*hz0;
2175 if( fabs(h100) > 0.0 &&
2176 cx100 > 9.0e5 && cy100 > 9.0e5 && cz100 > 9.0e5) return;
2177
2178 double h101 = hx1*hy0*hz1;
2179 if( fabs(h101) > 0.0 &&
2180 cx101 > 9.0e5 && cy101 > 9.0e5 && cz101 > 9.0e5) return;
2181
2182 double h110 = hx1*hy1*hz0;
2183 if( fabs(h110) > 0.0 &&
2184 cx110 > 9.0e5 && cy110 > 9.0e5 && cz110 > 9.0e5) return;
2185
2186 double h111 = hx1*hy1*hz1;
2187 if( fabs(h111) > 0.0 &&
2188 cx111 > 9.0e5 && cy111 > 9.0e5 && cz111 > 9.0e5) return;
2189
2190 bf(0) = ( cx000*h000 + cx001*h001 + cx010*h010 + cx011*h011 +
2191 cx100*h100 + cx101*h101 + cx110*h110 + cx111*h111);
2192 bf(1) = ( cy000*h000 + cy001*h001 + cy010*h010 + cy011*h011 +
2193 cy100*h100 + cy101*h101 + cy110*h110 + cy111*h111 );
2194 bf(2) = ( cz000*h000 + cz001*h001 + cz010*h010 + cz011*h011 +
2195 cz100*h100 + cz101*h101 + cz110*h110 + cz111*h111 );
2196 return;
2197}
2198
2199//=============================================================================
2200// routine to fill the field vector
2201//=============================================================================
2202void MagneticFieldSvc::fieldGrid_TE (const HepPoint3D& r,
2203 HepVector3D& bf ) const {
2204
2205 bf[0] = 0.0;
2206 bf[1] = 0.0;
2207 bf[2] = 0.0;
2208
2209 /// Linear interpolated field
2210 double z = r.z() - m_zOffSet_TE;
2211 if( fabs(z) < m_min_FL_TE[2] || fabs(z) > m_max_FL_TE[2] ) return;
2212 double x = r.x();
2213 if( fabs(x) < m_min_FL_TE[0] || fabs(x) > m_max_FL_TE[0] ) return;
2214 double y = r.y();
2215 if( fabs(y) < m_min_FL_TE[1] || fabs(y) > m_max_FL_TE[1] ) return;
2216 int i = int( (fabs(x)-m_min_FL_TE[0])/m_Dxyz_TE[0]);
2217 int j = int( (fabs(y)-m_min_FL_TE[1])/m_Dxyz_TE[1] );
2218 int k = int( (fabs(z)-m_min_FL_TE[2])/m_Dxyz_TE[2] );
2219
2220 int ijk000 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*k + j ) + i );
2221 int ijk001 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*(k+1) + j ) + i );
2222 int ijk010 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*k + j + 1 ) + i );
2223 int ijk011 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*(k+1) + j + 1) + i );
2224 int ijk100 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*k + j) + i + 1 );
2225 int ijk101 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*(k+1) + j) + i + 1 );
2226 int ijk110 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*k + j + 1) + i + 1 );
2227 int ijk111 = 3*( m_Nxyz_TE[0]*( m_Nxyz_TE[1]*(k+1) + j + 1 ) + i + 1 );
2228
2229 // auxiliary variables defined at the vertices of the cube that
2230 // contains the (x, y, z) point where the field is interpolated
2231/* std::cout<<"x,y,z: "<<x<<","<<y<<","<<z<<std::endl;
2232 std::cout<<"point1(x,y,z): "<<m_P[ijk000]<<","<<m_P[ijk000+1]<<","<<m_P[ijk000+2]<<std::endl;
2233 std::cout<<"point2(x,y,z): "<<m_P[ijk001]<<","<<m_P[ijk001+1]<<","<<m_P[ijk001+2]<<std::endl;
2234 std::cout<<"point3(x,y,z): "<<m_P[ijk010]<<","<<m_P[ijk010+1]<<","<<m_P[ijk010+2]<<std::endl;
2235 std::cout<<"point4(x,y,z): "<<m_P[ijk011]<<","<<m_P[ijk011+1]<<","<<m_P[ijk011+2]<<std::endl;
2236 std::cout<<"point5(x,y,z): "<<m_P[ijk100]<<","<<m_P[ijk100+1]<<","<<m_P[ijk100+2]<<std::endl;
2237 std::cout<<"point6(x,y,z): "<<m_P[ijk101]<<","<<m_P[ijk101+1]<<","<<m_P[ijk101+2]<<std::endl;
2238 std::cout<<"point7(x,y,z): "<<m_P[ijk110]<<","<<m_P[ijk110+1]<<","<<m_P[ijk110+2]<<std::endl;
2239 std::cout<<"point8(x,y,z): "<<m_P[ijk111]<<","<<m_P[ijk111+1]<<","<<m_P[ijk111+2]<<std::endl;
2240 */
2241 double cx000 = m_Q_TE[ ijk000 ];
2242 double cx001 = m_Q_TE[ ijk001 ];
2243 double cx010 = m_Q_TE[ ijk010 ];
2244 double cx011 = m_Q_TE[ ijk011 ];
2245 double cx100 = m_Q_TE[ ijk100 ];
2246 double cx101 = m_Q_TE[ ijk101 ];
2247 double cx110 = m_Q_TE[ ijk110 ];
2248 double cx111 = m_Q_TE[ ijk111 ];
2249 double cy000 = m_Q_TE[ ijk000+1 ];
2250 double cy001 = m_Q_TE[ ijk001+1 ];
2251 double cy010 = m_Q_TE[ ijk010+1 ];
2252 double cy011 = m_Q_TE[ ijk011+1 ];
2253 double cy100 = m_Q_TE[ ijk100+1 ];
2254 double cy101 = m_Q_TE[ ijk101+1 ];
2255 double cy110 = m_Q_TE[ ijk110+1 ];
2256 double cy111 = m_Q_TE[ ijk111+1 ];
2257 double cz000 = m_Q_TE[ ijk000+2 ];
2258 double cz001 = m_Q_TE[ ijk001+2 ];
2259 double cz010 = m_Q_TE[ ijk010+2 ];
2260 double cz011 = m_Q_TE[ ijk011+2 ];
2261 double cz100 = m_Q_TE[ ijk100+2 ];
2262 double cz101 = m_Q_TE[ ijk101+2 ];
2263 double cz110 = m_Q_TE[ ijk110+2 ];
2264 double cz111 = m_Q_TE[ ijk111+2 ];
2265 double hx1 = ( fabs(x)-m_min_FL_TE[0]-i*m_Dxyz_TE[0] )/m_Dxyz_TE[0];
2266 double hy1 = ( fabs(y)-m_min_FL_TE[1]-j*m_Dxyz_TE[1] )/m_Dxyz_TE[1];
2267 double hz1 = ( fabs(z)-m_min_FL_TE[2]-k*m_Dxyz_TE[2] )/m_Dxyz_TE[2];
2268 double hx0 = 1.0-hx1;
2269 double hy0 = 1.0-hy1;
2270 double hz0 = 1.0-hz1;
2271 double h000 = hx0*hy0*hz0;
2272 if( fabs(h000) > 0.0 &&
2273 cx000 > 9.0e5 && cy000 > 9.0e5 && cz000 > 9.0e5) return;
2274
2275 double h001 = hx0*hy0*hz1;
2276 if( fabs(h001) > 0.0 &&
2277 cx001 > 9.0e5 && cy001 > 9.0e5 && cz001 > 9.0e5) return;
2278
2279 double h010 = hx0*hy1*hz0;
2280 if( fabs(h010) > 0.0 &&
2281 cx010 > 9.0e5 && cy010 > 9.0e5 && cz010 > 9.0e5) return;
2282
2283 double h011 = hx0*hy1*hz1;
2284 if( fabs(h011) > 0.0 &&
2285 cx011 > 9.0e5 && cy011 > 9.0e5 && cz011 > 9.0e5) return;
2286
2287 double h100 = hx1*hy0*hz0;
2288 if( fabs(h100) > 0.0 &&
2289 cx100 > 9.0e5 && cy100 > 9.0e5 && cz100 > 9.0e5) return;
2290
2291 double h101 = hx1*hy0*hz1;
2292 if( fabs(h101) > 0.0 &&
2293 cx101 > 9.0e5 && cy101 > 9.0e5 && cz101 > 9.0e5) return;
2294
2295 double h110 = hx1*hy1*hz0;
2296 if( fabs(h110) > 0.0 &&
2297 cx110 > 9.0e5 && cy110 > 9.0e5 && cz110 > 9.0e5) return;
2298
2299 double h111 = hx1*hy1*hz1;
2300 if( fabs(h111) > 0.0 &&
2301 cx111 > 9.0e5 && cy111 > 9.0e5 && cz111 > 9.0e5) return;
2302
2303 bf(0) = ( cx000*h000 + cx001*h001 + cx010*h010 + cx011*h011 +
2304 cx100*h100 + cx101*h101 + cx110*h110 + cx111*h111);
2305 bf(1) = ( cy000*h000 + cy001*h001 + cy010*h010 + cy011*h011 +
2306 cy100*h100 + cy101*h101 + cy110*h110 + cy111*h111 );
2307 bf(2) = ( cz000*h000 + cz001*h001 + cz010*h010 + cz011*h011 +
2308 cz100*h100 + cz101*h101 + cz110*h110 + cz111*h111 );
2309
2310
2311 if( r.x() < 0. && r.y() >= 0. && r.z() > 0. ){
2312 bf(0) = -bf(0);
2313 }
2314 else if( r.x() <= 0. && r.y() < 0. && r.z() > 0. ){
2315 bf(0) = -bf(0);
2316 bf(1) = -bf(1);
2317 }
2318 else if( r.x() > 0. && r.y() < 0. && r.z() > 0. ){
2319 bf(1) = -bf(1);
2320 }
2321
2322 else if( r.x() >= 0. && r.y() > 0. && r.z() < 0. ){
2323 bf(0) = -bf(0);
2324 bf(1) = -bf(1);
2325 }
2326 else if( r.x() < 0. && r.y() >= 0. && r.z() < 0. ){
2327 bf(1) = -bf(1);
2328 }
2329 else if( r.x() <= 0. && r.y() < 0. && r.z() < 0. ){
2330 bf(0) = bf(0);
2331 bf(1) = bf(1);
2332 }
2333 else if( r.x() > 0. && r.y() <= 0. && r.z() < 0. ){
2334 bf(0) = -bf(0);
2335 }
2336
2337 return;
2338}
HepVector3D
HepGeom::Vector3D< double > HepVector3D
Definition Analysis/VertexFit/VertexFit-00-03-07/src/BField.cxx:13
sin
double sin(const BesAngle a)
Definition BesAngle.h:210
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
sprintf
sprintf(cut,"kal_costheta0_em>-0.93&&kal_costheta0_em<0.93&&kal_pxy0_em>=0.05+%d*0.1&&kal_pxy0_em<0.15+%d*0.1&&NGch>=2", j, j)
x
Double_t x[10]
Definition DataBase/tau_mode.c:57
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition Gam4pikp.cxx:37
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition MagneticFieldSvc.h:11
HepVector3D
HepGeom::Vector3D< double > HepVector3D
Definition MagneticFieldSvc.h:12
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:13
NULL
#define NULL
Definition VertexExtrapolate.h:14
MagneticFieldSvc::m_connect_run
FieldDBUtil::ConnectionDB * m_connect_run
Definition MagneticFieldSvc.h:204
MagneticFieldSvc::m_mapMagnetInfo
std::map< int, std::vector< double > > m_mapMagnetInfo
Definition MagneticFieldSvc.h:205
MagneticFieldSvc::getReferField
virtual double getReferField()
Definition MagneticFieldSvc.cxx:2063
MagneticFieldSvc::fieldVector
virtual StatusCode fieldVector(const HepPoint3D &xyz, HepVector3D &fvec) const
Definition MagneticFieldSvc.cxx:1080
MagneticFieldSvc::~MagneticFieldSvc
virtual ~MagneticFieldSvc()
Virtual destructor.
Definition MagneticFieldSvc.cxx:113
MagneticFieldSvc::finalize
virtual StatusCode finalize()
Finalise the service.
Definition MagneticFieldSvc.cxx:695
MagneticFieldSvc::ifRealField
virtual bool ifRealField() const
Definition MagneticFieldSvc.cxx:2076
MagneticFieldSvc::beamEnergy
std::vector< double > beamEnergy
Definition MagneticFieldSvc.h:207
MagneticFieldSvc::initialize
virtual StatusCode initialize()
Initialise the service (Inherited Service overrides)
Definition MagneticFieldSvc.cxx:140
MagneticFieldSvc::m_mapBeamEnergy
std::map< int, std::vector< double > > m_mapBeamEnergy
Definition MagneticFieldSvc.h:206
MagneticFieldSvc::handle
void handle(const Incident &)
Definition MagneticFieldSvc.cxx:727
MagneticFieldSvc::uniFieldVector
virtual StatusCode uniFieldVector(const HepPoint3D &xyz, HepVector3D &fvec) const
Definition MagneticFieldSvc.cxx:2015
MagneticFieldSvc::MagneticFieldSvc
MagneticFieldSvc(const std::string &name, ISvcLocator *svc)
Definition MagneticFieldSvc.cxx:48
MagneticFieldSvc::current
std::vector< double > current
Definition MagneticFieldSvc.h:208
MagneticFieldSvc::m_eventSvc
IDataProviderSvc * m_eventSvc
Definition MagneticFieldSvc.h:201
y
double y[1000]
Definition draw_charge_space_coarse_II.cxx:8
const
EvtCyclic3::c_str
char * c_str(Index i)
Definition EvtCyclic3.cc:252
b
const double b
Definition slope.cxx:9
pi
const float pi
Definition vector3.h:133