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