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