BOSS 8.0.0
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/DataObject.h"
4# include "GaudiKernel/IDataProviderSvc.h"
5# include "GaudiKernel/IIncidentSvc.h"
6# include "GaudiKernel/ISvcLocator.h"
7# include "GaudiKernel/Incident.h"
8# include "GaudiKernel/MsgStream.h"
9# include "GaudiKernel/SmartDataPtr.h"
10
11# include "MagneticFieldSvc/IBesMagFieldSvc.h"
12#endif
13
14#include "MagneticFieldSvc.h"
15#include "MucMagneticField.h"
16
17#include "CLHEP/Units/PhysicalConstants.h"
18
19#ifndef BEAN
20# include "EventModel/EventHeader.h"
21# include "EventModel/EventModel.h"
22#endif
23
24#include <cstdlib>
25#include <cstring>
26#include <fstream>
27#include <iostream>
28using namespace std;
29using namespace CLHEP;
30
31/** @file MagneticFieldSvc.cxx
32 * Implementation of MagneticFieldSvc class
33 *
34 */
35
36#ifndef BEAN
37// Instantiation of a static factory class used by clients to create
38// instances of this service
39// static SvcFactory<MagneticFieldSvc> s_factory;
40// const ISvcFactory& MagneticFieldSvcFactory = s_factory;
41
42//=============================================================================
43// Standard constructor, initializes variables
44//=============================================================================
45DECLARE_COMPONENT( MagneticFieldSvc )
46MagneticFieldSvc::MagneticFieldSvc( const std::string& name, ISvcLocator* svc )
47 : base_class( name, svc ) {
48 declareProperty( "TurnOffField", m_turnOffField = false );
49 declareProperty( "FieldMapFile", m_filename );
50 declareProperty( "GridDistance", m_gridDistance = 5 );
51 declareProperty( "RunMode", m_runmode = 2 );
52 declareProperty( "IfRealField", m_ifRealField = true );
53 declareProperty( "OutLevel", m_outlevel = 1 );
54 declareProperty( "Scale", m_scale = 1.0 );
55 declareProperty( "UniField", m_uniField = false );
56
57 declareProperty( "Cur_SCQ1_55", m_Cur_SCQ1_55 = 349.4 );
58 declareProperty( "Cur_SCQ1_89", m_Cur_SCQ1_89 = 426.2 );
59 declareProperty( "Cur_SCQ2_10", m_Cur_SCQ2_10 = 474.2 );
60
61 declareProperty( "UseDBFlag", m_useDB = true );
62 declareProperty( "ReadOneTime", m_readOneTime = false );
63 declareProperty( "RunFrom", m_runFrom = 8093 );
64 declareProperty( "RunTo", m_runTo = 9025 );
65
66 m_Mucfield = new MucMagneticField();
67 if ( !m_Mucfield ) cout << "error: can not get MucMagneticField pointer" << endl;
68
69 m_zfield = -1.0 * tesla;
70
71 if ( getenv( "MAGNETICFIELDDATAROOT" ) )
72 { path = std::string( getenv( "MAGNETICFIELDDATAROOT" ) ); }
73 else { std::cerr << "Couldn't find MAGNETICFIELDDATAROOT" << std::endl; }
74}
75
76#else
77// -------------------------- BEAN ------------------------------------
78MagneticFieldSvc* MagneticFieldSvc::m_field = 0;
79
80//=============================================================================
81// Standard constructor, initializes variables
82//=============================================================================
83MagneticFieldSvc::MagneticFieldSvc() {
84 // use functions instead of "declareProperty"
85 m_gridDistance = 5;
86 m_runmode = 2;
87 m_ifRealField = true;
88 m_outlevel = 1;
89
90 m_Cur_SCQ1_55 = 349.4;
91 m_Cur_SCQ1_89 = 426.2;
92 m_Cur_SCQ2_10 = 474.2;
93
94 m_useDB = true;
95
96 path = "";
97 m_Mucfield = 0;
98
99 m_zfield = -1.0 * tesla;
100}
101#endif
102
103//=============================================================================
104// Standard destructor
105//=============================================================================
106MagneticFieldSvc::~MagneticFieldSvc() {
107 if ( m_Mucfield ) delete m_Mucfield;
108}
109
110//=============================================================================
111// Initialize
112//=============================================================================
113#ifdef BEAN
114bool init_mucMagneticField() {
115 // reinstall MucMagneticField probably with a new path
116 if ( m_Mucfield )
117 {
118 if ( m_Mucfield->getPath() == path ) return true;
119 delete m_Mucfield;
120 }
121
122 m_Mucfield = new ( std::nothrow ) MucMagneticField( path );
123 if ( !m_Mucfield )
124 {
125 cout << "error: can not get MucMagneticField pointer" << endl;
126 return false;
127 }
128
129 return true;
130}
131#endif
132
133StatusCode MagneticFieldSvc::initialize() {
134#ifndef BEAN
135 MsgStream log( msgSvc(), name() );
136 StatusCode status = Service::initialize();
137 former_m_filename_TE = "First Run";
138 former_m_filename = "First Run";
139 // cout<<"former_m_filename_TE is:"<<former_m_filename_TE<<":former_m_filename
140 // is:"<<former_m_filename<<endl;
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() ) { incsvc->addListener( this, "NewRun", priority ); }
146
147 status = serviceLocator()->service( "EventDataSvc", m_eventSvc, true );
148 if ( status.isFailure() )
149 {
150 log << MSG::ERROR << "Unable to find EventDataSvc " << endmsg;
151 return status;
152 }
153#else
154 if ( !init_mucMagneticField() ) return false;
155 StatusCode status = true;
156#endif
157 if ( m_useDB == false )
158 {
159 if ( m_gridDistance == 5 )
160 {
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 {
168 m_Q.reserve( 27082 );
169 m_P.reserve( 27082 );
170 m_Q_TE.reserve( 23833 );
171 m_P_TE.reserve( 23833 );
172 }
173
174 m_filename = path;
175 m_filename_TE = path;
176 if ( m_gridDistance == 10 )
177 {
178 m_filename_TE += std::string( "/dat/TEMap10cm.dat" );
179 if ( m_runmode == 2 )
180 m_filename += std::string( "/dat/2008-04-03BESIII-field-Mode2.dat" );
181 else if ( m_runmode == 4 )
182 m_filename += std::string( "/dat/2008-04-03BESIII-field-Mode4.dat" );
183 else if ( m_runmode == 6 )
184 m_filename += std::string( "/dat/2008-04-03BESIII-field-Mode6.dat" );
185 else if ( m_runmode == 7 )
186 m_filename += std::string( "/dat/2008-04-03BESIII-field-Mode7.dat" );
187 else
188 {
189 m_filename += std::string( "/dat/2007-08-28BESIII-field.dat" );
190 std::cout << "WARNING: YOU ARE USING OLD FIELD MAP!!!!" << std::endl;
191 }
192 }
193 if ( m_gridDistance == 5 )
194 {
195 m_filename_TE += std::string( "/dat/TEMap5cm.dat" );
196 if ( m_runmode == 1 )
197 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode1.dat" );
198 else if ( m_runmode == 2 )
199 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode2.dat" );
200 else if ( m_runmode == 3 )
201 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode3.dat" );
202 else if ( m_runmode == 4 )
203 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode4.dat" );
204 else if ( m_runmode == 5 )
205 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode5.dat" );
206 else if ( m_runmode == 6 )
207 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode6.dat" );
208 else if ( m_runmode == 7 )
209 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode7.dat" );
210 else if ( m_runmode == 8 )
211 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode8.dat" );
212 else
213 {
214 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode2.dat" );
215 std::cout << "WARNING: NO RUN MODE, YOU ARE USING DEFAULT FIELD MAP!!!!" << std::endl;
216 }
217 }
218 cout << "*** Read field map: " << m_filename << endl;
219 cout << "*** Read field map " << m_filename_TE << endl;
220
221 status = parseFile();
222 status = parseFile_TE();
223#ifndef BEAN
224 if ( status.isSuccess() )
225 {
226 log << MSG::DEBUG << "Magnetic field parsed successfully" << endmsg;
227#else
228 if ( status )
229 {
230 cout << "Magnetic field parsed successfully" << endl;
231#endif
232
233 for ( int iC = 0; iC < 2; ++iC )
234 {
235 m_min_FL[iC] = -90.0 * cm;
236 m_max_FL[iC] = m_min_FL[iC] + ( m_Nxyz[iC] - 1 ) * m_Dxyz[iC];
237 // for tof and emc
238 m_min_FL_TE[iC] = 0.0 * cm;
239 m_max_FL_TE[iC] = m_min_FL_TE[iC] + ( m_Nxyz_TE[iC] - 1 ) * m_Dxyz_TE[iC];
240 } // iC
241 m_min_FL[2] = -120.0 * cm;
242 m_max_FL[2] = m_min_FL[2] + ( m_Nxyz[2] - 1 ) * m_Dxyz[2];
243 // for tof and emc
244 m_min_FL_TE[2] = 0.0 * cm;
245 m_max_FL_TE[2] = m_min_FL_TE[2] + ( m_Nxyz_TE[2] - 1 ) * m_Dxyz_TE[2];
246 }
247 else
248 {
249#ifndef BEAN
250 log << MSG::DEBUG << "Magnetic field parse failled" << endmsg;
251#else
252 cout << "Magnetic field parse failled" << endl;
253#endif
254 }
255 }
256 else
257 {
258 m_connect_run = new FieldDBUtil::ConnectionDB();
259 if ( m_readOneTime == true )
260 {
261 bool e = m_connect_run->getReadSC_MagnetInfo( m_mapMagnetInfo, m_runFrom, m_runTo );
262 e = m_connect_run->getBeamEnergy( m_mapBeamEnergy, m_runFrom, m_runTo );
263 }
264#ifndef BEAN
265 if ( !m_connect_run )
266 { log << MSG::ERROR << "Could not open connection to database" << endmsg; }
267#endif
268 }
269 return status;
270}
271
272#ifndef BEAN
273void MagneticFieldSvc::init_params( std::vector<double> current,
274 std::vector<double> beamEnergy, int runNo ) {
275 MsgStream log( msgSvc(), name() );
276#else
277void init_params( int runNo ) {
278 if ( !init_mucMagneticField() )
279 {
280 cerr << " STOP! " << endl;
281 exit( 1 );
282 }
283#endif
284
285 m_Q.clear();
286 m_P.clear();
287 m_Q_1.clear();
288 m_P_1.clear();
289 m_Q_2.clear();
290 m_P_2.clear();
291 m_P_TE.clear();
292 m_Q_TE.clear();
293
294 if ( m_gridDistance == 5 )
295 {
296 m_Q.reserve( 201250 );
297 m_P.reserve( 201250 );
298 m_Q_1.reserve( 201250 );
299 m_P_1.reserve( 201250 );
300 m_Q_2.reserve( 201250 );
301 m_P_2.reserve( 201250 );
302 m_Q_TE.reserve( 176608 );
303 m_P_TE.reserve( 176608 );
304 }
305 if ( m_gridDistance == 10 )
306 {
307 m_Q.reserve( 27082 );
308 m_P.reserve( 27082 );
309 m_Q_1.reserve( 27082 );
310 m_P_1.reserve( 27082 );
311 m_Q_2.reserve( 27082 );
312 m_P_2.reserve( 27082 );
313 m_Q_TE.reserve( 23833 );
314 m_P_TE.reserve( 23833 );
315 }
316
317 char BPR_PRB[200];
318 char BER_PRB[200];
319 sprintf( BPR_PRB, "%f ", beamEnergy[0] );
320 sprintf( BER_PRB, "%f ", beamEnergy[1] );
321 setenv( "BEPCII_INFO.BPR_PRB", BPR_PRB, 1 );
322 setenv( "BEPCII_INFO.BER_PRB", BER_PRB, 1 );
323 int ssm_curr = int( current[0] );
324 double scql_curr = current[1];
325 double scqr_curr = current[2];
326 // cout<<"curent is:"<<ssm_curr<<"scql_curr is:"<<scql_curr<<"scqr_curr
327 // is:"<<scqr_currendl;
328
329#ifndef BEAN
330 log << MSG::INFO << "SSM current: " << current[0] << " SCQL current: " << current[1]
331 << " SCQR current: " << current[2] << " in Run " << runNo << endmsg;
332 log << MSG::INFO << "beamEnergy is:" << beamEnergy[0] << " BER_PRB:" << beamEnergy[1]
333 << " in Run " << runNo << endmsg;
334#else
335
336 cout << "SSM current: " << current[0] << " SCQL current: " << current[1]
337 << " SCQR current: " << current[2] << " in Run " << runNo << endl;
338 log << MSG::INFO << "beamEnergy is:" << beamEnergy[0] << " BER_PRB:" << beamEnergy[1]
339 << " in Run " << runNo << endmsg;
340#endif
341
342 // int ssm_curr = 3369;
343 // double scql_curr = 426.2;
344 // double scqr_curr = 426.2;
345 // for the energy less than 1.89GeV
346 if ( ( ( ssm_curr >= 3367 ) && ( ssm_curr <= 3370 ) &&
347 ( ( scql_curr + scqr_curr ) / 2 < m_Cur_SCQ1_89 ) ) )
348 {
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 << endmsg;
354#else
355 cout << "Select field map: " << m_filename << endl;
356#endif
357
358 if ( ( ( former_m_filename == "First Run" ) || ( former_m_filename != m_filename ) ) &&
359 ( m_readOneTime == true ) )
360 { former_m_filename = m_filename; }
361 else if ( m_readOneTime == false ) {}
362 StatusCode status = parseFile();
363#ifndef BEAN
364 if ( !status.isSuccess() )
365 { log << MSG::DEBUG << "Magnetic field parse failled" << endmsg; }
366#else
367 if ( !status ) { cout << "Magnetic field parse failled" << endl; }
368#endif
369
370 m_Q_1 = m_Q;
371 m_P_1 = m_P;
372
373 m_Q.clear();
374 m_P.clear();
375 if ( m_gridDistance == 5 )
376 {
377 m_Q.reserve( 201250 );
378 m_P.reserve( 201250 );
379 }
380 if ( m_gridDistance == 10 )
381 {
382 m_Q.reserve( 27082 );
383 m_P.reserve( 27082 );
384 }
385
386 m_filename = path;
387 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode3.dat" );
388#ifndef BEAN
389 log << MSG::INFO << "Select field map: " << m_filename << endmsg;
390#else
391 cout << "Select field map: " << m_filename << endl;
392#endif
393
394 if ( ( ( former_m_filename == "First Run" ) || ( former_m_filename != m_filename ) ) &&
395 ( m_readOneTime == true ) )
396 { former_m_filename = m_filename; }
397 else if ( m_readOneTime == false ) {}
398 status = parseFile();
399#ifndef BEAN
400 if ( !status.isSuccess() )
401 { log << MSG::DEBUG << "Magnetic field parse failled" << endmsg; }
402#else
403 if ( !status ) { cout << "Magnetic field parse failled" << endl; }
404#endif
405
406 m_Q_2 = m_Q;
407 m_P_2 = m_P;
408
409 m_Q.clear();
410 m_P.clear();
411 if ( m_gridDistance == 5 )
412 {
413 m_Q.reserve( 201250 );
414 m_P.reserve( 201250 );
415 }
416 if ( m_gridDistance == 10 )
417 {
418 m_Q.reserve( 27082 );
419 m_P.reserve( 27082 );
420 }
421 // check
422 if ( m_Q_2.size() != m_Q_1.size() )
423 {
424#ifndef BEAN
425 log << MSG::FATAL << "The two field maps used in this run are wrong!!!" << endmsg;
426#else
427 cout << "The two field maps used in this run are wrong!!!" << endl;
428#endif
429 exit( 1 );
430 }
431
432 for ( unsigned int iQ = 0; iQ < m_Q_1.size(); iQ++ )
433 {
434 double fieldvalue = ( m_Q_1[iQ] - m_Q_2[iQ] ) / ( m_Cur_SCQ1_55 - m_Cur_SCQ1_89 ) *
435 ( ( scql_curr + scqr_curr ) / 2 - m_Cur_SCQ1_89 ) +
436 m_Q_2[iQ];
437 m_Q.push_back( fieldvalue );
438 m_P.push_back( m_P_2[iQ] );
439 }
440 }
441 // for the energy greater than 1.89GeV
442 if ( ( ( ssm_curr >= 3367 ) && ( ssm_curr <= 3370 ) &&
443 ( ( scql_curr + scqr_curr ) / 2 >= m_Cur_SCQ1_89 ) ) )
444 {
445 m_zfield = -1.0 * tesla;
446 m_filename = path;
447 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode3.dat" );
448#ifndef BEAN
449 log << MSG::INFO << "Select field map: " << m_filename << endmsg;
450#else
451 cout << "Select field map: " << m_filename << endl;
452#endif
453 if ( ( ( former_m_filename == "First Run" ) || ( former_m_filename != m_filename ) ) &&
454 ( m_readOneTime == true ) )
455 { former_m_filename = m_filename; }
456 else if ( m_readOneTime == false ) {}
457 StatusCode status = parseFile();
458#ifndef BEAN
459 if ( !status.isSuccess() )
460 { log << MSG::DEBUG << "Magnetic field parse failled" << endmsg; }
461#else
462 if ( !status ) { cout << "Magnetic field parse failled" << endl; }
463#endif
464
465 m_Q_1 = m_Q;
466 m_P_1 = m_P;
467
468 m_Q.clear();
469 m_P.clear();
470 if ( m_gridDistance == 5 )
471 {
472 m_Q.reserve( 201250 );
473 m_P.reserve( 201250 );
474 }
475 if ( m_gridDistance == 10 )
476 {
477 m_Q.reserve( 27082 );
478 m_P.reserve( 27082 );
479 }
480
481 m_filename = path;
482 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode4.dat" );
483#ifndef BEAN
484 log << MSG::INFO << "Select field map: " << m_filename << endmsg;
485#else
486 cout << "Select field map: " << m_filename << endl;
487#endif
488
489 if ( ( ( former_m_filename == "First Run" ) || ( former_m_filename != m_filename ) ) &&
490 ( m_readOneTime == true ) )
491 { former_m_filename = m_filename; }
492 else if ( m_readOneTime == false ) {}
493 status = parseFile();
494#ifndef BEAN
495 if ( !status.isSuccess() )
496 { log << MSG::DEBUG << "Magnetic field parse failled" << endmsg; }
497#else
498 if ( !status ) { cout << "Magnetic field parse failled" << endl; }
499#endif
500
501 m_Q_2 = m_Q;
502 m_P_2 = m_P;
503
504 m_Q.clear();
505 m_P.clear();
506 if ( m_gridDistance == 5 )
507 {
508 m_Q.reserve( 201250 );
509 m_P.reserve( 201250 );
510 }
511 if ( m_gridDistance == 10 )
512 {
513 m_Q.reserve( 27082 );
514 m_P.reserve( 27082 );
515 }
516
517 // check
518 if ( m_Q_2.size() != m_Q_1.size() )
519 {
520#ifndef BEAN
521
522 log << MSG::FATAL << "The two field maps used in this run are wrong!!!" << endmsg;
523
524#else
525
526 cout << "The two field maps used in this run are wrong!!!" << endl;
527
528#endif
529 exit( 1 );
530 }
531
532 for ( unsigned int iQ = 0; iQ < m_Q_1.size(); iQ++ )
533 {
534 double fieldvalue = ( m_Q_1[iQ] - m_Q_2[iQ] ) / ( m_Cur_SCQ1_89 - m_Cur_SCQ2_10 ) *
535 ( ( scql_curr + scqr_curr ) / 2 - m_Cur_SCQ2_10 ) +
536 m_Q_2[iQ];
537 m_Q.push_back( fieldvalue );
538 m_P.push_back( m_P_2[iQ] );
539 }
540 }
541 if ( ( ssm_curr >= 3033 ) && ( ssm_curr <= 3035 ) )
542 {
543 m_zfield = -0.9 * tesla;
544 m_filename = path;
545 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode7.dat" );
546#ifndef BEAN
547 log << MSG::INFO << "Select field map: " << m_filename << endmsg;
548#else
549 cout << "Select field map: " << m_filename << endl;
550#endif
551
552 if ( ( ( former_m_filename == "First Run" ) || ( former_m_filename != m_filename ) ) &&
553 ( m_readOneTime == true ) )
554 { former_m_filename = m_filename; }
555 else if ( m_readOneTime == false ) {}
556 StatusCode status = parseFile();
557#ifndef BEAN
558 if ( !status.isSuccess() )
559 { log << MSG::DEBUG << "Magnetic field parse failled" << endmsg; }
560#else
561 if ( !status ) { cout << "Magnetic field parse failled" << endl; }
562#endif
563
564 m_Q_1 = m_Q;
565 m_P_1 = m_P;
566
567 m_Q.clear();
568 m_P.clear();
569 if ( m_gridDistance == 5 )
570 {
571 m_Q.reserve( 201250 );
572 m_P.reserve( 201250 );
573 }
574 if ( m_gridDistance == 10 )
575 {
576 m_Q.reserve( 27082 );
577 m_P.reserve( 27082 );
578 }
579
580 m_filename = path;
581 m_filename += std::string( "/dat/2008-05-27BESIII-field-Mode8.dat" );
582#ifndef BEAN
583 log << MSG::INFO << "Select field map: " << m_filename << endmsg;
584#else
585 cout << "Select field map: " << m_filename << endl;
586#endif
587
588 if ( ( ( former_m_filename == "First Run" ) || ( former_m_filename != m_filename ) ) &&
589 ( m_readOneTime == true ) )
590 { former_m_filename = m_filename; }
591 else if ( m_readOneTime == false ) {}
592 status = parseFile();
593#ifndef BEAN
594 if ( !status.isSuccess() )
595 { log << MSG::DEBUG << "Magnetic field parse failled" << endmsg; }
596#else
597 if ( !status ) { cout << "Magnetic field parse failled" << endl; }
598#endif
599
600 m_Q_2 = m_Q;
601 m_P_2 = m_P;
602
603 m_Q.clear();
604 m_P.clear();
605 if ( m_gridDistance == 5 )
606 {
607 m_Q.reserve( 201250 );
608 m_P.reserve( 201250 );
609 }
610 if ( m_gridDistance == 10 )
611 {
612 m_Q.reserve( 27082 );
613 m_P.reserve( 27082 );
614 }
615 // check
616 if ( m_Q_2.size() != m_Q_1.size() )
617 {
618#ifndef BEAN
619 log << MSG::FATAL << "The two field maps used in this run are wrong!!!" << endmsg;
620#else
621 cout << "The two field maps used in this run are wrong!!!" << endl;
622#endif
623 exit( 1 );
624 }
625
626 for ( unsigned int iQ = 0; iQ < m_Q_1.size(); iQ++ )
627 {
628 double fieldvalue = ( m_Q_1[iQ] - m_Q_2[iQ] ) / ( m_Cur_SCQ1_55 - m_Cur_SCQ1_89 ) *
629 ( ( scql_curr + scqr_curr ) / 2 - m_Cur_SCQ1_89 ) +
630 m_Q_2[iQ];
631 m_Q.push_back( fieldvalue );
632 m_P.push_back( m_P_2[iQ] );
633 }
634 }
635 if ( ( !( ( ssm_curr >= 3367 ) && ( ssm_curr <= 3370 ) ) &&
636 !( ( ssm_curr >= 3033 ) && ( ssm_curr <= 3035 ) ) ) ||
637 scql_curr == 0 || scqr_curr == 0 )
638 {
639#ifndef BEAN
640 log << MSG::FATAL << "The current of run " << runNo
641 << " is abnormal in database, please check it! " << endmsg;
642 log << MSG::FATAL << "The current of SSM is " << ssm_curr << ", SCQR is " << scqr_curr
643 << ", SCQL is " << scql_curr << endmsg;
644#else
645 cout << "The current of run " << runNo << " is abnormal in database, please check it! "
646 << endl;
647 cout << "The current of SSM is " << ssm_curr << ", SCQR is " << scqr_curr << ", SCQL is "
648 << scql_curr << endl;
649#endif
650 exit( 1 );
651 }
652 if ( m_Q.size() == 0 )
653 {
654#ifndef BEAN
655 log << MSG::FATAL
656 << "This size of field map vector is ZERO, please check the current of run " << runNo
657 << " in database!" << endmsg;
658#else
659 cout << "This size of field map vector is ZERO,"
660 << " please check the current of run " << runNo << " in database!" << endl;
661#endif
662 exit( 1 );
663 }
664
665 m_filename_TE = path;
666 if ( m_gridDistance == 10 ) m_filename_TE += std::string( "/dat/TEMap10cm.dat" );
667 if ( m_gridDistance == 5 ) m_filename_TE += std::string( "/dat/TEMap5cm.dat" );
668#ifndef BEAN
669 log << MSG::INFO << "Select field map: " << m_filename_TE << endmsg;
670#else
671 cout << "Select field map: " << m_filename_TE << endl;
672#endif
673 if ( ( ( former_m_filename_TE == "First Run" ) ||
674 ( former_m_filename_TE != m_filename_TE ) ) &&
675 ( m_readOneTime == true ) )
676 { former_m_filename_TE = m_filename_TE; }
677 else if ( m_readOneTime == false ) {}
678 StatusCode status = parseFile_TE();
679#ifndef BEAN
680 if ( !status.isSuccess() ) { log << MSG::DEBUG << "Magnetic field parse failled" << endmsg; }
681#else
682 if ( !status ) { cout << "Magnetic field parse failled" << endl; }
683#endif
684
685 for ( int iC = 0; iC < 2; ++iC )
686 {
687 m_min_FL[iC] = -90.0 * cm;
688 m_max_FL[iC] = m_min_FL[iC] + ( m_Nxyz[iC] - 1 ) * m_Dxyz[iC];
689 // for tof and emc
690 m_min_FL_TE[iC] = 0.0 * cm;
691 m_max_FL_TE[iC] = m_min_FL_TE[iC] + ( m_Nxyz_TE[iC] - 1 ) * m_Dxyz_TE[iC];
692 } // iC
693 m_min_FL[2] = -120.0 * cm;
694 m_max_FL[2] = m_min_FL[2] + ( m_Nxyz[2] - 1 ) * m_Dxyz[2];
695 // for tof and emc
696 m_min_FL_TE[2] = 0.0 * cm;
697 m_max_FL_TE[2] = m_min_FL_TE[2] + ( m_Nxyz_TE[2] - 1 ) * m_Dxyz_TE[2];
698}
699
700#ifndef BEAN
701//=============================================================================
702// Finalize
703//=============================================================================
704StatusCode MagneticFieldSvc::finalize() {
705 MsgStream log( msgSvc(), name() );
706 // if(m_connect_run) delete m_connect_run;
707 StatusCode status = Service::finalize();
708
709 if ( status.isSuccess() ) log << MSG::INFO << "Service finalized successfully" << endmsg;
710 return status;
711}
712
713//=============================================================================
714// QueryInterface
715//=============================================================================
716/*StatusCode MagneticFieldSvc::queryInterface( const InterfaceID& riid,
717 void** ppvInterface )
718{
719 StatusCode sc = StatusCode::FAILURE;
720 if ( ppvInterface ) {
721 *ppvInterface = 0;
722
723 if ( riid == IID_IMagneticFieldSvc ) {
724 *ppvInterface = static_cast<IMagneticFieldSvc*>(this);
725 sc = StatusCode::SUCCESS;
726 addRef();
727 }
728 else
729 sc = Service::queryInterface( riid, ppvInterface );
730 }
731 return sc;
732}
733*/
734void MagneticFieldSvc::handle( const Incident& inc ) {
735 MsgStream log( msgSvc(), name() );
736 log << MSG::DEBUG << "handle: " << inc.type() << endmsg;
737 if ( inc.type() != "NewRun" ) { return; }
738 log << MSG::DEBUG << "Begin New Runcc" << endmsg;
739
740 SmartDataPtr<Event::EventHeader> eventHeader( m_eventSvc, "/Event/EventHeader" );
741 int new_run = eventHeader->runNumber();
742 if ( new_run < 0 ) new_run = -new_run;
743
744 if ( ( m_useDB == true ) && ( m_readOneTime == false ) )
745 {
746 using FieldDBUtil::ConnectionDB;
747 ConnectionDB::eRet e = m_connect_run->getReadSC_MagnetInfo( current, std::abs( new_run ) );
748 e = m_connect_run->getBeamEnergy( beamEnergy, std::abs( new_run ) );
749 init_params( current, beamEnergy, new_run );
750 }
751 else if ( ( m_useDB == true ) && ( m_readOneTime == true ) )
752 {
753 std::vector<double> beamEnergy;
754 std::vector<double> current( 3, 0. );
755 // if (m_mapBeamEnergy == NULL ||m_mapBeamEnergy.size() == 0))
756 if ( m_mapBeamEnergy.size() == 0 )
757 {
758 cout << "Run:" << new_run << " BeamEnergy is NULL, please check!!" << endl;
759 exit( 1 );
760 }
761 beamEnergy.push_back( ( m_mapBeamEnergy[new_run] )[0] );
762 beamEnergy.push_back( ( m_mapBeamEnergy[new_run] )[1] );
763 // if (m_mapMagnetInfo.size()<1||m_mapMagnetInfo.isEmpty())
764 if ( m_mapMagnetInfo.size() == 0 )
765 {
766 cout << "Run:" << new_run << " MagnetInfo is NULL, please check!!" << endl;
767 exit( 1 );
768 }
769 current[0] = ( m_mapMagnetInfo[new_run] )[0];
770 current[1] = ( m_mapMagnetInfo[new_run] )[1];
771 current[2] = ( m_mapMagnetInfo[new_run] )[2];
772 init_params( current, beamEnergy, new_run );
773 }
774}
775
776#else
777void MagneticFieldSvc::handle( int new_run ) {
778 static int save_run = 0;
779 if ( new_run == save_run ) return;
780
781 cout << "Begin New Run " << new_run << endl;
782 save_run = new_run;
783 if ( m_useDB == true ) { init_params( current, beamEnergy, new_run ); }
784}
785#endif
786
787// ---------------------------------------------------------------------------
788// Routine: parseFile
789// Purpose: Parses the file and fill a magnetic field vector
790// ---------------------------------------------------------------------------
791StatusCode MagneticFieldSvc::parseFile() {
792#ifndef BEAN
793 StatusCode sc = StatusCode::FAILURE;
794
795 MsgStream log( msgSvc(), name() );
796#else
797 StatusCode sc = false;
798#endif
799
800 char line[255];
801 std::ifstream infile( m_filename.c_str() );
802 if ( infile )
803 {
804#ifndef BEAN
805 sc = StatusCode::SUCCESS;
806#else
807 sc = true;
808#endif
809
810 // Skip the header till PARAMETER
811 do {
812 infile.getline( line, 255 );
813 } while ( line[0] != 'P' );
814
815 // Get the PARAMETER
816 std::string sPar[2];
817 char* token = strtok( line, " " );
818 int ip = 0;
819 do {
820 if ( token )
821 {
822 sPar[ip] = token;
823 token = strtok( NULL, " " );
824 }
825 else continue;
826 ip++;
827 } while ( token != NULL );
828 long int npar = atoi( sPar[1].c_str() );
829
830 // Skip the header till GEOMETRY
831 do {
832 infile.getline( line, 255 );
833 } while ( line[0] != 'G' );
834
835 // Skip any comment before GEOMETRY
836 do {
837 infile.getline( line, 255 );
838 } while ( line[0] != '#' );
839
840 // Get the GEOMETRY
841 infile.getline( line, 255 );
842 std::string sGeom[7];
843 token = strtok( line, " " );
844 int ig = 0;
845 do {
846 if ( token )
847 {
848 sGeom[ig] = token;
849 token = strtok( NULL, " " );
850 }
851 else continue;
852 ig++;
853 } while ( token != NULL );
854
855 // Grid dimensions are given in cm in CDF file. Convert to CLHEP units
856 m_Dxyz[0] = atof( sGeom[0].c_str() ) * cm;
857 m_Dxyz[1] = atof( sGeom[1].c_str() ) * cm;
858 m_Dxyz[2] = atof( sGeom[2].c_str() ) * cm;
859 m_Nxyz[0] = atoi( sGeom[3].c_str() );
860 m_Nxyz[1] = atoi( sGeom[4].c_str() );
861 m_Nxyz[2] = atoi( sGeom[5].c_str() );
862 m_zOffSet = atof( sGeom[6].c_str() ) * cm;
863 // Number of lines with data to be read
864 long int nlines = ( npar - 7 ) / 3;
865
866 // Check number of lines with data read in the loop
867 long int ncheck = 0;
868
869 // Skip comments and fill a vector of magnetic components for the
870 // x, y and z positions given in GEOMETRY
871
872 while ( infile )
873 {
874 // parse each line of the file,
875 // comment lines begin with '#' in the cdf file
876 infile.getline( line, 255 );
877 if ( line[0] == '#' ) continue;
878 std::string gridx, gridy, gridz, sFx, sFy, sFz;
879 char* token = strtok( line, " " );
880 if ( token )
881 {
882 gridx = token;
883 token = strtok( NULL, " " );
884 }
885 else continue;
886 if ( token )
887 {
888 gridy = token;
889 token = strtok( NULL, " " );
890 }
891 else continue;
892 if ( token )
893 {
894 gridz = token;
895 token = strtok( NULL, " " );
896 }
897 else continue;
898 if ( token )
899 {
900 sFx = token;
901 token = strtok( NULL, " " );
902 }
903 else continue;
904 if ( token )
905 {
906 sFy = token;
907 token = strtok( NULL, " " );
908 }
909 else continue;
910 if ( token )
911 {
912 sFz = token;
913 token = strtok( NULL, " " );
914 }
915 else continue;
916 if ( token != NULL ) continue;
917 // Grid position
918 double gx = atof( gridx.c_str() ) * m;
919 double gy = atof( gridy.c_str() ) * m;
920 double gz = atof( gridz.c_str() ) * m;
921 // Field values are given in tesla in CDF file. Convert to CLHEP units
922 double fx = atof( sFx.c_str() ) * tesla;
923 double fy = atof( sFy.c_str() ) * tesla;
924 double fz = atof( sFz.c_str() ) * tesla;
925 // for debug
926 if ( m_outlevel == 0 )
927 {
928#ifndef BEAN
929 log << MSG::DEBUG << "grid x, y, z is " << gx << ", " << gy << ", " << gz << endmsg;
930 log << MSG::DEBUG << "field x, y, z is " << fx << ", " << fy << ", " << fz << endmsg;
931#else
932 cout << "grid x, y, z is " << gx << ", " << gy << ", " << gz << endl;
933 cout << "field x, y, z is " << fx << ", " << fy << ", " << fz << endl;
934#endif
935 }
936 // Fill the postion to the vector
937 m_P.push_back( gx );
938 m_P.push_back( gy );
939 m_P.push_back( gz );
940 // Add the magnetic field components of each point to
941 // sequentialy in a vector
942 m_Q.push_back( fx );
943 m_Q.push_back( fy );
944 m_Q.push_back( fz );
945 // counts after reading and filling to match the number of lines
946 ncheck++;
947 }
948 infile.close();
949 if ( nlines != ncheck )
950 {
951#ifndef BEAN
952 log << MSG::ERROR << "nline is " << nlines << "; ncheck is " << ncheck
953 << " Number of points in field map does not match!" << endmsg;
954 return StatusCode::FAILURE;
955#else
956 cout << "nline is " << nlines << "; ncheck is " << ncheck
957 << " Number of points in field map does not match!" << endl;
958 return false;
959#endif
960 }
961 }
962 else
963 {
964#ifndef BEAN
965 log << MSG::ERROR << __FILE__ ":" << __LINE__
966 << ": Unable to open magnetic field file : " << m_filename << endmsg;
967#else
968 cout << "Unable to open magnetic field file : " << m_filename << endl;
969#endif
970 }
971
972 return sc;
973}
974
975// ---------------------------------------------------------------------------
976// Routine: parseFile_TE
977// Purpose: Parses the file and fill a magnetic field vector for tof and emc
978// ---------------------------------------------------------------------------
979StatusCode MagneticFieldSvc::parseFile_TE() {
980#ifndef BEAN
981 StatusCode sc = StatusCode::FAILURE;
982
983 MsgStream log( msgSvc(), name() );
984#else
985 StatusCode sc = false;
986#endif
987
988 char line[255];
989 std::ifstream infile( m_filename_TE.c_str() );
990
991 if ( infile )
992 {
993#ifndef BEAN
994 sc = StatusCode::SUCCESS;
995#else
996 sc = true;
997#endif
998
999 // Skip the header till PARAMETER
1000 do {
1001 infile.getline( line, 255 );
1002 } while ( line[0] != 'P' );
1003
1004 // Get the PARAMETER
1005 std::string sPar[2];
1006 char* token = strtok( line, " " );
1007 int ip = 0;
1008 do {
1009 if ( token )
1010 {
1011 sPar[ip] = token;
1012 token = strtok( NULL, " " );
1013 }
1014 else continue;
1015 ip++;
1016 } while ( token != NULL );
1017 long int npar = atoi( sPar[1].c_str() );
1018
1019 // Skip the header till GEOMETRY
1020 do {
1021 infile.getline( line, 255 );
1022 } while ( line[0] != 'G' );
1023
1024 // Skip any comment before GEOMETRY
1025 do {
1026 infile.getline( line, 255 );
1027 } while ( line[0] != '#' );
1028
1029 // Get the GEOMETRY
1030 infile.getline( line, 255 );
1031 std::string sGeom[7];
1032 token = strtok( line, " " );
1033 int ig = 0;
1034 do {
1035 if ( token )
1036 {
1037 sGeom[ig] = token;
1038 token = strtok( NULL, " " );
1039 }
1040 else continue;
1041 ig++;
1042 } while ( token != NULL );
1043
1044 // Grid dimensions are given in cm in CDF file. Convert to CLHEP units
1045 m_Dxyz_TE[0] = atof( sGeom[0].c_str() ) * cm;
1046 m_Dxyz_TE[1] = atof( sGeom[1].c_str() ) * cm;
1047 m_Dxyz_TE[2] = atof( sGeom[2].c_str() ) * cm;
1048 m_Nxyz_TE[0] = atoi( sGeom[3].c_str() );
1049 m_Nxyz_TE[1] = atoi( sGeom[4].c_str() );
1050 m_Nxyz_TE[2] = atoi( sGeom[5].c_str() );
1051 m_zOffSet_TE = atof( sGeom[6].c_str() ) * cm;
1052 // Number of lines with data to be read
1053 long int nlines = ( npar - 7 ) / 3;
1054
1055 // Check number of lines with data read in the loop
1056 long int ncheck = 0;
1057
1058 // Skip comments and fill a vector of magnetic components for the
1059 // x, y and z positions given in GEOMETRY
1060
1061 while ( infile )
1062 {
1063 // parse each line of the file,
1064 // comment lines begin with '#' in the cdf file
1065 infile.getline( line, 255 );
1066 if ( line[0] == '#' ) continue;
1067 std::string gridx, gridy, gridz, sFx, sFy, sFz;
1068 char* token = strtok( line, " " );
1069 if ( token )
1070 {
1071 gridx = token;
1072 token = strtok( NULL, " " );
1073 }
1074 else continue;
1075 if ( token )
1076 {
1077 gridy = token;
1078 token = strtok( NULL, " " );
1079 }
1080 else continue;
1081 if ( token )
1082 {
1083 gridz = token;
1084 token = strtok( NULL, " " );
1085 }
1086 else continue;
1087 if ( token )
1088 {
1089 sFx = token;
1090 token = strtok( NULL, " " );
1091 }
1092 else continue;
1093 if ( token )
1094 {
1095 sFy = token;
1096 token = strtok( NULL, " " );
1097 }
1098 else continue;
1099 if ( token )
1100 {
1101 sFz = token;
1102 token = strtok( NULL, " " );
1103 }
1104 else continue;
1105 if ( token != NULL ) continue;
1106 // Grid position
1107 double gx = atof( gridx.c_str() ) * m;
1108 double gy = atof( gridy.c_str() ) * m;
1109 double gz = atof( gridz.c_str() ) * m;
1110 // Field values are given in tesla in CDF file. Convert to CLHEP units
1111 double fx = atof( sFx.c_str() ) * tesla;
1112 double fy = atof( sFy.c_str() ) * tesla;
1113 double fz = atof( sFz.c_str() ) * tesla;
1114 // for debug
1115 if ( m_outlevel == 0 )
1116 {
1117#ifndef BEAN
1118 log << MSG::DEBUG << "grid x, y, z is " << gx << ", " << gy << ", " << gz << endmsg;
1119 log << MSG::DEBUG << "field x, y, z is " << fx << ", " << fy << ", " << fz << endmsg;
1120#else
1121 cout << "grid x, y, z is " << gx << ", " << gy << ", " << gz << endl;
1122 cout << "field x, y, z is " << fx << ", " << fy << ", " << fz << endl;
1123#endif
1124 }
1125 // Fill the postion to the vector
1126 m_P_TE.push_back( gx );
1127 m_P_TE.push_back( gy );
1128 m_P_TE.push_back( gz );
1129 // Add the magnetic field components of each point to
1130 // sequentialy in a vector
1131 m_Q_TE.push_back( fx );
1132 m_Q_TE.push_back( fy );
1133 m_Q_TE.push_back( fz );
1134 // counts after reading and filling to match the number of lines
1135 ncheck++;
1136 }
1137 infile.close();
1138 if ( nlines != ncheck )
1139 {
1140#ifndef BEAN
1141 log << MSG::ERROR << "nline is " << nlines << "; ncheck is " << ncheck
1142 << " Number of points in field map does not match!" << endmsg;
1143 return StatusCode::FAILURE;
1144#else
1145 cout << "nline is " << nlines << "; ncheck is " << ncheck
1146 << " Number of points in field map does not match!" << endl;
1147 return false;
1148#endif
1149 }
1150 }
1151 else
1152 {
1153#ifndef BEAN
1154 log << MSG::ERROR << __FILE__ << ":" << __LINE__
1155 << ": Unable to open magnetic field file : " << m_filename_TE << endmsg;
1156#else
1157 cout << "Unable to open magnetic field file : " << m_filename_TE << endl;
1158#endif
1159 }
1160
1161 return sc;
1162}
1163//=============================================================================
1164// FieldVector: find the magnetic field value at a given point in space
1165//=============================================================================
1166StatusCode MagneticFieldSvc::fieldVector( const HepPoint3D& newr, HepVector3D& newb ) const {
1167
1168 if ( m_turnOffField == true )
1169 {
1170 newb[0] = 0.;
1171 newb[1] = 0.;
1172 newb[2] = 0.;
1173#ifndef BEAN
1174 return StatusCode::SUCCESS;
1175#else
1176 return true;
1177#endif
1178 }
1179
1180 // wll added 2012-08-27
1181 if ( m_uniField )
1182 {
1183 uniFieldVector( newr, newb );
1184#ifndef BEAN
1185 return StatusCode::SUCCESS;
1186#else
1187 return true;
1188#endif
1189 }
1190
1191 // reference frames defined by simulation and SCM are different. In simulation: x--south to
1192 // north, y--down to up, but in SCM: x--north to south, y--down to up
1193 double new_x = -newr.x();
1194 double new_y = newr.y();
1195 double new_z = -newr.z();
1196 HepPoint3D r( new_x, new_y, new_z );
1197 HepVector3D b;
1198 b[0] = 0.0 * tesla;
1199 b[1] = 0.0 * tesla;
1200 b[2] = 0.0 * tesla;
1201 // This routine is now dummy. Was previously converting to/from CLHEP units
1202 if ( -2.1 * m < r.z() && r.z() < 2.1 * m && -1.8 * m < r.x() && r.x() < 1.8 * m &&
1203 -1.8 * m < r.y() && r.y() < 1.8 * m )
1204 {
1205 if ( -1.2 * m < r.z() && r.z() < 1.2 * m &&
1206 0 * m <= std::sqrt( r.x() * r.x() + r.y() * r.y() ) &&
1207 std::sqrt( r.x() * r.x() + r.y() * r.y() ) < 0.9 * m )
1208 // if(std::abs(r.z())<1.2*m&&std::abs(r.x())<=0.9*m&&std::abs(r.y())<=0.9*m)
1209 { this->fieldGrid( r, b ); }
1210 else { this->fieldGrid_TE( r, b ); }
1211 }
1212 if ( ( fabs( r.z() ) <= 1970 * mm && sqrt( r.x() * r.x() + r.y() * r.y() ) >= 1740 * mm ) ||
1213 ( fabs( r.z() ) >= 2050 * mm ) )
1214 {
1215 HepPoint3D mr;
1216 HepVector3D tb;
1217 int part = 0, layer = 0, mat = 0;
1218 double theta;
1219 bool ifbar = false;
1220 bool ifend = false;
1221 if ( r.x() != 0. )
1222 {
1223 theta = atan2( fabs( r.y() ), fabs( r.x() ) );
1224 if ( 0 <= theta && theta < pi / 8 )
1225 {
1226 mr[0] = fabs( r.x() );
1227 mr[1] = -fabs( r.y() );
1228 mr[2] = fabs( r.z() );
1229 if ( mr[2] <= 1970 * mm && 1740 * mm <= mr[0] && mr[0] <= 2620 * mm )
1230 {
1231 part = 1;
1232 if ( 1740 * mm <= mr[0] && mr[0] < 1770 * mm )
1233 {
1234 layer = 0;
1235 mat = 0;
1236 }
1237 if ( 1770 * mm <= mr[0] && mr[0] < 1810 * mm )
1238 {
1239 layer = 0;
1240 mat = 1;
1241 }
1242 if ( 1810 * mm <= mr[0] && mr[0] < 1840 * mm )
1243 {
1244 layer = 1;
1245 mat = 0;
1246 }
1247 if ( 1840 * mm <= mr[0] && mr[0] < 1880 * mm )
1248 {
1249 layer = 1;
1250 mat = 1;
1251 }
1252 if ( 1880 * mm <= mr[0] && mr[0] < 1910 * mm )
1253 {
1254 layer = 2;
1255 mat = 0;
1256 }
1257 if ( 1910 * mm <= mr[0] && mr[0] < 1950 * mm )
1258 {
1259 layer = 2;
1260 mat = 1;
1261 }
1262 if ( 1950 * mm <= mr[0] && mr[0] < 1990 * mm )
1263 {
1264 layer = 3;
1265 mat = 0;
1266 }
1267 if ( 1990 * mm <= mr[0] && mr[0] < 2030 * mm )
1268 {
1269 layer = 3;
1270 mat = 1;
1271 }
1272 if ( 2030 * mm <= mr[0] && mr[0] < 2070 * mm )
1273 {
1274 layer = 4;
1275 mat = 0;
1276 }
1277 if ( 2070 * mm <= mr[0] && mr[0] < 2110 * mm )
1278 {
1279 layer = 4;
1280 mat = 1;
1281 }
1282 if ( 2110 * mm <= mr[0] && mr[0] < 2190 * mm )
1283 {
1284 layer = 5;
1285 mat = 0;
1286 }
1287 if ( 2190 * mm <= mr[0] && mr[0] < 2230 * mm )
1288 {
1289 layer = 5;
1290 mat = 1;
1291 }
1292 if ( 2230 * mm <= mr[0] && mr[0] < 2310 * mm )
1293 {
1294 layer = 6;
1295 mat = 0;
1296 }
1297 if ( 2310 * mm <= mr[0] && mr[0] < 2350 * mm )
1298 {
1299 layer = 6;
1300 mat = 1;
1301 }
1302 if ( 2350 * mm <= mr[0] && mr[0] < 2430 * mm )
1303 {
1304 layer = 7;
1305 mat = 0;
1306 }
1307 if ( 2430 * mm <= mr[0] && mr[0] < 2470 * mm )
1308 {
1309 layer = 7;
1310 mat = 1;
1311 }
1312 if ( 2470 * mm <= mr[0] && mr[0] <= 2620 * mm )
1313 {
1314 layer = 8;
1315 mat = 0;
1316 }
1317 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1318 b[0] = tb[0];
1319 b[1] = -tb[1];
1320 b[2] = tb[2];
1321 ifbar = true;
1322 }
1323 if ( 2050 * mm <= mr[2] && mr[2] < 2090 * mm && 1034 * mm <= mr[0] &&
1324 mr[0] <= 2500 * mm )
1325 {
1326 part = 0;
1327 layer = 0;
1328 mat = 0;
1329 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1330 b[0] = tb[0];
1331 b[1] = -tb[1];
1332 b[2] = tb[2];
1333 ifend = true;
1334 }
1335 if ( 2090 * mm <= mr[2] && mr[2] < 2130 * mm && 1067 * mm <= mr[0] &&
1336 mr[0] <= 2500 * mm )
1337 {
1338 part = 0;
1339 layer = 0;
1340 mat = 1;
1341 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1342 b[0] = tb[0];
1343 b[1] = -tb[1];
1344 b[2] = tb[2];
1345 ifend = true;
1346 }
1347 if ( 2130 * mm <= mr[2] && mr[2] < 2170 * mm && 1067 * mm <= mr[0] &&
1348 mr[0] <= 2500 * mm )
1349 {
1350 part = 0;
1351 layer = 1;
1352 mat = 0;
1353 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1354 b[0] = tb[0];
1355 b[1] = -tb[1];
1356 b[2] = tb[2];
1357 ifend = true;
1358 }
1359 if ( 2170 * mm <= mr[2] && mr[2] < 2210 * mm && 1100 * mm <= mr[0] &&
1360 mr[0] <= 2500 * mm )
1361 {
1362 part = 0;
1363 layer = 1;
1364 mat = 1;
1365 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1366 b[0] = tb[0];
1367 b[1] = -tb[1];
1368 b[2] = tb[2];
1369 ifend = true;
1370 }
1371 if ( 2210 * mm <= mr[2] && mr[2] < 2240 * mm && 1100 * mm < mr[0] &&
1372 mr[0] <= 2500 * mm )
1373 {
1374 part = 0;
1375 layer = 2;
1376 mat = 0;
1377 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1378 b[0] = tb[0];
1379 b[1] = -tb[1];
1380 b[2] = tb[2];
1381 ifend = true;
1382 }
1383 if ( 2240 * mm <= mr[2] && mr[2] < 2280 * mm && 1133 * mm <= mr[0] &&
1384 mr[0] <= 2500 * mm )
1385 {
1386 part = 0;
1387 layer = 2;
1388 mat = 1;
1389 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1390 b[0] = tb[0];
1391 b[1] = -tb[1];
1392 b[2] = tb[2];
1393 ifend = true;
1394 }
1395 if ( 2280 * mm <= mr[2] && mr[2] < 2310 * mm && 1133 * mm <= mr[0] &&
1396 mr[0] <= 2500 * mm )
1397 {
1398 part = 0;
1399 layer = 3;
1400 mat = 0;
1401 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1402 b[0] = tb[0];
1403 b[1] = -tb[1];
1404 b[2] = tb[2];
1405 ifend = true;
1406 }
1407 if ( 2310 * mm <= mr[2] && mr[2] < 2350 * mm && 1167 * mm <= mr[0] &&
1408 mr[0] <= 2500 * mm )
1409 {
1410 part = 0;
1411 layer = 3;
1412 mat = 1;
1413 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1414 b[0] = tb[0];
1415 b[1] = -tb[1];
1416 b[2] = tb[2];
1417 ifend = true;
1418 }
1419 if ( 2350 * mm <= mr[2] && mr[2] < 2380 * mm && 1167 * mm <= mr[0] &&
1420 mr[0] <= 2500 * mm )
1421 {
1422 part = 0;
1423 layer = 4;
1424 mat = 0;
1425 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1426 b[0] = tb[0];
1427 b[1] = -tb[1];
1428 b[2] = tb[2];
1429 ifend = true;
1430 }
1431 if ( 2380 * mm <= mr[2] && mr[2] < 2420 * mm && 1203 * mm <= mr[0] &&
1432 mr[0] <= 2500 * mm )
1433 {
1434 part = 0;
1435 layer = 4;
1436 mat = 1;
1437 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1438 b[0] = tb[0];
1439 b[1] = -tb[1];
1440 b[2] = tb[2];
1441 ifend = true;
1442 }
1443 if ( 2420 * mm <= mr[2] && mr[2] < 2470 * mm && 1203 * mm <= mr[0] &&
1444 mr[0] <= 2500 * mm )
1445 {
1446 part = 0;
1447 layer = 5;
1448 mat = 0;
1449 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1450 b[0] = tb[0];
1451 b[1] = -tb[1];
1452 b[2] = tb[2];
1453 ifend = true;
1454 }
1455 if ( 2470 * mm <= mr[2] && mr[2] < 2510 * mm && 1241 * mm <= mr[0] &&
1456 mr[0] <= 2500 * mm )
1457 {
1458 part = 0;
1459 layer = 5;
1460 mat = 1;
1461 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1462 b[0] = tb[0];
1463 b[1] = -tb[1];
1464 b[2] = tb[2];
1465 ifend = true;
1466 }
1467 if ( 2510 * mm <= mr[2] && mr[2] < 2590 * mm && 1241 * mm <= mr[0] &&
1468 mr[0] <= 2500 * mm )
1469 {
1470 part = 0;
1471 layer = 6;
1472 mat = 0;
1473 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1474 b[0] = tb[0];
1475 b[1] = -tb[1];
1476 b[2] = tb[2];
1477 ifend = true;
1478 }
1479 if ( 2590 * mm <= mr[2] && mr[2] < 2630 * mm && 1302 * mm <= mr[0] &&
1480 mr[0] <= 2500 * mm )
1481 {
1482 part = 0;
1483 layer = 6;
1484 mat = 1;
1485 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1486 b[0] = tb[0];
1487 b[1] = -tb[1];
1488 b[2] = tb[2];
1489 ifend = true;
1490 }
1491 if ( 2630 * mm <= mr[2] && mr[2] < 2710 * mm && 1302 * mm <= mr[0] &&
1492 mr[0] <= 2500 * mm )
1493 {
1494 part = 0;
1495 layer = 7;
1496 mat = 0;
1497 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1498 b[0] = tb[0];
1499 b[1] = -tb[1];
1500 b[2] = tb[2];
1501 ifend = true;
1502 }
1503 if ( 2710 * mm <= mr[2] && mr[2] < 2750 * mm && 1362 * mm <= mr[0] &&
1504 mr[0] <= 2500 * mm )
1505 {
1506 part = 0;
1507 layer = 7;
1508 mat = 1;
1509 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1510 b[0] = tb[0];
1511 b[1] = -tb[1];
1512 b[2] = tb[2];
1513 ifend = true;
1514 }
1515 if ( 2750 * mm <= mr[2] && mr[2] <= 2800 * mm && 1302 * mm <= mr[0] &&
1516 mr[0] <= 2500 * mm )
1517 {
1518 part = 0;
1519 layer = 8;
1520 mat = 0;
1521 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1522 b[0] = tb[0];
1523 b[1] = -tb[1];
1524 b[2] = tb[2];
1525 ifend = true;
1526 }
1527 }
1528 if ( pi / 8 <= theta && theta < pi / 4 )
1529 {
1530 mr[0] = fabs( r.x() ) * cos( pi / 4 ) + fabs( r.y() ) * sin( pi / 4 );
1531 mr[1] = -fabs( r.x() ) * sin( pi / 4 ) + fabs( r.y() ) * cos( pi / 4 );
1532 mr[2] = fabs( r.z() );
1533 if ( mr[2] <= 1970 * mm && 1740 * mm <= mr[0] && mr[0] <= 2620 * mm )
1534 {
1535 part = 1;
1536 if ( 1740 * mm <= mr[0] && mr[0] < 1770 * mm )
1537 {
1538 layer = 0;
1539 mat = 0;
1540 }
1541 if ( 1770 * mm <= mr[0] && mr[0] < 1810 * mm )
1542 {
1543 layer = 0;
1544 mat = 1;
1545 }
1546 if ( 1810 * mm <= mr[0] && mr[0] < 1840 * mm )
1547 {
1548 layer = 1;
1549 mat = 0;
1550 }
1551 if ( 1840 * mm <= mr[0] && mr[0] < 1880 * mm )
1552 {
1553 layer = 1;
1554 mat = 1;
1555 }
1556 if ( 1880 * mm <= mr[0] && mr[0] < 1910 * mm )
1557 {
1558 layer = 2;
1559 mat = 0;
1560 }
1561 if ( 1910 * mm <= mr[0] && mr[0] < 1950 * mm )
1562 {
1563 layer = 2;
1564 mat = 1;
1565 }
1566 if ( 1950 * mm <= mr[0] && mr[0] < 1990 * mm )
1567 {
1568 layer = 3;
1569 mat = 0;
1570 }
1571 if ( 1990 * mm <= mr[0] && mr[0] < 2030 * mm )
1572 {
1573 layer = 3;
1574 mat = 1;
1575 }
1576 if ( 2030 * mm <= mr[0] && mr[0] < 2070 * mm )
1577 {
1578 layer = 4;
1579 mat = 0;
1580 }
1581 if ( 2070 * mm <= mr[0] && mr[0] < 2110 * mm )
1582 {
1583 layer = 4;
1584 mat = 1;
1585 }
1586 if ( 2110 * mm <= mr[0] && mr[0] < 2190 * mm )
1587 {
1588 layer = 5;
1589 mat = 0;
1590 }
1591 if ( 2190 * mm <= mr[0] && mr[0] < 2230 * mm )
1592 {
1593 layer = 5;
1594 mat = 1;
1595 }
1596 if ( 2230 * mm <= mr[0] && mr[0] < 2310 * mm )
1597 {
1598 layer = 6;
1599 mat = 0;
1600 }
1601 if ( 2310 * mm <= mr[0] && mr[0] < 2350 * mm )
1602 {
1603 layer = 6;
1604 mat = 1;
1605 }
1606 if ( 2350 * mm <= mr[0] && mr[0] < 2430 * mm )
1607 {
1608 layer = 7;
1609 mat = 0;
1610 }
1611 if ( 2430 * mm <= mr[0] && mr[0] < 2470 * mm )
1612 {
1613 layer = 7;
1614 mat = 1;
1615 }
1616 if ( 2470 * mm <= mr[0] && mr[0] <= 2620 * mm )
1617 {
1618 layer = 8;
1619 mat = 0;
1620 }
1621 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1622 b[0] = tb[0] * cos( pi / 4 ) - tb[1] * sin( pi / 4 );
1623 b[1] = tb[0] * sin( pi / 4 ) + tb[1] * cos( pi / 4 );
1624 b[2] = tb[2];
1625 ifbar = true;
1626 }
1627 if ( 2050 * mm <= mr[2] && mr[2] < 2090 * mm && 1034 * mm <= mr[0] &&
1628 mr[0] <= 2500 * mm )
1629 {
1630 part = 0;
1631 layer = 0;
1632 mat = 0;
1633 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1634 b[0] = tb[0] * cos( pi / 4 ) - tb[1] * sin( pi / 4 );
1635 b[1] = tb[0] * sin( pi / 4 ) + tb[1] * cos( pi / 4 );
1636 b[2] = tb[2];
1637 ifend = true;
1638 }
1639 if ( 2090 * mm <= mr[2] && mr[2] < 2130 * mm && 1067 * mm <= mr[0] &&
1640 mr[0] <= 2500 * mm )
1641 {
1642 part = 0;
1643 layer = 0;
1644 mat = 1;
1645 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1646 b[0] = tb[0] * cos( pi / 4 ) - tb[1] * sin( pi / 4 );
1647 b[1] = tb[0] * sin( pi / 4 ) + tb[1] * cos( pi / 4 );
1648 b[2] = tb[2];
1649 ifend = true;
1650 }
1651 if ( 2130 * mm <= mr[2] && mr[2] < 2170 * mm && 1067 * mm <= mr[0] &&
1652 mr[0] <= 2500 * mm )
1653 {
1654 part = 0;
1655 layer = 1;
1656 mat = 0;
1657 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1658 b[0] = tb[0] * cos( pi / 4 ) - tb[1] * sin( pi / 4 );
1659 b[1] = tb[0] * sin( pi / 4 ) + tb[1] * cos( pi / 4 );
1660 b[2] = tb[2];
1661 ifend = true;
1662 }
1663 if ( 2170 * mm <= mr[2] && mr[2] < 2210 * mm && 1100 * mm <= mr[0] &&
1664 mr[0] <= 2500 * mm )
1665 {
1666 part = 0;
1667 layer = 1;
1668 mat = 1;
1669 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1670 b[0] = tb[0] * cos( pi / 4 ) - tb[1] * sin( pi / 4 );
1671 b[1] = tb[0] * sin( pi / 4 ) + tb[1] * cos( pi / 4 );
1672 b[2] = tb[2];
1673 ifend = true;
1674 }
1675 if ( 2210 * mm <= mr[2] && mr[2] < 2240 * mm && 1100 * mm < mr[0] &&
1676 mr[0] <= 2500 * mm )
1677 {
1678 part = 0;
1679 layer = 2;
1680 mat = 0;
1681 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1682 b[0] = tb[0] * cos( pi / 4 ) - tb[1] * sin( pi / 4 );
1683 b[1] = tb[0] * sin( pi / 4 ) + tb[1] * cos( pi / 4 );
1684 b[2] = tb[2];
1685 ifend = true;
1686 }
1687 if ( 2240 * mm <= mr[2] && mr[2] < 2280 * mm && 1133 * mm <= mr[0] &&
1688 mr[0] <= 2500 * mm )
1689 {
1690 part = 0;
1691 layer = 2;
1692 mat = 1;
1693 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1694 b[0] = tb[0] * cos( pi / 4 ) - tb[1] * sin( pi / 4 );
1695 b[1] = tb[0] * sin( pi / 4 ) + tb[1] * cos( pi / 4 );
1696 b[2] = tb[2];
1697 ifend = true;
1698 }
1699 if ( 2280 * mm <= mr[2] && mr[2] < 2310 * mm && 1133 * mm <= mr[0] &&
1700 mr[0] <= 2500 * mm )
1701 {
1702 part = 0;
1703 layer = 3;
1704 mat = 0;
1705 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1706 b[0] = tb[0] * cos( pi / 4 ) - tb[1] * sin( pi / 4 );
1707 b[1] = tb[0] * sin( pi / 4 ) + tb[1] * cos( pi / 4 );
1708 b[2] = tb[2];
1709 ifend = true;
1710 }
1711 if ( 2310 * mm <= mr[2] && mr[2] < 2350 * mm && 1167 * mm <= mr[0] &&
1712 mr[0] <= 2500 * mm )
1713 {
1714 part = 0;
1715 layer = 3;
1716 mat = 1;
1717 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1718 b[0] = tb[0] * cos( pi / 4 ) - tb[1] * sin( pi / 4 );
1719 b[1] = tb[0] * sin( pi / 4 ) + tb[1] * cos( pi / 4 );
1720 b[2] = tb[2];
1721 ifend = true;
1722 }
1723 if ( 2350 * mm <= mr[2] && mr[2] < 2380 * mm && 1167 * mm <= mr[0] &&
1724 mr[0] <= 2500 * mm )
1725 {
1726 part = 0;
1727 layer = 4;
1728 mat = 0;
1729 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1730 b[0] = tb[0] * cos( pi / 4 ) - tb[1] * sin( pi / 4 );
1731 b[1] = tb[0] * sin( pi / 4 ) + tb[1] * cos( pi / 4 );
1732 b[2] = tb[2];
1733 ifend = true;
1734 }
1735 if ( 2380 * mm <= mr[2] && mr[2] < 2420 * mm && 1203 * mm <= mr[0] &&
1736 mr[0] <= 2500 * mm )
1737 {
1738 part = 0;
1739 layer = 4;
1740 mat = 1;
1741 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1742 b[0] = tb[0] * cos( pi / 4 ) - tb[1] * sin( pi / 4 );
1743 b[1] = tb[0] * sin( pi / 4 ) + tb[1] * cos( pi / 4 );
1744 b[2] = tb[2];
1745 ifend = true;
1746 }
1747 if ( 2420 * mm <= mr[2] && mr[2] < 2470 * mm && 1203 * mm <= mr[0] &&
1748 mr[0] <= 2500 * mm )
1749 {
1750 part = 0;
1751 layer = 5;
1752 mat = 0;
1753 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1754 b[0] = tb[0] * cos( pi / 4 ) - tb[1] * sin( pi / 4 );
1755 b[1] = tb[0] * sin( pi / 4 ) + tb[1] * cos( pi / 4 );
1756 b[2] = tb[2];
1757 ifend = true;
1758 }
1759 if ( 2470 * mm <= mr[2] && mr[2] < 2510 * mm && 1241 * mm <= mr[0] &&
1760 mr[0] <= 2500 * mm )
1761 {
1762 part = 0;
1763 layer = 5;
1764 mat = 1;
1765 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1766 b[0] = tb[0] * cos( pi / 4 ) - tb[1] * sin( pi / 4 );
1767 b[1] = tb[0] * sin( pi / 4 ) + tb[1] * cos( pi / 4 );
1768 b[2] = tb[2];
1769 ifend = true;
1770 }
1771 if ( 2510 * mm <= mr[2] && mr[2] < 2590 * mm && 1241 * mm <= mr[0] &&
1772 mr[0] <= 2500 * mm )
1773 {
1774 part = 0;
1775 layer = 6;
1776 mat = 0;
1777 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1778 b[0] = tb[0] * cos( pi / 4 ) - tb[1] * sin( pi / 4 );
1779 b[1] = tb[0] * sin( pi / 4 ) + tb[1] * cos( pi / 4 );
1780 b[2] = tb[2];
1781 ifend = true;
1782 }
1783 if ( 2590 * mm <= mr[2] && mr[2] < 2630 * mm && 1302 * mm <= mr[0] &&
1784 mr[0] <= 2500 * mm )
1785 {
1786 part = 0;
1787 layer = 6;
1788 mat = 1;
1789 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1790 b[0] = tb[0] * cos( pi / 4 ) - tb[1] * sin( pi / 4 );
1791 b[1] = tb[0] * sin( pi / 4 ) + tb[1] * cos( pi / 4 );
1792 b[2] = tb[2];
1793 ifend = true;
1794 }
1795 if ( 2630 * mm <= mr[2] && mr[2] < 2710 * mm && 1302 * mm <= mr[0] &&
1796 mr[0] <= 2500 * mm )
1797 {
1798 part = 0;
1799 layer = 7;
1800 mat = 0;
1801 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1802 b[0] = tb[0] * cos( pi / 4 ) - tb[1] * sin( pi / 4 );
1803 b[1] = tb[0] * sin( pi / 4 ) + tb[1] * cos( pi / 4 );
1804 b[2] = tb[2];
1805 ifend = true;
1806 }
1807 if ( 2710 * mm <= mr[2] && mr[2] < 2750 * mm && 1362 * mm <= mr[0] &&
1808 mr[0] <= 2500 * mm )
1809 {
1810 part = 0;
1811 layer = 7;
1812 mat = 1;
1813 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1814 b[0] = tb[0] * cos( pi / 4 ) - tb[1] * sin( pi / 4 );
1815 b[1] = tb[0] * sin( pi / 4 ) + tb[1] * cos( pi / 4 );
1816 b[2] = tb[2];
1817 ifend = true;
1818 }
1819 if ( 2750 * mm <= mr[2] && mr[2] <= 2800 * mm && 1302 * mm <= mr[0] &&
1820 mr[0] <= 2500 * mm )
1821 {
1822 part = 0;
1823 layer = 8;
1824 mat = 0;
1825 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1826 b[0] = tb[0] * cos( pi / 4 ) - tb[1] * sin( pi / 4 );
1827 b[1] = tb[0] * sin( pi / 4 ) + tb[1] * cos( pi / 4 );
1828 b[2] = tb[2];
1829 ifend = true;
1830 }
1831 }
1832 if ( pi / 4 <= theta && theta < 3 * pi / 8 )
1833 {
1834 mr[0] = fabs( r.x() ) * cos( pi / 4 ) + fabs( r.y() ) * sin( pi / 4 );
1835 mr[1] = fabs( r.x() ) * sin( pi / 4 ) - fabs( r.y() ) * cos( pi / 4 );
1836 mr[2] = fabs( r.z() );
1837 if ( mr[2] <= 1970 * mm && 1740 * mm <= mr[0] && mr[0] <= 2620 * mm )
1838 {
1839 part = 1;
1840 if ( 1740 * mm <= mr[0] && mr[0] < 1770 * mm )
1841 {
1842 layer = 0;
1843 mat = 0;
1844 }
1845 if ( 1770 * mm <= mr[0] && mr[0] < 1810 * mm )
1846 {
1847 layer = 0;
1848 mat = 1;
1849 }
1850 if ( 1810 * mm <= mr[0] && mr[0] < 1840 * mm )
1851 {
1852 layer = 1;
1853 mat = 0;
1854 }
1855 if ( 1840 * mm <= mr[0] && mr[0] < 1880 * mm )
1856 {
1857 layer = 1;
1858 mat = 1;
1859 }
1860 if ( 1880 * mm <= mr[0] && mr[0] < 1910 * mm )
1861 {
1862 layer = 2;
1863 mat = 0;
1864 }
1865 if ( 1910 * mm <= mr[0] && mr[0] < 1950 * mm )
1866 {
1867 layer = 2;
1868 mat = 1;
1869 }
1870 if ( 1950 * mm <= mr[0] && mr[0] < 1990 * mm )
1871 {
1872 layer = 3;
1873 mat = 0;
1874 }
1875 if ( 1990 * mm <= mr[0] && mr[0] < 2030 * mm )
1876 {
1877 layer = 3;
1878 mat = 1;
1879 }
1880 if ( 2030 * mm <= mr[0] && mr[0] < 2070 * mm )
1881 {
1882 layer = 4;
1883 mat = 0;
1884 }
1885 if ( 2070 * mm <= mr[0] && mr[0] < 2110 * mm )
1886 {
1887 layer = 4;
1888 mat = 1;
1889 }
1890 if ( 2110 * mm <= mr[0] && mr[0] < 2190 * mm )
1891 {
1892 layer = 5;
1893 mat = 0;
1894 }
1895 if ( 2190 * mm <= mr[0] && mr[0] < 2230 * mm )
1896 {
1897 layer = 5;
1898 mat = 1;
1899 }
1900 if ( 2230 * mm <= mr[0] && mr[0] < 2310 * mm )
1901 {
1902 layer = 6;
1903 mat = 0;
1904 }
1905 if ( 2310 * mm <= mr[0] && mr[0] < 2350 * mm )
1906 {
1907 layer = 6;
1908 mat = 1;
1909 }
1910 if ( 2350 * mm <= mr[0] && mr[0] < 2430 * mm )
1911 {
1912 layer = 7;
1913 mat = 0;
1914 }
1915 if ( 2430 * mm <= mr[0] && mr[0] < 2470 * mm )
1916 {
1917 layer = 7;
1918 mat = 1;
1919 }
1920 if ( 2470 * mm <= mr[0] && mr[0] <= 2620 * mm )
1921 {
1922 layer = 8;
1923 mat = 0;
1924 }
1925 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1926 b[0] = tb[0] * cos( pi / 4 ) + tb[1] * sin( pi / 4 );
1927 b[1] = tb[0] * sin( pi / 4 ) - tb[1] * cos( pi / 4 );
1928 b[2] = tb[2];
1929 ifbar = true;
1930 }
1931 if ( 2050 * mm <= mr[2] && mr[2] < 2090 * mm && 1034 * mm <= mr[0] &&
1932 mr[0] <= 2500 * mm )
1933 {
1934 part = 0;
1935 layer = 0;
1936 mat = 0;
1937 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1938 b[0] = tb[0] * cos( pi / 4 ) + tb[1] * sin( pi / 4 );
1939 b[1] = tb[0] * sin( pi / 4 ) - tb[1] * cos( pi / 4 );
1940 b[2] = tb[2];
1941 ifend = true;
1942 }
1943 if ( 2090 * mm <= mr[2] && mr[2] < 2130 * mm && 1067 * mm <= mr[0] &&
1944 mr[0] <= 2500 * mm )
1945 {
1946 part = 0;
1947 layer = 0;
1948 mat = 1;
1949 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1950 b[0] = tb[0] * cos( pi / 4 ) + tb[1] * sin( pi / 4 );
1951 b[1] = tb[0] * sin( pi / 4 ) - tb[1] * cos( pi / 4 );
1952 b[2] = tb[2];
1953 ifend = true;
1954 }
1955 if ( 2130 * mm <= mr[2] && mr[2] < 2170 * mm && 1067 * mm <= mr[0] &&
1956 mr[0] <= 2500 * mm )
1957 {
1958 part = 0;
1959 layer = 1;
1960 mat = 0;
1961 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1962 b[0] = tb[0] * cos( pi / 4 ) + tb[1] * sin( pi / 4 );
1963 b[1] = tb[0] * sin( pi / 4 ) - tb[1] * cos( pi / 4 );
1964 b[2] = tb[2];
1965 ifend = true;
1966 }
1967 if ( 2170 * mm <= mr[2] && mr[2] < 2210 * mm && 1100 * mm <= mr[0] &&
1968 mr[0] <= 2500 * mm )
1969 {
1970 part = 0;
1971 layer = 1;
1972 mat = 1;
1973 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1974 b[0] = tb[0] * cos( pi / 4 ) + tb[1] * sin( pi / 4 );
1975 b[1] = tb[0] * sin( pi / 4 ) - tb[1] * cos( pi / 4 );
1976 b[2] = tb[2];
1977 ifend = true;
1978 }
1979 if ( 2210 * mm <= mr[2] && mr[2] < 2240 * mm && 1100 * mm < mr[0] &&
1980 mr[0] <= 2500 * mm )
1981 {
1982 part = 0;
1983 layer = 2;
1984 mat = 0;
1985 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1986 b[0] = tb[0] * cos( pi / 4 ) + tb[1] * sin( pi / 4 );
1987 b[1] = tb[0] * sin( pi / 4 ) - tb[1] * cos( pi / 4 );
1988 b[2] = tb[2];
1989 ifend = true;
1990 }
1991 if ( 2240 * mm <= mr[2] && mr[2] < 2280 * mm && 1133 * mm <= mr[0] &&
1992 mr[0] <= 2500 * mm )
1993 {
1994 part = 0;
1995 layer = 2;
1996 mat = 1;
1997 m_Mucfield->getMucField( part, layer, mat, mr, tb );
1998 b[0] = tb[0] * cos( pi / 4 ) + tb[1] * sin( pi / 4 );
1999 b[1] = tb[0] * sin( pi / 4 ) - tb[1] * cos( pi / 4 );
2000 b[2] = tb[2];
2001 ifend = true;
2002 }
2003 if ( 2280 * mm <= mr[2] && mr[2] < 2310 * mm && 1133 * mm <= mr[0] &&
2004 mr[0] <= 2500 * mm )
2005 {
2006 part = 0;
2007 layer = 3;
2008 mat = 0;
2009 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2010 b[0] = tb[0] * cos( pi / 4 ) + tb[1] * sin( pi / 4 );
2011 b[1] = tb[0] * sin( pi / 4 ) - tb[1] * cos( pi / 4 );
2012 b[2] = tb[2];
2013 ifend = true;
2014 }
2015 if ( 2310 * mm <= mr[2] && mr[2] < 2350 * mm && 1167 * mm <= mr[0] &&
2016 mr[0] <= 2500 * mm )
2017 {
2018 part = 0;
2019 layer = 3;
2020 mat = 1;
2021 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2022 b[0] = tb[0] * cos( pi / 4 ) + tb[1] * sin( pi / 4 );
2023 b[1] = tb[0] * sin( pi / 4 ) - tb[1] * cos( pi / 4 );
2024 b[2] = tb[2];
2025 ifend = true;
2026 }
2027 if ( 2350 * mm <= mr[2] && mr[2] < 2380 * mm && 1167 * mm <= mr[0] &&
2028 mr[0] <= 2500 * mm )
2029 {
2030 part = 0;
2031 layer = 4;
2032 mat = 0;
2033 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2034 b[0] = tb[0] * cos( pi / 4 ) + tb[1] * sin( pi / 4 );
2035 b[1] = tb[0] * sin( pi / 4 ) - tb[1] * cos( pi / 4 );
2036 b[2] = tb[2];
2037 ifend = true;
2038 }
2039 if ( 2380 * mm <= mr[2] && mr[2] < 2420 * mm && 1203 * mm <= mr[0] &&
2040 mr[0] <= 2500 * mm )
2041 {
2042 part = 0;
2043 layer = 4;
2044 mat = 1;
2045 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2046 b[0] = tb[0] * cos( pi / 4 ) + tb[1] * sin( pi / 4 );
2047 b[1] = tb[0] * sin( pi / 4 ) - tb[1] * cos( pi / 4 );
2048 b[2] = tb[2];
2049 ifend = true;
2050 }
2051 if ( 2420 * mm <= mr[2] && mr[2] < 2470 * mm && 1203 * mm <= mr[0] &&
2052 mr[0] <= 2500 * mm )
2053 {
2054 part = 0;
2055 layer = 5;
2056 mat = 0;
2057 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2058 b[0] = tb[0] * cos( pi / 4 ) + tb[1] * sin( pi / 4 );
2059 b[1] = tb[0] * sin( pi / 4 ) - tb[1] * cos( pi / 4 );
2060 b[2] = tb[2];
2061 ifend = true;
2062 }
2063 if ( 2470 * mm <= mr[2] && mr[2] < 2510 * mm && 1241 * mm <= mr[0] &&
2064 mr[0] <= 2500 * mm )
2065 {
2066 part = 0;
2067 layer = 5;
2068 mat = 1;
2069 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2070 b[0] = tb[0] * cos( pi / 4 ) + tb[1] * sin( pi / 4 );
2071 b[1] = tb[0] * sin( pi / 4 ) - tb[1] * cos( pi / 4 );
2072 b[2] = tb[2];
2073 ifend = true;
2074 }
2075 if ( 2510 * mm <= mr[2] && mr[2] < 2590 * mm && 1241 * mm <= mr[0] &&
2076 mr[0] <= 2500 * mm )
2077 {
2078 part = 0;
2079 layer = 6;
2080 mat = 0;
2081 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2082 b[0] = tb[0] * cos( pi / 4 ) + tb[1] * sin( pi / 4 );
2083 b[1] = tb[0] * sin( pi / 4 ) - tb[1] * cos( pi / 4 );
2084 b[2] = tb[2];
2085 ifend = true;
2086 }
2087 if ( 2590 * mm <= mr[2] && mr[2] < 2630 * mm && 1302 * mm <= mr[0] &&
2088 mr[0] <= 2500 * mm )
2089 {
2090 part = 0;
2091 layer = 6;
2092 mat = 1;
2093 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2094 b[0] = tb[0] * cos( pi / 4 ) + tb[1] * sin( pi / 4 );
2095 b[1] = tb[0] * sin( pi / 4 ) - tb[1] * cos( pi / 4 );
2096 b[2] = tb[2];
2097 ifend = true;
2098 }
2099 if ( 2630 * mm <= mr[2] && mr[2] < 2710 * mm && 1302 * mm <= mr[0] &&
2100 mr[0] <= 2500 * mm )
2101 {
2102 part = 0;
2103 layer = 7;
2104 mat = 0;
2105 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2106 b[0] = tb[0] * cos( pi / 4 ) + tb[1] * sin( pi / 4 );
2107 b[1] = tb[0] * sin( pi / 4 ) - tb[1] * cos( pi / 4 );
2108 b[2] = tb[2];
2109 ifend = true;
2110 }
2111 if ( 2710 * mm <= mr[2] && mr[2] < 2750 * mm && 1362 * mm <= mr[0] &&
2112 mr[0] <= 2500 * mm )
2113 {
2114 part = 0;
2115 layer = 7;
2116 mat = 1;
2117 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2118 b[0] = tb[0] * cos( pi / 4 ) + tb[1] * sin( pi / 4 );
2119 b[1] = tb[0] * sin( pi / 4 ) - tb[1] * cos( pi / 4 );
2120 b[2] = tb[2];
2121 ifend = true;
2122 }
2123 if ( 2750 * mm <= mr[2] && mr[2] <= 2800 * mm && 1302 * mm <= mr[0] &&
2124 mr[0] <= 2500 * mm )
2125 {
2126 part = 0;
2127 layer = 8;
2128 mat = 0;
2129 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2130 b[0] = tb[0] * cos( pi / 4 ) + tb[1] * sin( pi / 4 );
2131 b[1] = tb[0] * sin( pi / 4 ) - tb[1] * cos( pi / 4 );
2132 b[2] = tb[2];
2133 ifend = true;
2134 }
2135 }
2136 if ( 3 * pi / 8 <= theta && theta < pi / 2 )
2137 {
2138 mr[0] = fabs( r.y() );
2139 mr[1] = -fabs( r.x() );
2140 mr[2] = fabs( r.z() );
2141 if ( mr[2] <= 1970 * mm && 1740 * mm <= mr[0] && mr[0] <= 2620 * mm )
2142 {
2143 part = 1;
2144 if ( 1740 * mm <= mr[0] && mr[0] < 1770 * mm )
2145 {
2146 layer = 0;
2147 mat = 0;
2148 }
2149 if ( 1770 * mm <= mr[0] && mr[0] < 1810 * mm )
2150 {
2151 layer = 0;
2152 mat = 1;
2153 }
2154 if ( 1810 * mm <= mr[0] && mr[0] < 1840 * mm )
2155 {
2156 layer = 1;
2157 mat = 0;
2158 }
2159 if ( 1840 * mm <= mr[0] && mr[0] < 1880 * mm )
2160 {
2161 layer = 1;
2162 mat = 1;
2163 }
2164 if ( 1880 * mm <= mr[0] && mr[0] < 1910 * mm )
2165 {
2166 layer = 2;
2167 mat = 0;
2168 }
2169 if ( 1910 * mm <= mr[0] && mr[0] < 1950 * mm )
2170 {
2171 layer = 2;
2172 mat = 1;
2173 }
2174 if ( 1950 * mm <= mr[0] && mr[0] < 1990 * mm )
2175 {
2176 layer = 3;
2177 mat = 0;
2178 }
2179 if ( 1990 * mm <= mr[0] && mr[0] < 2030 * mm )
2180 {
2181 layer = 3;
2182 mat = 1;
2183 }
2184 if ( 2030 * mm <= mr[0] && mr[0] < 2070 * mm )
2185 {
2186 layer = 4;
2187 mat = 0;
2188 }
2189 if ( 2070 * mm <= mr[0] && mr[0] < 2110 * mm )
2190 {
2191 layer = 4;
2192 mat = 1;
2193 }
2194 if ( 2110 * mm <= mr[0] && mr[0] < 2190 * mm )
2195 {
2196 layer = 5;
2197 mat = 0;
2198 }
2199 if ( 2190 * mm <= mr[0] && mr[0] < 2230 * mm )
2200 {
2201 layer = 5;
2202 mat = 1;
2203 }
2204 if ( 2230 * mm <= mr[0] && mr[0] < 2310 * mm )
2205 {
2206 layer = 6;
2207 mat = 0;
2208 }
2209 if ( 2310 * mm <= mr[0] && mr[0] < 2350 * mm )
2210 {
2211 layer = 6;
2212 mat = 1;
2213 }
2214 if ( 2350 * mm <= mr[0] && mr[0] < 2430 * mm )
2215 {
2216 layer = 7;
2217 mat = 0;
2218 }
2219 if ( 2430 * mm <= mr[0] && mr[0] < 2470 * mm )
2220 {
2221 layer = 7;
2222 mat = 1;
2223 }
2224 if ( 2470 * mm <= mr[0] && mr[0] <= 2620 * mm )
2225 {
2226 layer = 8;
2227 mat = 0;
2228 }
2229 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2230 b[0] = -tb[1];
2231 b[1] = tb[0];
2232 b[2] = tb[2];
2233 ifbar = true;
2234 }
2235 if ( 2050 * mm <= mr[2] && mr[2] < 2090 * mm && 1034 * mm <= mr[0] &&
2236 mr[0] <= 2500 * mm )
2237 {
2238 part = 0;
2239 layer = 0;
2240 mat = 0;
2241 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2242 b[0] = -tb[1];
2243 b[1] = tb[0];
2244 b[2] = tb[2];
2245 ifend = true;
2246 }
2247 if ( 2090 * mm <= mr[2] && mr[2] < 2130 * mm && 1067 * mm <= mr[0] &&
2248 mr[0] <= 2500 * mm )
2249 {
2250 part = 0;
2251 layer = 0;
2252 mat = 1;
2253 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2254 b[0] = -tb[1];
2255 b[1] = tb[0];
2256 b[2] = tb[2];
2257 ifend = true;
2258 }
2259 if ( 2130 * mm <= mr[2] && mr[2] < 2170 * mm && 1067 * mm <= mr[0] &&
2260 mr[0] <= 2500 * mm )
2261 {
2262 part = 0;
2263 layer = 1;
2264 mat = 0;
2265 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2266 b[0] = -tb[1];
2267 b[1] = tb[0];
2268 b[2] = tb[2];
2269 ifend = true;
2270 }
2271 if ( 2170 * mm <= mr[2] && mr[2] < 2210 * mm && 1100 * mm <= mr[0] &&
2272 mr[0] <= 2500 * mm )
2273 {
2274 part = 0;
2275 layer = 1;
2276 mat = 1;
2277 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2278 b[0] = -tb[1];
2279 b[1] = tb[0];
2280 b[2] = tb[2];
2281 ifend = true;
2282 }
2283 if ( 2210 * mm <= mr[2] && mr[2] < 2240 * mm && 1100 * mm < mr[0] &&
2284 mr[0] <= 2500 * mm )
2285 {
2286 part = 0;
2287 layer = 2;
2288 mat = 0;
2289 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2290 b[0] = -tb[1];
2291 b[1] = tb[0];
2292 b[2] = tb[2];
2293 ifend = true;
2294 }
2295 if ( 2240 * mm <= mr[2] && mr[2] < 2280 * mm && 1133 * mm <= mr[0] &&
2296 mr[0] <= 2500 * mm )
2297 {
2298 part = 0;
2299 layer = 2;
2300 mat = 1;
2301 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2302 b[0] = -tb[1];
2303 b[1] = tb[0];
2304 b[2] = tb[2];
2305 ifend = true;
2306 }
2307 if ( 2280 * mm <= mr[2] && mr[2] < 2310 * mm && 1133 * mm <= mr[0] &&
2308 mr[0] <= 2500 * mm )
2309 {
2310 part = 0;
2311 layer = 3;
2312 mat = 0;
2313 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2314 b[0] = -tb[1];
2315 b[1] = tb[0];
2316 b[2] = tb[2];
2317 ifend = true;
2318 }
2319 if ( 2310 * mm <= mr[2] && mr[2] < 2350 * mm && 1167 * mm <= mr[0] &&
2320 mr[0] <= 2500 * mm )
2321 {
2322 part = 0;
2323 layer = 3;
2324 mat = 1;
2325 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2326 b[0] = -tb[1];
2327 b[1] = tb[0];
2328 b[2] = tb[2];
2329 ifend = true;
2330 }
2331 if ( 2350 * mm <= mr[2] && mr[2] < 2380 * mm && 1167 * mm <= mr[0] &&
2332 mr[0] <= 2500 * mm )
2333 {
2334 part = 0;
2335 layer = 4;
2336 mat = 0;
2337 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2338 b[0] = -tb[1];
2339 b[1] = tb[0];
2340 b[2] = tb[2];
2341 ifend = true;
2342 }
2343 if ( 2380 * mm <= mr[2] && mr[2] < 2420 * mm && 1203 * mm <= mr[0] &&
2344 mr[0] <= 2500 * mm )
2345 {
2346 part = 0;
2347 layer = 4;
2348 mat = 1;
2349 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2350 b[0] = -tb[1];
2351 b[1] = tb[0];
2352 b[2] = tb[2];
2353 ifend = true;
2354 }
2355 if ( 2420 * mm <= mr[2] && mr[2] < 2470 * mm && 1203 * mm <= mr[0] &&
2356 mr[0] <= 2500 * mm )
2357 {
2358 part = 0;
2359 layer = 5;
2360 mat = 0;
2361 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2362 b[0] = -tb[1];
2363 b[1] = tb[0];
2364 b[2] = tb[2];
2365 ifend = true;
2366 }
2367 if ( 2470 * mm <= mr[2] && mr[2] < 2510 * mm && 1241 * mm <= mr[0] &&
2368 mr[0] <= 2500 * mm )
2369 {
2370 part = 0;
2371 layer = 5;
2372 mat = 1;
2373 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2374 b[0] = -tb[1];
2375 b[1] = tb[0];
2376 b[2] = tb[2];
2377 ifend = true;
2378 }
2379 if ( 2510 * mm <= mr[2] && mr[2] < 2590 * mm && 1241 * mm <= mr[0] &&
2380 mr[0] <= 2500 * mm )
2381 {
2382 part = 0;
2383 layer = 6;
2384 mat = 0;
2385 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2386 b[0] = -tb[1];
2387 b[1] = tb[0];
2388 b[2] = tb[2];
2389 ifend = true;
2390 }
2391 if ( 2590 * mm <= mr[2] && mr[2] < 2630 * mm && 1302 * mm <= mr[0] &&
2392 mr[0] <= 2500 * mm )
2393 {
2394 part = 0;
2395 layer = 6;
2396 mat = 1;
2397 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2398 b[0] = -tb[1];
2399 b[1] = tb[0];
2400 b[2] = tb[2];
2401 ifend = true;
2402 }
2403 if ( 2630 * mm <= mr[2] && mr[2] < 2710 * mm && 1302 * mm <= mr[0] &&
2404 mr[0] <= 2500 * mm )
2405 {
2406 part = 0;
2407 layer = 7;
2408 mat = 0;
2409 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2410 b[0] = -tb[1];
2411 b[1] = tb[0];
2412 b[2] = tb[2];
2413 ifend = true;
2414 }
2415 if ( 2710 * mm <= mr[2] && mr[2] < 2750 * mm && 1362 * mm <= mr[0] &&
2416 mr[0] <= 2500 * mm )
2417 {
2418 part = 0;
2419 layer = 7;
2420 mat = 1;
2421 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2422 b[0] = -tb[1];
2423 b[1] = tb[0];
2424 b[2] = tb[2];
2425 ifend = true;
2426 }
2427 if ( 2750 * mm <= mr[2] && mr[2] <= 2800 * mm && 1302 * mm <= mr[0] &&
2428 mr[0] <= 2500 * mm )
2429 {
2430 part = 0;
2431 layer = 8;
2432 mat = 0;
2433 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2434 b[0] = -tb[1];
2435 b[1] = tb[0];
2436 b[2] = tb[2];
2437 ifend = true;
2438 }
2439 }
2440 }
2441 if ( r.x() == 0. )
2442 {
2443 mr[0] = fabs( r.y() );
2444 mr[1] = -fabs( r.x() );
2445 mr[2] = fabs( r.z() );
2446 if ( mr[2] <= 1970 * mm && 1740 * mm <= mr[0] && mr[0] <= 2620 * mm )
2447 {
2448 part = 1;
2449 if ( 1740 * mm <= mr[0] && mr[0] < 1770 * mm )
2450 {
2451 layer = 0;
2452 mat = 0;
2453 }
2454 if ( 1770 * mm <= mr[0] && mr[0] < 1810 * mm )
2455 {
2456 layer = 0;
2457 mat = 1;
2458 }
2459 if ( 1810 * mm <= mr[0] && mr[0] < 1840 * mm )
2460 {
2461 layer = 1;
2462 mat = 0;
2463 }
2464 if ( 1840 * mm <= mr[0] && mr[0] < 1880 * mm )
2465 {
2466 layer = 1;
2467 mat = 1;
2468 }
2469 if ( 1880 * mm <= mr[0] && mr[0] < 1910 * mm )
2470 {
2471 layer = 2;
2472 mat = 0;
2473 }
2474 if ( 1910 * mm <= mr[0] && mr[0] < 1950 * mm )
2475 {
2476 layer = 2;
2477 mat = 1;
2478 }
2479 if ( 1950 * mm <= mr[0] && mr[0] < 1990 * mm )
2480 {
2481 layer = 3;
2482 mat = 0;
2483 }
2484 if ( 1990 * mm <= mr[0] && mr[0] < 2030 * mm )
2485 {
2486 layer = 3;
2487 mat = 1;
2488 }
2489 if ( 2030 * mm <= mr[0] && mr[0] < 2070 * mm )
2490 {
2491 layer = 4;
2492 mat = 0;
2493 }
2494 if ( 2070 * mm <= mr[0] && mr[0] < 2110 * mm )
2495 {
2496 layer = 4;
2497 mat = 1;
2498 }
2499 if ( 2110 * mm <= mr[0] && mr[0] < 2190 * mm )
2500 {
2501 layer = 5;
2502 mat = 0;
2503 }
2504 if ( 2190 * mm <= mr[0] && mr[0] < 2230 * mm )
2505 {
2506 layer = 5;
2507 mat = 1;
2508 }
2509 if ( 2230 * mm <= mr[0] && mr[0] < 2310 * mm )
2510 {
2511 layer = 6;
2512 mat = 0;
2513 }
2514 if ( 2310 * mm <= mr[0] && mr[0] < 2350 * mm )
2515 {
2516 layer = 6;
2517 mat = 1;
2518 }
2519 if ( 2350 * mm <= mr[0] && mr[0] < 2430 * mm )
2520 {
2521 layer = 7;
2522 mat = 0;
2523 }
2524 if ( 2430 * mm <= mr[0] && mr[0] < 2470 * mm )
2525 {
2526 layer = 7;
2527 mat = 1;
2528 }
2529 if ( 2470 * mm <= mr[0] && mr[0] <= 2620 * mm )
2530 {
2531 layer = 8;
2532 mat = 0;
2533 }
2534 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2535 b[0] = -tb[1];
2536 b[1] = tb[0];
2537 b[2] = tb[2];
2538 ifbar = true;
2539 }
2540 if ( 2050 * mm <= mr[2] && mr[2] < 2090 * mm && 1034 * mm <= mr[0] &&
2541 mr[0] <= 2500 * mm )
2542 {
2543 part = 0;
2544 layer = 0;
2545 mat = 0;
2546 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2547 b[0] = -tb[1];
2548 b[1] = tb[0];
2549 b[2] = tb[2];
2550 ifend = true;
2551 }
2552 if ( 2090 * mm <= mr[2] && mr[2] < 2130 * mm && 1067 * mm <= mr[0] &&
2553 mr[0] <= 2500 * mm )
2554 {
2555 part = 0;
2556 layer = 0;
2557 mat = 1;
2558 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2559 b[0] = -tb[1];
2560 b[1] = tb[0];
2561 b[2] = tb[2];
2562 ifend = true;
2563 }
2564 if ( 2130 * mm <= mr[2] && mr[2] < 2170 * mm && 1067 * mm <= mr[0] &&
2565 mr[0] <= 2500 * mm )
2566 {
2567 part = 0;
2568 layer = 1;
2569 mat = 0;
2570 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2571 b[0] = -tb[1];
2572 b[1] = tb[0];
2573 b[2] = tb[2];
2574 ifend = true;
2575 }
2576 if ( 2170 * mm <= mr[2] && mr[2] < 2210 * mm && 1100 * mm <= mr[0] &&
2577 mr[0] <= 2500 * mm )
2578 {
2579 part = 0;
2580 layer = 1;
2581 mat = 1;
2582 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2583 b[0] = -tb[1];
2584 b[1] = tb[0];
2585 b[2] = tb[2];
2586 ifend = true;
2587 }
2588 if ( 2210 * mm <= mr[2] && mr[2] < 2240 * mm && 1100 * mm < mr[0] && mr[0] <= 2500 * mm )
2589 {
2590 part = 0;
2591 layer = 2;
2592 mat = 0;
2593 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2594 b[0] = -tb[1];
2595 b[1] = tb[0];
2596 b[2] = tb[2];
2597 ifend = true;
2598 }
2599 if ( 2240 * mm <= mr[2] && mr[2] < 2280 * mm && 1133 * mm <= mr[0] &&
2600 mr[0] <= 2500 * mm )
2601 {
2602 part = 0;
2603 layer = 2;
2604 mat = 1;
2605 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2606 b[0] = -tb[1];
2607 b[1] = tb[0];
2608 b[2] = tb[2];
2609 ifend = true;
2610 }
2611 if ( 2280 * mm <= mr[2] && mr[2] < 2310 * mm && 1133 * mm <= mr[0] &&
2612 mr[0] <= 2500 * mm )
2613 {
2614 part = 0;
2615 layer = 3;
2616 mat = 0;
2617 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2618 b[0] = -tb[1];
2619 b[1] = tb[0];
2620 b[2] = tb[2];
2621 ifend = true;
2622 }
2623 if ( 2310 * mm <= mr[2] && mr[2] < 2350 * mm && 1167 * mm <= mr[0] &&
2624 mr[0] <= 2500 * mm )
2625 {
2626 part = 0;
2627 layer = 3;
2628 mat = 1;
2629 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2630 b[0] = -tb[1];
2631 b[1] = tb[0];
2632 b[2] = tb[2];
2633 ifend = true;
2634 }
2635 if ( 2350 * mm <= mr[2] && mr[2] < 2380 * mm && 1167 * mm <= mr[0] &&
2636 mr[0] <= 2500 * mm )
2637 {
2638 part = 0;
2639 layer = 4;
2640 mat = 0;
2641 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2642 b[0] = -tb[1];
2643 b[1] = tb[0];
2644 b[2] = tb[2];
2645 ifend = true;
2646 }
2647 if ( 2380 * mm <= mr[2] && mr[2] < 2420 * mm && 1203 * mm <= mr[0] &&
2648 mr[0] <= 2500 * mm )
2649 {
2650 part = 0;
2651 layer = 4;
2652 mat = 1;
2653 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2654 b[0] = -tb[1];
2655 b[1] = tb[0];
2656 b[2] = tb[2];
2657 ifend = true;
2658 }
2659 if ( 2420 * mm <= mr[2] && mr[2] < 2470 * mm && 1203 * mm <= mr[0] &&
2660 mr[0] <= 2500 * mm )
2661 {
2662 part = 0;
2663 layer = 5;
2664 mat = 0;
2665 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2666 b[0] = -tb[1];
2667 b[1] = tb[0];
2668 b[2] = tb[2];
2669 ifend = true;
2670 }
2671 if ( 2470 * mm <= mr[2] && mr[2] < 2510 * mm && 1241 * mm <= mr[0] &&
2672 mr[0] <= 2500 * mm )
2673 {
2674 part = 0;
2675 layer = 5;
2676 mat = 1;
2677 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2678 b[0] = -tb[1];
2679 b[1] = tb[0];
2680 b[2] = tb[2];
2681 ifend = true;
2682 }
2683 if ( 2510 * mm <= mr[2] && mr[2] < 2590 * mm && 1241 * mm <= mr[0] &&
2684 mr[0] <= 2500 * mm )
2685 {
2686 part = 0;
2687 layer = 6;
2688 mat = 0;
2689 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2690 b[0] = -tb[1];
2691 b[1] = tb[0];
2692 b[2] = tb[2];
2693 ifend = true;
2694 }
2695 if ( 2590 * mm <= mr[2] && mr[2] < 2630 * mm && 1302 * mm <= mr[0] &&
2696 mr[0] <= 2500 * mm )
2697 {
2698 part = 0;
2699 layer = 6;
2700 mat = 1;
2701 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2702 b[0] = -tb[1];
2703 b[1] = tb[0];
2704 b[2] = tb[2];
2705 ifend = true;
2706 }
2707 if ( 2630 * mm <= mr[2] && mr[2] < 2710 * mm && 1302 * mm <= mr[0] &&
2708 mr[0] <= 2500 * mm )
2709 {
2710 part = 0;
2711 layer = 7;
2712 mat = 0;
2713 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2714 b[0] = -tb[1];
2715 b[1] = tb[0];
2716 b[2] = tb[2];
2717 ifend = true;
2718 }
2719 if ( 2710 * mm <= mr[2] && mr[2] < 2750 * mm && 1362 * mm <= mr[0] &&
2720 mr[0] <= 2500 * mm )
2721 {
2722 part = 0;
2723 layer = 7;
2724 mat = 1;
2725 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2726 b[0] = -tb[1];
2727 b[1] = tb[0];
2728 b[2] = tb[2];
2729 ifend = true;
2730 }
2731 if ( 2750 * mm <= mr[2] && mr[2] <= 2800 * mm && 1302 * mm <= mr[0] &&
2732 mr[0] <= 2500 * mm )
2733 {
2734 part = 0;
2735 layer = 8;
2736 mat = 0;
2737 m_Mucfield->getMucField( part, layer, mat, mr, tb );
2738 b[0] = -tb[1];
2739 b[1] = tb[0];
2740 b[2] = tb[2];
2741 ifend = true;
2742 }
2743 }
2744 if ( ifbar == true || ifend == true )
2745 {
2746 if ( r.x() < 0. && r.y() >= 0. && r.z() > 0. ) { b[0] = -b[0]; }
2747 else if ( r.x() <= 0. && r.y() < 0. && r.z() > 0. )
2748 {
2749 b[0] = -b[0];
2750 b[1] = -b[1];
2751 }
2752 else if ( r.x() > 0. && r.y() < 0. && r.z() > 0. ) { b[1] = -b[1]; }
2753 else if ( r.x() >= 0. && r.y() > 0. && r.z() < 0. )
2754 {
2755 b[0] = -b[0];
2756 b[1] = -b[1];
2757 }
2758 else if ( r.x() < 0. && r.y() >= 0. && r.z() < 0. ) { b[1] = -b[1]; }
2759 else if ( r.x() <= 0. && r.y() < 0. && r.z() < 0. )
2760 {
2761 b[0] = b[0];
2762 b[1] = b[1];
2763 }
2764 else if ( r.x() > 0. && r.y() <= 0. && r.z() < 0. ) { b[0] = -b[0]; }
2765 }
2766 }
2767
2768 // reference frames defined by simulation and SCM are different.
2769 newb[0] = -b[0] * m_scale;
2770 newb[1] = b[1] * m_scale;
2771 newb[2] = -b[2] * m_scale;
2772 /* 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)
2773 { return StatusCode::SUCCESS;
2774 }
2775 else {
2776 newb[0] = newb[0] - 0.10*tesla;
2777 newb[1] = newb[1] + 0.10*tesla;
2778 newb[2] = newb[2] - 0.10*tesla;
2779 }*/
2780
2781 // newb[0] = b[0];
2782 // newb[1] = b[1];
2783 // newb[2] = b[2];
2784
2785#ifndef BEAN
2786 return StatusCode::SUCCESS;
2787#else
2788 return true;
2789#endif
2790}
2791
2792StatusCode MagneticFieldSvc::uniFieldVector( const HepPoint3D& r, HepVector3D& b ) const {
2793 if ( m_runmode == 8 || m_runmode == 7 )
2794 {
2795 if ( -1.5 * m <= r.z() && r.z() <= 1.5 * m &&
2796 0 * m <= std::sqrt( r.x() * r.x() + r.y() * r.y() ) &&
2797 std::sqrt( r.x() * r.x() + r.y() * r.y() ) <= 1.5 * m )
2798 {
2799 b[0] = 0.0;
2800 b[1] = 0.0;
2801 b[2] = -0.9 * tesla;
2802 }
2803 else
2804 {
2805 b[0] = 0.0;
2806 b[1] = 0.0;
2807 b[2] = 0.0;
2808 }
2809 }
2810 else
2811 {
2812 if ( -1.5 * m <= r.z() && r.z() <= 1.5 * m &&
2813 0 * m <= std::sqrt( r.x() * r.x() + r.y() * r.y() ) &&
2814 std::sqrt( r.x() * r.x() + r.y() * r.y() ) <= 1.5 * m )
2815 {
2816 b[0] = 0.0;
2817 b[1] = 0.0;
2818 b[2] = -1.0 * tesla;
2819 }
2820 else
2821 {
2822 b[0] = 0.0;
2823 b[1] = 0.0;
2824 b[2] = 0.0;
2825 }
2826 }
2827
2828 if ( m_turnOffField == true )
2829 {
2830 b[0] = 0.;
2831 b[1] = 0.;
2832 b[2] = 0.;
2833 }
2834 // yzhang add 2012-04-25
2835 b[0] *= m_scale;
2836 b[1] *= m_scale;
2837 b[2] *= m_scale;
2838#ifndef BEAN
2839 return StatusCode::SUCCESS;
2840#else
2841 return true;
2842#endif
2843}
2844
2845double MagneticFieldSvc::getReferField() {
2846 if ( m_useDB == false )
2847 {
2848 if ( m_runmode == 8 || m_runmode == 7 ) m_zfield = -0.9 * tesla;
2849 else m_zfield = -1.0 * tesla;
2850 }
2851
2852 if ( m_turnOffField == true ) { m_zfield = 0.; }
2853 return m_zfield * m_scale;
2854}
2855
2856bool MagneticFieldSvc::ifRealField() const { return m_ifRealField; }
2857
2858//=============================================================================
2859// routine to fill the field vector
2860//=============================================================================
2861void MagneticFieldSvc::fieldGrid( const HepPoint3D& r, HepVector3D& bf ) const {
2862
2863 bf[0] = 0.0;
2864 bf[1] = 0.0;
2865 bf[2] = 0.0;
2866
2867 /// Linear interpolated field
2868 double z = r.z() - m_zOffSet;
2869 if ( z < m_min_FL[2] || z > m_max_FL[2] ) return;
2870 double x = r.x();
2871 if ( x < m_min_FL[0] || x > m_max_FL[0] ) return;
2872 double y = r.y();
2873 if ( y < m_min_FL[1] || y > m_max_FL[1] ) return;
2874 int i = int( ( x - m_min_FL[0] ) / m_Dxyz[0] );
2875 int j = int( ( y - m_min_FL[1] ) / m_Dxyz[1] );
2876 int k = int( ( z - m_min_FL[2] ) / m_Dxyz[2] );
2877
2878 int ijk000 = 3 * ( m_Nxyz[0] * ( m_Nxyz[1] * k + j ) + i );
2879 int ijk001 = 3 * ( m_Nxyz[0] * ( m_Nxyz[1] * ( k + 1 ) + j ) + i );
2880 int ijk010 = 3 * ( m_Nxyz[0] * ( m_Nxyz[1] * k + j + 1 ) + i );
2881 int ijk011 = 3 * ( m_Nxyz[0] * ( m_Nxyz[1] * ( k + 1 ) + j + 1 ) + i );
2882 int ijk100 = 3 * ( m_Nxyz[0] * ( m_Nxyz[1] * k + j ) + i + 1 );
2883 int ijk101 = 3 * ( m_Nxyz[0] * ( m_Nxyz[1] * ( k + 1 ) + j ) + i + 1 );
2884 int ijk110 = 3 * ( m_Nxyz[0] * ( m_Nxyz[1] * k + j + 1 ) + i + 1 );
2885 int ijk111 = 3 * ( m_Nxyz[0] * ( m_Nxyz[1] * ( k + 1 ) + j + 1 ) + i + 1 );
2886
2887 // auxiliary variables defined at the vertices of the cube that
2888 // contains the (x, y, z) point where the field is interpolated
2889 /* std::cout<<"x,y,z: "<<x<<","<<y<<","<<z<<std::endl;
2890 std::cout<<"point1(x,y,z):
2891 "<<m_P[ijk000]<<","<<m_P[ijk000+1]<<","<<m_P[ijk000+2]<<std::endl;
2892 std::cout<<"point2(x,y,z):
2893 "<<m_P[ijk001]<<","<<m_P[ijk001+1]<<","<<m_P[ijk001+2]<<std::endl;
2894 std::cout<<"point3(x,y,z):
2895 "<<m_P[ijk010]<<","<<m_P[ijk010+1]<<","<<m_P[ijk010+2]<<std::endl;
2896 std::cout<<"point4(x,y,z):
2897 "<<m_P[ijk011]<<","<<m_P[ijk011+1]<<","<<m_P[ijk011+2]<<std::endl;
2898 std::cout<<"point5(x,y,z):
2899 "<<m_P[ijk100]<<","<<m_P[ijk100+1]<<","<<m_P[ijk100+2]<<std::endl;
2900 std::cout<<"point6(x,y,z):
2901 "<<m_P[ijk101]<<","<<m_P[ijk101+1]<<","<<m_P[ijk101+2]<<std::endl;
2902 std::cout<<"point7(x,y,z):
2903 "<<m_P[ijk110]<<","<<m_P[ijk110+1]<<","<<m_P[ijk110+2]<<std::endl;
2904 std::cout<<"point8(x,y,z):
2905 "<<m_P[ijk111]<<","<<m_P[ijk111+1]<<","<<m_P[ijk111+2]<<std::endl;
2906
2907 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])
2908 std::cout<<"FATALERRORX****************************"<<std::endl;
2909 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])
2910 std::cout<<"FATALERRORY***************************"<<std::endl;
2911 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])
2912 std::cout<<"FATALERRORZ****************************"<<std::endl;
2913 */
2914 double cx000 = m_Q[ijk000];
2915 double cx001 = m_Q[ijk001];
2916 double cx010 = m_Q[ijk010];
2917 double cx011 = m_Q[ijk011];
2918 double cx100 = m_Q[ijk100];
2919 double cx101 = m_Q[ijk101];
2920 double cx110 = m_Q[ijk110];
2921 double cx111 = m_Q[ijk111];
2922 double cy000 = m_Q[ijk000 + 1];
2923 double cy001 = m_Q[ijk001 + 1];
2924 double cy010 = m_Q[ijk010 + 1];
2925 double cy011 = m_Q[ijk011 + 1];
2926 double cy100 = m_Q[ijk100 + 1];
2927 double cy101 = m_Q[ijk101 + 1];
2928 double cy110 = m_Q[ijk110 + 1];
2929 double cy111 = m_Q[ijk111 + 1];
2930 double cz000 = m_Q[ijk000 + 2];
2931 double cz001 = m_Q[ijk001 + 2];
2932 double cz010 = m_Q[ijk010 + 2];
2933 double cz011 = m_Q[ijk011 + 2];
2934 double cz100 = m_Q[ijk100 + 2];
2935 double cz101 = m_Q[ijk101 + 2];
2936 double cz110 = m_Q[ijk110 + 2];
2937 double cz111 = m_Q[ijk111 + 2];
2938 double hx1 = ( x - m_min_FL[0] - i * m_Dxyz[0] ) / m_Dxyz[0];
2939 double hy1 = ( y - m_min_FL[1] - j * m_Dxyz[1] ) / m_Dxyz[1];
2940 double hz1 = ( z - m_min_FL[2] - k * m_Dxyz[2] ) / m_Dxyz[2];
2941 double hx0 = 1.0 - hx1;
2942 double hy0 = 1.0 - hy1;
2943 double hz0 = 1.0 - hz1;
2944 double h000 = hx0 * hy0 * hz0;
2945 if ( fabs( h000 ) > 0.0 && cx000 > 9.0e5 && cy000 > 9.0e5 && cz000 > 9.0e5 ) return;
2946
2947 double h001 = hx0 * hy0 * hz1;
2948 if ( fabs( h001 ) > 0.0 && cx001 > 9.0e5 && cy001 > 9.0e5 && cz001 > 9.0e5 ) return;
2949
2950 double h010 = hx0 * hy1 * hz0;
2951 if ( fabs( h010 ) > 0.0 && cx010 > 9.0e5 && cy010 > 9.0e5 && cz010 > 9.0e5 ) return;
2952
2953 double h011 = hx0 * hy1 * hz1;
2954 if ( fabs( h011 ) > 0.0 && cx011 > 9.0e5 && cy011 > 9.0e5 && cz011 > 9.0e5 ) return;
2955
2956 double h100 = hx1 * hy0 * hz0;
2957 if ( fabs( h100 ) > 0.0 && cx100 > 9.0e5 && cy100 > 9.0e5 && cz100 > 9.0e5 ) return;
2958
2959 double h101 = hx1 * hy0 * hz1;
2960 if ( fabs( h101 ) > 0.0 && cx101 > 9.0e5 && cy101 > 9.0e5 && cz101 > 9.0e5 ) return;
2961
2962 double h110 = hx1 * hy1 * hz0;
2963 if ( fabs( h110 ) > 0.0 && cx110 > 9.0e5 && cy110 > 9.0e5 && cz110 > 9.0e5 ) return;
2964
2965 double h111 = hx1 * hy1 * hz1;
2966 if ( fabs( h111 ) > 0.0 && cx111 > 9.0e5 && cy111 > 9.0e5 && cz111 > 9.0e5 ) return;
2967
2968 bf( 0 ) = ( cx000 * h000 + cx001 * h001 + cx010 * h010 + cx011 * h011 + cx100 * h100 +
2969 cx101 * h101 + cx110 * h110 + cx111 * h111 );
2970 bf( 1 ) = ( cy000 * h000 + cy001 * h001 + cy010 * h010 + cy011 * h011 + cy100 * h100 +
2971 cy101 * h101 + cy110 * h110 + cy111 * h111 );
2972 bf( 2 ) = ( cz000 * h000 + cz001 * h001 + cz010 * h010 + cz011 * h011 + cz100 * h100 +
2973 cz101 * h101 + cz110 * h110 + cz111 * h111 );
2974 return;
2975}
2976
2977//=============================================================================
2978// routine to fill the field vector
2979//=============================================================================
2980void MagneticFieldSvc::fieldGrid_TE( const HepPoint3D& r, HepVector3D& bf ) const {
2981
2982 bf[0] = 0.0;
2983 bf[1] = 0.0;
2984 bf[2] = 0.0;
2985
2986 /// Linear interpolated field
2987 double z = r.z() - m_zOffSet_TE;
2988 if ( fabs( z ) < m_min_FL_TE[2] || fabs( z ) > m_max_FL_TE[2] ) return;
2989 double x = r.x();
2990 if ( fabs( x ) < m_min_FL_TE[0] || fabs( x ) > m_max_FL_TE[0] ) return;
2991 double y = r.y();
2992 if ( fabs( y ) < m_min_FL_TE[1] || fabs( y ) > m_max_FL_TE[1] ) return;
2993 int i = int( ( fabs( x ) - m_min_FL_TE[0] ) / m_Dxyz_TE[0] );
2994 int j = int( ( fabs( y ) - m_min_FL_TE[1] ) / m_Dxyz_TE[1] );
2995 int k = int( ( fabs( z ) - m_min_FL_TE[2] ) / m_Dxyz_TE[2] );
2996
2997 int ijk000 = 3 * ( m_Nxyz_TE[0] * ( m_Nxyz_TE[1] * k + j ) + i );
2998 int ijk001 = 3 * ( m_Nxyz_TE[0] * ( m_Nxyz_TE[1] * ( k + 1 ) + j ) + i );
2999 int ijk010 = 3 * ( m_Nxyz_TE[0] * ( m_Nxyz_TE[1] * k + j + 1 ) + i );
3000 int ijk011 = 3 * ( m_Nxyz_TE[0] * ( m_Nxyz_TE[1] * ( k + 1 ) + j + 1 ) + i );
3001 int ijk100 = 3 * ( m_Nxyz_TE[0] * ( m_Nxyz_TE[1] * k + j ) + i + 1 );
3002 int ijk101 = 3 * ( m_Nxyz_TE[0] * ( m_Nxyz_TE[1] * ( k + 1 ) + j ) + i + 1 );
3003 int ijk110 = 3 * ( m_Nxyz_TE[0] * ( m_Nxyz_TE[1] * k + j + 1 ) + i + 1 );
3004 int ijk111 = 3 * ( m_Nxyz_TE[0] * ( m_Nxyz_TE[1] * ( k + 1 ) + j + 1 ) + i + 1 );
3005
3006 // auxiliary variables defined at the vertices of the cube that
3007 // contains the (x, y, z) point where the field is interpolated
3008 /* std::cout<<"x,y,z: "<<x<<","<<y<<","<<z<<std::endl;
3009 std::cout<<"point1(x,y,z):
3010 "<<m_P[ijk000]<<","<<m_P[ijk000+1]<<","<<m_P[ijk000+2]<<std::endl;
3011 std::cout<<"point2(x,y,z):
3012 "<<m_P[ijk001]<<","<<m_P[ijk001+1]<<","<<m_P[ijk001+2]<<std::endl;
3013 std::cout<<"point3(x,y,z):
3014 "<<m_P[ijk010]<<","<<m_P[ijk010+1]<<","<<m_P[ijk010+2]<<std::endl;
3015 std::cout<<"point4(x,y,z):
3016 "<<m_P[ijk011]<<","<<m_P[ijk011+1]<<","<<m_P[ijk011+2]<<std::endl;
3017 std::cout<<"point5(x,y,z):
3018 "<<m_P[ijk100]<<","<<m_P[ijk100+1]<<","<<m_P[ijk100+2]<<std::endl;
3019 std::cout<<"point6(x,y,z):
3020 "<<m_P[ijk101]<<","<<m_P[ijk101+1]<<","<<m_P[ijk101+2]<<std::endl;
3021 std::cout<<"point7(x,y,z):
3022 "<<m_P[ijk110]<<","<<m_P[ijk110+1]<<","<<m_P[ijk110+2]<<std::endl;
3023 std::cout<<"point8(x,y,z):
3024 "<<m_P[ijk111]<<","<<m_P[ijk111+1]<<","<<m_P[ijk111+2]<<std::endl;
3025 */
3026 double cx000 = m_Q_TE[ijk000];
3027 double cx001 = m_Q_TE[ijk001];
3028 double cx010 = m_Q_TE[ijk010];
3029 double cx011 = m_Q_TE[ijk011];
3030 double cx100 = m_Q_TE[ijk100];
3031 double cx101 = m_Q_TE[ijk101];
3032 double cx110 = m_Q_TE[ijk110];
3033 double cx111 = m_Q_TE[ijk111];
3034 double cy000 = m_Q_TE[ijk000 + 1];
3035 double cy001 = m_Q_TE[ijk001 + 1];
3036 double cy010 = m_Q_TE[ijk010 + 1];
3037 double cy011 = m_Q_TE[ijk011 + 1];
3038 double cy100 = m_Q_TE[ijk100 + 1];
3039 double cy101 = m_Q_TE[ijk101 + 1];
3040 double cy110 = m_Q_TE[ijk110 + 1];
3041 double cy111 = m_Q_TE[ijk111 + 1];
3042 double cz000 = m_Q_TE[ijk000 + 2];
3043 double cz001 = m_Q_TE[ijk001 + 2];
3044 double cz010 = m_Q_TE[ijk010 + 2];
3045 double cz011 = m_Q_TE[ijk011 + 2];
3046 double cz100 = m_Q_TE[ijk100 + 2];
3047 double cz101 = m_Q_TE[ijk101 + 2];
3048 double cz110 = m_Q_TE[ijk110 + 2];
3049 double cz111 = m_Q_TE[ijk111 + 2];
3050 double hx1 = ( fabs( x ) - m_min_FL_TE[0] - i * m_Dxyz_TE[0] ) / m_Dxyz_TE[0];
3051 double hy1 = ( fabs( y ) - m_min_FL_TE[1] - j * m_Dxyz_TE[1] ) / m_Dxyz_TE[1];
3052 double hz1 = ( fabs( z ) - m_min_FL_TE[2] - k * m_Dxyz_TE[2] ) / m_Dxyz_TE[2];
3053 double hx0 = 1.0 - hx1;
3054 double hy0 = 1.0 - hy1;
3055 double hz0 = 1.0 - hz1;
3056 double h000 = hx0 * hy0 * hz0;
3057 if ( fabs( h000 ) > 0.0 && cx000 > 9.0e5 && cy000 > 9.0e5 && cz000 > 9.0e5 ) return;
3058
3059 double h001 = hx0 * hy0 * hz1;
3060 if ( fabs( h001 ) > 0.0 && cx001 > 9.0e5 && cy001 > 9.0e5 && cz001 > 9.0e5 ) return;
3061
3062 double h010 = hx0 * hy1 * hz0;
3063 if ( fabs( h010 ) > 0.0 && cx010 > 9.0e5 && cy010 > 9.0e5 && cz010 > 9.0e5 ) return;
3064
3065 double h011 = hx0 * hy1 * hz1;
3066 if ( fabs( h011 ) > 0.0 && cx011 > 9.0e5 && cy011 > 9.0e5 && cz011 > 9.0e5 ) return;
3067
3068 double h100 = hx1 * hy0 * hz0;
3069 if ( fabs( h100 ) > 0.0 && cx100 > 9.0e5 && cy100 > 9.0e5 && cz100 > 9.0e5 ) return;
3070
3071 double h101 = hx1 * hy0 * hz1;
3072 if ( fabs( h101 ) > 0.0 && cx101 > 9.0e5 && cy101 > 9.0e5 && cz101 > 9.0e5 ) return;
3073
3074 double h110 = hx1 * hy1 * hz0;
3075 if ( fabs( h110 ) > 0.0 && cx110 > 9.0e5 && cy110 > 9.0e5 && cz110 > 9.0e5 ) return;
3076
3077 double h111 = hx1 * hy1 * hz1;
3078 if ( fabs( h111 ) > 0.0 && cx111 > 9.0e5 && cy111 > 9.0e5 && cz111 > 9.0e5 ) return;
3079
3080 bf( 0 ) = ( cx000 * h000 + cx001 * h001 + cx010 * h010 + cx011 * h011 + cx100 * h100 +
3081 cx101 * h101 + cx110 * h110 + cx111 * h111 );
3082 bf( 1 ) = ( cy000 * h000 + cy001 * h001 + cy010 * h010 + cy011 * h011 + cy100 * h100 +
3083 cy101 * h101 + cy110 * h110 + cy111 * h111 );
3084 bf( 2 ) = ( cz000 * h000 + cz001 * h001 + cz010 * h010 + cz011 * h011 + cz100 * h100 +
3085 cz101 * h101 + cz110 * h110 + cz111 * h111 );
3086
3087 if ( r.x() < 0. && r.y() >= 0. && r.z() > 0. ) { bf( 0 ) = -bf( 0 ); }
3088 else if ( r.x() <= 0. && r.y() < 0. && r.z() > 0. )
3089 {
3090 bf( 0 ) = -bf( 0 );
3091 bf( 1 ) = -bf( 1 );
3092 }
3093 else if ( r.x() > 0. && r.y() < 0. && r.z() > 0. ) { bf( 1 ) = -bf( 1 ); }
3094
3095 else if ( r.x() >= 0. && r.y() > 0. && r.z() < 0. )
3096 {
3097 bf( 0 ) = -bf( 0 );
3098 bf( 1 ) = -bf( 1 );
3099 }
3100 else if ( r.x() < 0. && r.y() >= 0. && r.z() < 0. ) { bf( 1 ) = -bf( 1 ); }
3101 else if ( r.x() <= 0. && r.y() < 0. && r.z() < 0. )
3102 {
3103 bf( 0 ) = bf( 0 );
3104 bf( 1 ) = bf( 1 );
3105 }
3106 else if ( r.x() > 0. && r.y() <= 0. && r.z() < 0. ) { bf( 0 ) = -bf( 0 ); }
3107
3108 return;
3109}
HepVector3D
HepGeom::Vector3D< double > HepVector3D
Definition Analysis/VertexFit/src/BField.cxx:13
DECLARE_COMPONENT
DECLARE_COMPONENT(BesBdkRc)
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition CosmicGenerator.h:36
sprintf
sprintf(cut, "kal_costheta0_em>-0.93&&kal_costheta0_em<0.93&&kal_pxy0_em>=0.05+%d*0.1&&kal_" "pxy0_em<0.15+%d*0.1&&NGch>=2", j, j)
x
Double_t x[10]
Definition DataBase/tau_mode.c:49
pi
double pi
Definition EvtPhokharaDef.hh:17
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition MagneticFieldSvc.h:11
HepVector3D
HepGeom::Vector3D< double > HepVector3D
Definition MagneticFieldSvc.h:12
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
MagneticFieldSvc::m_connect_run
FieldDBUtil::ConnectionDB * m_connect_run
Definition MagneticFieldSvc.h:199
MagneticFieldSvc::m_mapMagnetInfo
std::map< int, std::vector< double > > m_mapMagnetInfo
Definition MagneticFieldSvc.h:200
MagneticFieldSvc::getReferField
virtual double getReferField()
Definition MagneticFieldSvc.cxx:2845
MagneticFieldSvc::fieldVector
virtual StatusCode fieldVector(const HepPoint3D &xyz, HepVector3D &fvec) const
Definition MagneticFieldSvc.cxx:1166
MagneticFieldSvc::~MagneticFieldSvc
virtual ~MagneticFieldSvc()
Virtual destructor.
Definition MagneticFieldSvc.cxx:106
MagneticFieldSvc::finalize
virtual StatusCode finalize()
Finalise the service.
Definition MagneticFieldSvc.cxx:704
MagneticFieldSvc::ifRealField
virtual bool ifRealField() const
Definition MagneticFieldSvc.cxx:2856
MagneticFieldSvc::beamEnergy
std::vector< double > beamEnergy
Definition MagneticFieldSvc.h:202
MagneticFieldSvc::initialize
virtual StatusCode initialize()
Initialise the service (Inherited Service overrides).
Definition MagneticFieldSvc.cxx:133
MagneticFieldSvc::m_mapBeamEnergy
std::map< int, std::vector< double > > m_mapBeamEnergy
Definition MagneticFieldSvc.h:201
MagneticFieldSvc::handle
void handle(const Incident &)
Definition MagneticFieldSvc.cxx:734
MagneticFieldSvc::uniFieldVector
virtual StatusCode uniFieldVector(const HepPoint3D &xyz, HepVector3D &fvec) const
Definition MagneticFieldSvc.cxx:2792
MagneticFieldSvc::MagneticFieldSvc
MagneticFieldSvc(const std::string &name, ISvcLocator *svc)
Definition MagneticFieldSvc.cxx:46
MagneticFieldSvc::current
std::vector< double > current
Definition MagneticFieldSvc.h:203
MagneticFieldSvc::m_eventSvc
IDataProviderSvc * m_eventSvc
Definition MagneticFieldSvc.h:196
EvtCyclic3::c_str
char * c_str(Index i)
Definition EvtCyclic3.cc:240