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