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