BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
MagFieldReader.cxx
Go to the documentation of this file.
1#include "GaudiKernel/INTuple.h"
2#include "GaudiKernel/INTupleSvc.h"
3#include "GaudiKernel/MsgStream.h"
4#include "GaudiKernel/NTuple.h"
5#include "GaudiKernel/SmartDataPtr.h"
6
7#include "CLHEP/Geometry/Point3D.h"
8#include "CLHEP/Geometry/Vector3D.h"
9#include "CLHEP/Units/PhysicalConstants.h"
10
11#include "CLHEP/Random/RandFlat.h"
12#include "MagFieldReader.h"
13#include "MagneticFieldSvc/IBesMagFieldSvc.h"
14#include <fstream>
15#include <iostream>
16using namespace std;
17// Test with root
18#include "TAxis.h"
19#include <TCanvas.h>
20#include <TFile.h>
21#include <TGraph.h>
22#include <TGraph2D.h>
23#include <TStyle.h>
24
25#ifndef ENABLE_BACKWARDS_COMPATIBILITY
26// backwards compatibility will be enabled ONLY in CLHEP 1.9
27typedef HepGeom::Point3D<double> HepPoint3D;
28#endif
29#ifndef ENABLE_BACKWARDS_COMPATIBILITY
30typedef HepGeom::Vector3D<double> HepVector3D;
31#endif
32using namespace CLHEP;
33//-----------------------------------------------------------------------------
34// Implementation file for class : MagFieldReader
35//
36//-----------------------------------------------------------------------------
37// Declaration of the Algorithm Factory
38// static const AlgFactory<MagFieldReader> s_factory ;
39// const IAlgFactory& MagFieldReaderFactory = s_factory ;
40
41//=============================================================================
42// Standard constructor, initializes variables
43//=============================================================================
44MagFieldReader::MagFieldReader( const std::string& name, ISvcLocator* pSvcLocator )
45 : Algorithm( name, pSvcLocator ), m_pIMF( 0 ) {
46 declareProperty( "Zmin", m_zMin = -1200.0 * mm );
47 declareProperty( "Zmax", m_zMax = 1200.0 * mm );
48 declareProperty( "Step", m_step = 50.0 * mm );
49 declareProperty( "Xmin", m_xMin = -900.0 * mm );
50 declareProperty( "Xmax", m_xMax = 900.0 * mm );
51 declareProperty( "Ymin", m_yMin = -900.0 * mm );
52 declareProperty( "Ymax", m_yMax = 900.0 * mm );
53 declareProperty( "filename", m_filename = "rawmode3_out.dat" );
54}
55
56//=============================================================================
57// Initialisation. Check parameters
58//=============================================================================
59StatusCode MagFieldReader::initialize() {
60
61 MsgStream msg( msgSvc(), name() );
62
63 msg << MSG::INFO << "FieldReader intialize() has been called" << endmsg;
64 StatusCode status = service( "MagneticFieldSvc", m_pIMF, true );
65 if ( !status.isSuccess() )
66 {
67 msg << MSG::FATAL << "Unable to open Magnetic field service" << endmsg;
68 return StatusCode::FAILURE;
69 }
70
71 msg << MSG::DEBUG << " Book ntuple" << endmsg;
72
73 NTuplePtr nt1( ntupleSvc(), "FILE1/n1" );
74 if ( nt1 ) m_tuple1 = nt1;
75 else
76 {
77 m_tuple1 = ntupleSvc()->book( "FILE1/n1", CLID_ColumnWiseTuple, "Field" );
78 if ( m_tuple1 )
79 {
80 status = m_tuple1->addItem( "x", m_x );
81 status = m_tuple1->addItem( "y", m_y );
82 status = m_tuple1->addItem( "z", m_z );
83 status = m_tuple1->addItem( "r", m_r );
84 status = m_tuple1->addItem( "Bx", m_Bx );
85 status = m_tuple1->addItem( "By", m_By );
86 status = m_tuple1->addItem( "Bz", m_Bz );
87 status = m_tuple1->addItem( "SigmaBx", m_sigma_bx );
88 status = m_tuple1->addItem( "SigmaBy", m_sigma_by );
89 status = m_tuple1->addItem( "SigmaBz", m_sigma_bz );
90 }
91 else
92 {
93 msg << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple1 ) << endmsg;
94 return StatusCode::FAILURE;
95 }
96 }
97
98 NTuplePtr nt2( ntupleSvc(), "FILE1/n2" );
99 if ( nt2 ) m_tuple2 = nt2;
100 else
101 {
102 m_tuple2 = ntupleSvc()->book( "FILE1/n2", CLID_ColumnWiseTuple, "Field" );
103 if ( m_tuple2 )
104 {
105 status = m_tuple2->addItem( "x", m_x2 );
106 status = m_tuple2->addItem( "y", m_y2 );
107 status = m_tuple2->addItem( "z", m_z2 );
108 status = m_tuple2->addItem( "r", m_r2 );
109 status = m_tuple2->addItem( "Bx", m_Bx2 );
110 status = m_tuple2->addItem( "By", m_By2 );
111 status = m_tuple2->addItem( "Bz", m_Bz2 );
112 }
113 else
114 {
115 msg << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple2 ) << endmsg;
116 return StatusCode::FAILURE;
117 }
118 }
119
120 NTuplePtr nt3( ntupleSvc(), "FILE1/n3" );
121 if ( nt3 ) m_tuple3 = nt3;
122 else
123 {
124 m_tuple3 = ntupleSvc()->book( "FILE1/n3", CLID_ColumnWiseTuple, "Field" );
125 if ( m_tuple3 )
126 {
127 status = m_tuple3->addItem( "x", m_x3 );
128 status = m_tuple3->addItem( "y", m_y3 );
129 status = m_tuple3->addItem( "z", m_z3 );
130 status = m_tuple3->addItem( "r", m_r3 );
131 status = m_tuple3->addItem( "phi", m_phi3 );
132 status = m_tuple3->addItem( "Bx", m_Bx3 );
133 status = m_tuple3->addItem( "By", m_By3 );
134 status = m_tuple3->addItem( "Bz", m_Bz3 );
135 }
136 else
137 {
138 msg << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple3 ) << endmsg;
139 return StatusCode::FAILURE;
140 }
141 }
142
143 NTuplePtr nt4( ntupleSvc(), "FILE1/n4" );
144 if ( nt4 ) m_tuple4 = nt4;
145 else
146 {
147 m_tuple4 = ntupleSvc()->book( "FILE1/n4", CLID_ColumnWiseTuple, "Field" );
148 if ( m_tuple4 ) { status = m_tuple4->addItem( "time", m_time ); }
149 else
150 {
151 msg << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple4 ) << endmsg;
152 return StatusCode::FAILURE;
153 }
154 }
155
156 status = service( "BesTimerSvc", m_timersvc );
157 if ( status.isFailure() )
158 {
159 msg << MSG::ERROR << name() << ": Unable to locate BesTimer Service" << endmsg;
160 return StatusCode::FAILURE;
161 }
162 m_timer = m_timersvc->addItem( "Read field Time" );
163
164 msg << MSG::INFO << "MagFieldReader initialized" << endmsg;
165 return StatusCode::SUCCESS;
166}
167
168//=============================================================================
169// Main execution
170//=============================================================================
171StatusCode MagFieldReader::execute() {
172
173 MsgStream msg( msgSvc(), name() );
174
175 msg << MSG::DEBUG << "==> Execute MagFieldReader" << endmsg;
176 /*
177 //calculate the error of Bfield
178 std::ifstream infile( m_filename.c_str() );
179 if(!infile.good()) msg << MSG::FATAL << "Unable to read the file: " << m_filename <<
180 endmsg; int nLine = 0; char line[ 255 ]; while( infile ) {
181 // parse each line of the file,
182 // comment lines begin with '#' in the cdf file
183 infile.getline( line, 255 );
184 if ( line[0] == '#' ) continue;
185 std::string gridx, gridy, gridz, sFx, sFy, sFz;
186 char* token = strtok( line, " " );
187 if ( token ) { gridx = token; token = strtok( NULL, " " );} else continue;
188 if ( token ) { gridy = token; token = strtok( NULL, " " );} else continue;
189 if ( token ) { gridz = token; token = strtok( NULL, " " );} else continue;
190 if ( token ) { sFx = token; token = strtok( NULL, " " );} else continue;
191 if ( token ) { sFy = token; token = strtok( NULL, " " );} else continue;
192 if ( token ) { sFz = token; token = strtok( NULL, " " );} else continue;
193 if ( token != NULL ) continue;
194 nLine++;
195 double x,y,z,bx,by,bz;
196 //Grid position
197 x = atof( gridx.c_str() );
198 y = atof( gridy.c_str() );
199 z = atof( gridz.c_str() );
200 // Field values are given in tesla in CDF file. Convert to CLHEP units
201 bx = atof( sFx.c_str() );
202 by = atof( sFy.c_str() );
203 bz = atof( sFz.c_str() );
204 //std::cout<<"x,y,z: "<<x<<" "<<y<<" "<<z<<" "<<"bx,by,bz: "<<bx<<" "<<by<<"
205 "<<bz<<std::endl; x = -x*m; y = y*m; z = -z*m;
206
207 HepPoint3D r(x,y,z);
208 HepVector3D b;
209 m_pIMF->fieldVector(r,b);
210 m_Bx = b.x()/tesla;
211 m_By = b.y()/tesla;
212 m_Bz = b.z()/tesla;
213 m_sigma_bx = (b.x()/tesla + bx)*10000.;
214 m_sigma_by = (b.y()/tesla - by)*10000.;
215 m_sigma_bz = (b.z()/tesla + bz)*10000.;
216 m_r = std::sqrt(x/m*x/m+y/m*y/m);
217 m_x = x/m;
218 m_y = y/m;
219 m_z = z/m;
220 m_tuple1->write();
221 }
222 infile.close();
223 std::cout<<"Totally "<<nLine<<" in file "<<m_filename<<std::endl;
224
225 for(int k = 0; k < 5; k++) {
226 double rr;
227 if(k==0) rr = 0;
228 if(k==1) rr = 200;
229 if(k==2) rr = 400;
230 if(k==3) rr = 600;
231 if(k==4) rr = 800;
232 for(int zz = -1195; zz <1200; zz+=50) {
233 double bxt = 0.,byt= 0.,bzt = 0.;
234 for(int i = 0; i < 100000; i++) {
235 double phi = CLHEP::twopi*CLHEP::RandFlat::shoot();
236 double xx = rr*std::cos(phi);
237 double yy = rr*std::sin(phi);
238 HepPoint3D r(xx,yy,double(zz));
239 HepVector3D b;
240 m_pIMF->fieldVector(r,b);
241 bxt+=b.x()/tesla;
242 byt+=b.y()/tesla;
243 bzt+=b.z()/tesla;
244 }
245 m_z2 = zz;
246 m_r2 = rr;
247 m_Bx2 = bxt/100000;
248 m_By2 = byt/100000;
249 m_Bz2 = bzt/100000;
250 m_tuple2->write();
251 }
252 }
253
254 for(int k = 0; k < 3; k++) {
255 double zz;
256 if(k==0) zz = 0;
257 if(k==1) zz = -800;
258 if(k==2) zz = 800;
259 for(int rr = 200; rr <=600; rr+=400) {
260 for(double phi = 0; phi <= 2*3.14159; phi+=0.1745) {
261 //double phi = CLHEP::twopi*CLHEP::RandFlat::shoot();
262 double xx = rr*std::cos(phi);
263 double yy = rr*std::sin(phi);
264 HepPoint3D r(xx,yy,double(zz));
265 HepVector3D b;
266 m_pIMF->fieldVector(r,b);
267 m_Bx3 = b.x()/tesla;
268 m_By3 = b.y()/tesla;
269 m_Bz3 = b.z()/tesla;
270 m_phi3 = phi*180/CLHEP::pi;
271 m_x3 = xx;
272 m_y3 = yy;
273 m_z3 = zz;
274 m_r3 = rr;
275 m_tuple3->write();
276 }
277 }
278 }
279 */
280 m_timer->start();
281 for ( int i = 0; i < 20000; i++ )
282 {
283 double px, py, pz, bx, by, bz;
284 double max_x = 1.8 * m, min_x = -1.8 * m, max_y = 1.8 * m, min_y = -1.8 * m,
285 max_z = 2.1 * m, min_z = -2.1 * m;
286 px = min_x + ( max_x - min_x ) * RandFlat::shoot();
287 py = min_y + ( max_y - min_y ) * RandFlat::shoot();
288 pz = min_z + ( max_z - min_z ) * RandFlat::shoot();
289 HepPoint3D r( px, py, pz );
290 HepVector3D b;
291 m_pIMF->fieldVector( r, b );
292 bx = b.x() / tesla;
293 by = b.y() / tesla;
294 bz = b.z() / tesla;
295 }
296 m_timer->stop();
297 m_time = m_timer->elapsed();
298 m_tuple4->write();
299 /*
300 StatusCode sc = this->readField();
301 if( sc.isFailure() ) {
302 msg << MSG::FATAL << "Unable to execute MagFieldReader"
303 << endmsg;
304 return sc;
305 }
306 */
307 msg << MSG::INFO << "MagFieldReader executed" << endmsg;
308 return StatusCode::SUCCESS;
309}
310
311//=============================================================================
312// Finalize
313//=============================================================================
314StatusCode MagFieldReader::finalize() {
315
316 MsgStream msg( msgSvc(), name() );
317 msg << MSG::DEBUG << "==> Finalize MagFieldReader" << endmsg;
318
319 return StatusCode::SUCCESS;
320}
321
322//=============================================================================
323//=============================================================================
324// The MagFieldReader print out info messages
325// with the field value at different locations.
326//=============================================================================
327StatusCode MagFieldReader::readField() {
328
329 MsgStream msg( msgSvc(), name() );
330 msg << MSG::DEBUG << "m_pIMF = " << m_pIMF << endmsg;
331
332 double referfield = m_pIMF->getReferField() / tesla;
333 msg << MSG::DEBUG << "The reference field is " << referfield << " tesla" << endmsg;
334 bool ifrealfield = m_pIMF->ifRealField();
335 if ( ifrealfield ) msg << MSG::DEBUG << "Use the real field" << endmsg;
336 else msg << MSG::DEBUG << "Use the fake field" << endmsg;
337
338 HepVector3D B( 0.0, 0.0, 0.0 );
339 /* for ( double z = m_zMin; z <= m_zMax; z += m_step ) {
340 for( double y = m_yMin; y <= m_yMax; y += m_step ) {
341 for( double x = m_xMin; x <= m_xMax; x += m_step ) {
342 HepPoint3D P( x, y, z );
343 // get field at point P
344 m_pIMF->fieldVector( P, B );
345 // fill ntuple
346 m_x = P.x()/cm;
347 m_y = P.y()/cm;
348 m_z = P.z()/cm;
349 m_Bx = B.x()/tesla;
350 m_By = B.y()/tesla;
351 m_Bz = B.z()/tesla;
352 m_ntuple->write();
353 }
354 }
355 HepPoint3D P0( 0.0, 0.0, z );
356 m_pIMF->fieldVector( P0, B );
357 msg << MSG::DEBUG << "Magnetic Field at ("
358 << P0.x() << ", " << P0.y() << ", " << P0.z() << " ) = "
359 << (B.x())/tesla << ", "
360 << (B.y())/tesla << ", "
361 << (B.z())/tesla << " Tesla "
362 << endmsg;
363 }*/
364 // magnetic field check bz vs z(m_zMin,m_zMax), x=0,200,400,600,800mm, y=0
365 const int n1 = 241;
366 double x[n1], y[n1], y1[n1], y2[n1], y3[n1], y4[n1];
367 // magnetic field check bz,br,bendp,bphi vs z(0,m_zMax), x=50,100,200,400,600,800mm, y=0
368 const int n2 = 121;
369 double x1[n2], bz1[n2], bz2[n2], bz3[n2], bz4[n2], bz5[n2], bz6[n2];
370 double br1[n2], br2[n2], br3[n2], br4[n2], br5[n2], br6[n2];
371 double bp1[n2], bp2[n2], bp3[n2], bp4[n2], bp5[n2], bp6[n2];
372 double bphi1[n2], bphi2[n2], bphi3[n2], bphi4[n2], bphi5[n2], bphi6[n2];
373 // magnetic field check bx,by vs y(m_yMin,m_yMax)
374 const int n3 = 161;
375 double x2[n3], bx1[n3], bx2[n3], bx3[n3], bx4[n3], bx5[n3], bx6[n3];
376 double by1[n3], by2[n3], by3[n3], by4[n3], by5[n3], by6[n3];
377 // check globle field value bz vs x (-2.6m,2.6m),y=0,z=0.
378 const int n4 = 521;
379 double globle_x[n4], globle_bz[n4];
380 int i = 0;
381 double px = 0.0, py = 0.0, pz = 0.0;
382 HepPoint3D pos( px, py, pz );
383 for ( px = -2600; px <= 2600; px += 10 )
384 {
385 pos[0] = px;
386 pos[1] = 0;
387 pos[2] = 0;
388 m_pIMF->fieldVector( pos, B );
389 globle_x[i] = px / m;
390 globle_bz[i] = B.z() / tesla;
391 i++;
392 }
393 i = 0;
394 for ( pz = m_zMin; pz <= m_zMax; pz += 10 )
395 {
396 pos[0] = 0;
397 pos[1] = 0;
398 pos[2] = pz;
399 m_pIMF->fieldVector( pos, B );
400 x[i] = pz / m;
401 y[i] = B.z() / tesla;
402
403 pos[0] = 200;
404 pos[1] = 0;
405 pos[2] = pz;
406 m_pIMF->fieldVector( pos, B );
407 y1[i] = B.z() / tesla;
408
409 pos[0] = 400;
410 pos[1] = 0;
411 pos[2] = pz;
412 m_pIMF->fieldVector( pos, B );
413 y2[i] = B.z() / tesla;
414
415 pos[0] = 600;
416 pos[1] = 0;
417 pos[2] = pz;
418 m_pIMF->fieldVector( pos, B );
419 y3[i] = B.z() / tesla;
420
421 pos[0] = 800;
422 pos[1] = 0;
423 pos[2] = pz;
424 m_pIMF->fieldVector( pos, B );
425 y4[i] = B.z() / tesla;
426 i++;
427 }
428 int j = 0;
429 double tbx, tby, tbz, tbr, tbp, tbphi;
430 for ( pz = 0; pz <= m_zMax; pz += 10 )
431 {
432 pos[0] = 50;
433 pos[1] = 0;
434 pos[2] = pz;
435 m_pIMF->fieldVector( pos, B );
436 x1[j] = pz / m;
437 tbx = B.x() / tesla;
438 tby = B.y() / tesla;
439 tbz = B.z() / tesla;
440 tbr = tbx * pos[0] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] ) +
441 tby * pos[1] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] );
442 tbphi = tbx * pos[1] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] ) +
443 tby * pos[0] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] );
444 tbp = tbz - tbr * pz / sqrt( pos[0] * pos[0] + pos[1] * pos[1] );
445 bz1[j] = tbz;
446 br1[j] = tbr * 10000;
447 bp1[j] = tbp;
448 bphi1[j] = tbphi * 10000;
449
450 pos[0] = 100;
451 pos[1] = 0;
452 pos[2] = pz;
453 m_pIMF->fieldVector( pos, B );
454 tbx = B.x() / tesla;
455 tby = B.y() / tesla;
456 tbz = B.z() / tesla;
457 tbr = tbx * pos[0] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] ) +
458 tby * pos[1] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] );
459 tbphi = tbx * pos[1] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] ) +
460 tby * pos[0] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] );
461 tbp = tbz - tbr * pz / sqrt( pos[0] * pos[0] + pos[1] * pos[1] );
462 bz2[j] = tbz;
463 br2[j] = tbr * 10000;
464 bp2[j] = tbp;
465 bphi2[j] = tbphi * 10000;
466
467 pos[0] = 200;
468 pos[1] = 0;
469 pos[2] = pz;
470 m_pIMF->fieldVector( pos, B );
471 tbx = B.x() / tesla;
472 tby = B.y() / tesla;
473 tbz = B.z() / tesla;
474 tbr = tbx * pos[0] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] ) +
475 tby * pos[1] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] );
476 tbphi = tbx * pos[1] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] ) +
477 tby * pos[0] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] );
478 tbp = tbz - tbr * pz / sqrt( pos[0] * pos[0] + pos[1] * pos[1] );
479 bz3[j] = tbz;
480 br3[j] = tbr * 10000;
481 bp3[j] = tbp;
482 bphi3[j] = tbphi * 10000;
483
484 pos[0] = 400;
485 pos[1] = 0;
486 pos[2] = pz;
487 m_pIMF->fieldVector( pos, B );
488 tbx = B.x() / tesla;
489 tby = B.y() / tesla;
490 tbz = B.z() / tesla;
491 tbr = tbx * pos[0] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] ) +
492 tby * pos[1] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] );
493 tbphi = tbx * pos[1] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] ) +
494 tby * pos[0] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] );
495 tbp = tbz - tbr * pz / sqrt( pos[0] * pos[0] + pos[1] * pos[1] );
496 bz4[j] = tbz;
497 br4[j] = tbr * 10000;
498 bp4[j] = tbp;
499 bphi4[j] = tbphi * 10000;
500
501 pos[0] = 600;
502 pos[1] = 0;
503 pos[2] = pz;
504 m_pIMF->fieldVector( pos, B );
505 tbx = B.x() / tesla;
506 tby = B.y() / tesla;
507 tbz = B.z() / tesla;
508 tbr = tbx * pos[0] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] ) +
509 tby * pos[1] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] );
510 tbphi = tbx * pos[1] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] ) +
511 tby * pos[0] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] );
512 tbp = tbz - tbr * pz / sqrt( pos[0] * pos[0] + pos[1] * pos[1] );
513 bz5[j] = tbz;
514 br5[j] = tbr * 10000;
515 bp5[j] = tbp;
516 bphi5[j] = tbphi * 10000;
517
518 pos[0] = 800;
519 pos[1] = 0;
520 pos[2] = pz;
521 m_pIMF->fieldVector( pos, B );
522 tbx = B.x() / tesla;
523 tby = B.y() / tesla;
524 tbz = B.z() / tesla;
525 tbr = tbx * pos[0] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] ) +
526 tby * pos[1] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] );
527 tbphi = tbx * pos[1] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] ) +
528 tby * pos[0] / sqrt( pos[0] * pos[0] + pos[1] * pos[1] );
529 tbp = tbz - tbr * pz / sqrt( pos[0] * pos[0] + pos[1] * pos[1] );
530 bz6[j] = tbz;
531 br6[j] = tbr * 10000;
532 bp6[j] = tbp;
533 bphi6[j] = tbphi * 10000;
534 j++;
535 }
536 int l = 0;
537 for ( py = m_yMin; py <= m_yMax; py += 10 )
538 {
539 pos[0] = 0;
540 pos[1] = py;
541 pos[2] = 0;
542 m_pIMF->fieldVector( pos, B );
543 tbx = B.x() / tesla * 10000;
544 tby = B.y() / tesla * 10000;
545 tbz = B.z() / tesla * 10000;
546 x2[l] = py / m;
547 bx1[l] = tbx;
548 by1[l] = tby;
549
550 pos[0] = 100;
551 pos[1] = py;
552 pos[2] = 100;
553 m_pIMF->fieldVector( pos, B );
554 tbx = B.x() / tesla * 10000;
555 tby = B.y() / tesla * 10000;
556 tbz = B.z() / tesla * 10000;
557 bx2[l] = tbx;
558 by2[l] = tby;
559
560 pos[0] = 400;
561 pos[1] = py;
562 pos[2] = 400;
563 m_pIMF->fieldVector( pos, B );
564 tbx = B.x() / tesla * 10000;
565 tby = B.y() / tesla * 10000;
566 tbz = B.z() / tesla * 10000;
567 bx3[l] = tbx;
568 by3[l] = tby;
569 l++;
570 }
571
572 TGraph2D* dt = new TGraph2D();
573 int k = 0;
574 for ( pz = -3000; pz <= 3000; pz += 50 )
575 {
576 for ( py = -2700; py <= 2700; py += 50 )
577 {
578 pos[0] = 0;
579 pos[1] = py;
580 pos[2] = pz;
581 m_pIMF->fieldVector( pos, B );
582 tbx = B.x() / tesla;
583 tby = B.y() / tesla;
584 tbz = B.z() / tesla;
585 dt->SetPoint( k, pz / m, py / m, tbz );
586 k++;
587 }
588 }
589 TGraph2D* dt1 = new TGraph2D();
590 k = 0;
591 for ( pz = -3000; pz <= 3000; pz += 50 )
592 {
593 for ( py = -3000; py <= 3000; py += 50 )
594 {
595 pos[0] = 0;
596 pos[1] = py;
597 pos[2] = pz;
598 m_pIMF->fieldVector( pos, B );
599 tbx = B.x() / tesla;
600 tby = B.y() / tesla;
601 tbz = B.z() / tesla;
602 double btot = sqrt( tbx * tbx + tby * tby + tbz * tbz );
603 dt1->SetPoint( k, pz / m, py / m, btot );
604 k++;
605 }
606 }
607 TGraph2D* dt2 = new TGraph2D();
608 k = 0;
609 for ( pz = m_zMin; pz <= m_zMax; pz += 50 )
610 {
611 for ( py = m_yMin; py <= m_yMax; py += 50 )
612 {
613 pos[0] = 0;
614 pos[1] = py;
615 pos[2] = pz;
616 m_pIMF->fieldVector( pos, B );
617 tbx = B.x() / tesla;
618 tby = B.y() / tesla;
619 tbz = B.z() / tesla;
620 // double btot = sqrt(tbx*tbx+tby*tby+tbz*tbz);
621 dt2->SetPoint( k, pz / m, py / m, tbz );
622 k++;
623 }
624 }
625 gStyle->SetPalette( 1 );
626 gStyle->SetOptTitle( 1 );
627 gStyle->SetOptStat( 0 );
628
629 TFile* f1 = new TFile( "graph.root", "RECREATE" );
630 TGraph* gr1 = new TGraph( n1, x, y );
631 TGraph* gr2 = new TGraph( n1, x, y1 );
632 TGraph* gr3 = new TGraph( n1, x, y2 );
633 TGraph* gr4 = new TGraph( n1, x, y3 );
634 TGraph* gr5 = new TGraph( n1, x, y4 );
635
636 TGraph* g1 = new TGraph( n2, x1, bz1 );
637 TGraph* g2 = new TGraph( n2, x1, bz2 );
638 TGraph* g3 = new TGraph( n2, x1, bz3 );
639 TGraph* g4 = new TGraph( n2, x1, bz4 );
640 TGraph* g5 = new TGraph( n2, x1, bz5 );
641 TGraph* g6 = new TGraph( n2, x1, bz6 );
642
643 TGraph* g7 = new TGraph( n2, x1, br1 );
644 TGraph* g8 = new TGraph( n2, x1, br2 );
645 TGraph* g9 = new TGraph( n2, x1, br3 );
646 TGraph* g10 = new TGraph( n2, x1, br4 );
647 TGraph* g11 = new TGraph( n2, x1, br5 );
648 TGraph* g12 = new TGraph( n2, x1, br6 );
649
650 TGraph* g13 = new TGraph( n2, x1, bp1 );
651 TGraph* g14 = new TGraph( n2, x1, bp2 );
652 TGraph* g15 = new TGraph( n2, x1, bp3 );
653 TGraph* g16 = new TGraph( n2, x1, bp4 );
654 TGraph* g17 = new TGraph( n2, x1, bp5 );
655 TGraph* g18 = new TGraph( n2, x1, bp6 );
656
657 TGraph* g19 = new TGraph( n2, x1, bphi1 );
658 TGraph* g20 = new TGraph( n2, x1, bphi2 );
659 TGraph* g21 = new TGraph( n2, x1, bphi3 );
660 TGraph* g22 = new TGraph( n2, x1, bphi4 );
661 TGraph* g23 = new TGraph( n2, x1, bphi5 );
662 TGraph* g24 = new TGraph( n2, x1, bphi6 );
663
664 TGraph* g25 = new TGraph( n3, x2, bx1 );
665 TGraph* g26 = new TGraph( n3, x2, bx2 );
666 TGraph* g27 = new TGraph( n3, x2, bx3 );
667
668 TGraph* g28 = new TGraph( n3, x2, by1 );
669 TGraph* g29 = new TGraph( n3, x2, by2 );
670 TGraph* g30 = new TGraph( n3, x2, by3 );
671
672 TGraph* g31 = new TGraph( n4, globle_x, globle_bz );
673
674 TCanvas* c1 = new TCanvas( "c1", "Two Graphs", 200, 10, 600, 400 );
675 TCanvas* c2 = new TCanvas( "c2", "Two Graphs", 200, 10, 600, 400 );
676 c1->Divide( 2, 2 );
677 c2->Divide( 2, 2 );
678 c1->cd( 1 );
679 gr1->SetLineColor( 2 );
680 gr1->SetLineWidth( 2 );
681 gr1->Draw( "ALP" );
682 gr1->SetTitle( "bz vs z (y=0,x=0,200,400,600,800mm)" );
683 gr1->GetXaxis()->SetTitle( "m" );
684 gr1->GetYaxis()->SetTitle( "tesla" );
685 gr2->SetLineWidth( 2 );
686 gr2->SetLineColor( 3 );
687 gr2->Draw( "LP" );
688 gr3->SetLineWidth( 2 );
689 gr3->SetLineColor( 4 );
690 gr3->Draw( "LP" );
691 gr4->SetLineWidth( 2 );
692 gr4->SetLineColor( 5 );
693 gr4->Draw( "LP" );
694 gr5->SetLineWidth( 2 );
695 gr5->SetLineColor( 6 );
696 gr5->Draw( "LP" );
697
698 c1->cd( 2 );
699 g1->SetLineColor( 2 );
700 g1->SetLineWidth( 2 );
701 g1->Draw( "ALP" );
702 g1->SetTitle( "bz(red),bendp vs z (y=0,x=50,100,200,400,600,800mm)" );
703 g1->GetXaxis()->SetTitle( "m" );
704 g1->GetYaxis()->SetTitle( "tesla" );
705 g1->GetYaxis()->SetRangeUser( 0.2, 2 );
706 g2->SetLineWidth( 2 );
707 g2->SetLineColor( 2 );
708 g2->Draw( "LP" );
709 g3->SetLineWidth( 2 );
710 g3->SetLineColor( 2 );
711 g3->Draw( "LP" );
712 g4->SetLineWidth( 2 );
713 g4->SetLineColor( 2 );
714 g4->Draw( "LP" );
715 g5->SetLineWidth( 2 );
716 g5->SetLineColor( 2 );
717 g5->Draw( "LP" );
718 g6->SetLineColor( 2 );
719 g6->SetLineWidth( 2 );
720 g6->Draw( "LP" );
721
722 g13->SetLineWidth( 2 );
723 g13->SetLineColor( 4 );
724 g13->Draw( "LP" );
725 g14->SetLineWidth( 2 );
726 g14->SetLineColor( 4 );
727 g14->Draw( "LP" );
728 g15->SetLineWidth( 2 );
729 g15->SetLineColor( 4 );
730 g15->Draw( "LP" );
731 g16->SetLineWidth( 2 );
732 g16->SetLineColor( 4 );
733 g16->Draw( "LP" );
734 g17->SetLineWidth( 2 );
735 g17->SetLineColor( 4 );
736 g17->Draw( "LP" );
737 g18->SetLineWidth( 2 );
738 g18->SetLineColor( 4 );
739 g18->Draw( "LP" );
740
741 c1->cd( 3 );
742 g7->SetLineWidth( 2 );
743 g7->SetLineColor( 2 );
744 g7->Draw( "ALP" );
745 g7->SetTitle( "br vs z (y=0,x=50,100,200,400,600,800mm)" );
746 g7->GetXaxis()->SetTitle( "m" );
747 g7->GetYaxis()->SetTitle( "gauss" );
748 g7->GetYaxis()->SetRangeUser( -1100, 1000 );
749 g8->SetLineWidth( 2 );
750 g8->SetLineColor( 3 );
751 g8->Draw( "LP" );
752 g9->SetLineWidth( 2 );
753 g9->SetLineColor( 4 );
754 g9->Draw( "LP" );
755 g10->SetLineWidth( 2 );
756 g10->SetLineColor( 5 );
757 g10->Draw( "LP" );
758 g11->SetLineWidth( 2 );
759 g11->SetLineColor( 6 );
760 g11->Draw( "LP" );
761 g12->SetLineWidth( 2 );
762 g12->SetLineColor( 7 );
763 g12->Draw( "LP" );
764
765 c1->cd( 4 );
766 g19->SetLineWidth( 2 );
767 g19->SetLineColor( 2 );
768 g19->Draw( "ALP" );
769 g19->SetTitle( "bphi vs z (y=0,x=50,100,200,400,600,800mm)" );
770 g19->GetXaxis()->SetTitle( "m" );
771 g19->GetYaxis()->SetTitle( "gauss" );
772 g19->GetYaxis()->SetRangeUser( -1000, 200 );
773 g20->SetLineWidth( 2 );
774 g20->SetLineColor( 3 );
775 g20->Draw( "LP" );
776 g21->SetLineWidth( 2 );
777 g21->SetLineColor( 4 );
778 g21->Draw( "LP" );
779 g22->SetLineWidth( 2 );
780 g22->SetLineColor( 5 );
781 g22->Draw( "LP" );
782 g23->SetLineWidth( 2 );
783 g23->SetLineColor( 6 );
784 g23->Draw( "LP" );
785 g24->SetLineWidth( 2 );
786 g24->SetLineColor( 7 );
787 g24->Draw( "LP" );
788
789 TCanvas* c3 = new TCanvas( "c3", "Two Graphs", 200, 10, 600, 400 );
790 c3->cd( 1 );
791 // dt->Draw("surf2");
792 // dt->Draw("tril p3");
793 // dt->Draw("surf1z");
794 dt->Draw( "z sinusoidal" );
795 dt->SetTitle( "bz vs y,z (x=0)" );
796
797 TCanvas* c4 = new TCanvas( "c4", "Two Graphs", 200, 10, 600, 400 );
798 c4->cd( 1 );
799 dt1->Draw( "z sinusoidal" );
800 dt1->SetTitle( "btot vs y,z (x=0)" );
801
802 TCanvas* c5 = new TCanvas( "c5", "Two Graphs", 200, 10, 600, 400 );
803 c5->cd( 1 );
804 dt2->Draw( "z sinusoidal" );
805 dt2->SetTitle( "bz vs y,z (x=0)" );
806
807 c2->cd( 2 );
808 g25->SetLineWidth( 2 );
809 g25->SetLineColor( 2 );
810 g25->Draw( "ALP" );
811 g25->SetTitle( "bx(red),by vs y (x,z=0,100,400mm)" );
812 g25->GetXaxis()->SetTitle( "m" );
813 g25->GetYaxis()->SetTitle( "gauss" );
814 g25->GetYaxis()->SetRangeUser( -200, 300 );
815 g26->SetLineWidth( 2 );
816 g26->SetLineColor( 2 );
817 g26->Draw( "LP" );
818 g27->SetLineWidth( 2 );
819 g27->SetLineColor( 2 );
820 g27->Draw( "LP" );
821
822 g28->SetLineWidth( 2 );
823 g28->SetLineColor( 3 );
824 g28->Draw( "LP" );
825 g29->SetLineWidth( 2 );
826 g29->SetLineColor( 3 );
827 g29->Draw( "LP" );
828 g30->SetLineWidth( 2 );
829 g30->SetLineColor( 3 );
830 g30->Draw( "LP" );
831
832 c2->cd( 3 );
833 g31->SetLineWidth( 2 );
834 g31->SetLineColor( 2 );
835 g31->Draw( "ALP" );
836 g31->SetTitle( "bz vs x (y,z=0)" );
837 g31->GetXaxis()->SetTitle( "m" );
838 g31->GetYaxis()->SetTitle( "tesla" );
839
840 c1->Write();
841 c2->Write();
842 c3->Write();
843 c4->Write();
844 c5->Write();
845 /* HepVector3D B1(0.0,0.0,0.0);
846 int ii=0;
847 double posx,posy,posz;
848 ofstream outfile("tmp.txt",ios_base::app);
849 for ( double posz = -2800;posz <= 2800; posz += 10 ) {
850 for(double posy = -2650; posy <= 2650; posy += 10){
851 posx=1500; posy = 0;
852 HepPoint3D p1(posx,posy,posz);
853 m_pIMF->fieldVector( p1, B1 );
854 tbx=B1.x()/tesla;
855 tby=B1.y()/tesla;
856 tbz=B1.z()/tesla;
857 tbr=tbx*posx/sqrt(posx*posx+posy*posy)+tby*posy/sqrt(posx*posx+posy*posy);
858 tbp=tbz-tbr*posz/sqrt(posx*posx+posy*posy);
859 outfile<<posz/m<<" "<<tby<<endl;
860 ii++;
861 }
862 }
863 outfile.close();
864 */
865 // Return status code.
866 return StatusCode::SUCCESS;
867}
HepVector3D
HepGeom::Vector3D< double > HepVector3D
Definition Analysis/VertexFit/src/BField.cxx:13
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition CosmicGenerator.h:36
f1
TFile * f1
Definition DataBase/tau_mode.c:4
g1
TF1 * g1
Definition DataBase/tau_mode.c:58
x
Double_t x[10]
Definition DataBase/tau_mode.c:49
dt
TGraph2DErrors * dt
Definition McCor.cxx:41
n2
int n2
Definition SD0Tag.cxx:59
n1
int n1
Definition SD0Tag.cxx:58
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:51
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
IBesMagFieldSvc::fieldVector
virtual StatusCode fieldVector(const HepGeom::Point3D< double > &xyz, HepGeom::Vector3D< double > &fvec) const =0
IBesMagFieldSvc::getReferField
virtual double getReferField()=0
IBesMagFieldSvc::ifRealField
virtual bool ifRealField() const =0
MagFieldReader::MagFieldReader
MagFieldReader(const std::string &name, ISvcLocator *pSvcLocator)
Standard constructor.
Definition MagFieldReader.cxx:44
MagFieldReader::execute
virtual StatusCode execute()
Algorithm execution.
Definition MagFieldReader.cxx:171
MagFieldReader::finalize
virtual StatusCode finalize()
Algorithm finalization.
Definition MagFieldReader.cxx:314
MagFieldReader::initialize
virtual StatusCode initialize()
Destructor.
Definition MagFieldReader.cxx:59