00001 #include "GaudiKernel/AlgFactory.h"
00002 #include "GaudiKernel/MsgStream.h"
00003 #include "GaudiKernel/SmartDataPtr.h"
00004
00005 #include "GaudiKernel/NTuple.h"
00006 #include "GaudiKernel/INTuple.h"
00007 #include "GaudiKernel/INTupleSvc.h"
00008 #include "GaudiKernel/SmartDataPtr.h"
00009
00010 #include "CLHEP/Geometry/Vector3D.h"
00011 #include "CLHEP/Geometry/Point3D.h"
00012 #include "CLHEP/Units/PhysicalConstants.h"
00013
00014 #include "MagneticField/IMagneticFieldSvc.h"
00015 #include "MagneticField/MagFieldReader.h"
00016 #include <fstream>
00017 #include <iostream>
00018 #include "CLHEP/Random/RandFlat.h"
00019 using namespace std;
00020
00021 #include <TGraph.h>
00022 #include <TFile.h>
00023 #include <TCanvas.h>
00024 #include <TGraph2D.h>
00025 #include <TStyle.h>
00026 #include "TAxis.h"
00027
00028 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00029
00030 typedef HepGeom::Point3D<double> HepPoint3D;
00031 #endif
00032 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00033 typedef HepGeom::Vector3D<double> HepVector3D;
00034 #endif
00035 using namespace CLHEP;
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047 MagFieldReader::MagFieldReader( const std::string& name,
00048 ISvcLocator* pSvcLocator)
00049 : Algorithm ( name , pSvcLocator )
00050 , m_pIMF(0)
00051 {
00052 declareProperty("Zmin", m_zMin = -1200.0*mm);
00053 declareProperty("Zmax", m_zMax = 1200.0*mm);
00054 declareProperty("Step", m_step = 50.0*mm);
00055 declareProperty("Xmin", m_xMin = -900.0*mm);
00056 declareProperty("Xmax", m_xMax = 900.0*mm);
00057 declareProperty("Ymin", m_yMin = -900.0*mm);
00058 declareProperty("Ymax", m_yMax = 900.0*mm);
00059 declareProperty("filename", m_filename ="rawmode3_out.dat");
00060 }
00061
00062
00063
00064
00065 StatusCode MagFieldReader::initialize() {
00066
00067 MsgStream msg(msgSvc(), name());
00068
00069 msg << MSG::INFO << "FieldReader intialize() has been called" << endreq;
00070 StatusCode status = service("MagneticFieldSvc", m_pIMF, true);
00071 if( !status.isSuccess() ) {
00072 msg << MSG::FATAL << "Unable to open Magnetic field service"
00073 << endreq;
00074 return StatusCode::FAILURE;
00075 }
00076
00077 msg << MSG::DEBUG << " Book ntuple" << endreq;
00078
00079 NTuplePtr nt1(ntupleSvc(), "FILE1/n1");
00080 if(nt1) m_tuple1 = nt1;
00081 else {
00082 m_tuple1 = ntupleSvc()->book("FILE1/n1",CLID_ColumnWiseTuple,"Field");
00083 if( m_tuple1 ) {
00084 status = m_tuple1->addItem("x", m_x );
00085 status = m_tuple1->addItem("y", m_y );
00086 status = m_tuple1->addItem("z", m_z );
00087 status = m_tuple1->addItem("r", m_r );
00088 status = m_tuple1->addItem("Bx", m_Bx );
00089 status = m_tuple1->addItem("By", m_By );
00090 status = m_tuple1->addItem("Bz", m_Bz );
00091 status = m_tuple1->addItem("SigmaBx", m_sigma_bx );
00092 status = m_tuple1->addItem("SigmaBy", m_sigma_by );
00093 status = m_tuple1->addItem("SigmaBz", m_sigma_bz );
00094 }
00095 else {
00096 msg << MSG::ERROR << " Cannot book N-tuple:" <<long(m_tuple1)<< endreq;
00097 return StatusCode::FAILURE;
00098 }
00099 }
00100
00101 NTuplePtr nt2(ntupleSvc(), "FILE1/n2");
00102 if(nt2) m_tuple2 = nt2;
00103 else {
00104 m_tuple2 = ntupleSvc()->book("FILE1/n2",CLID_ColumnWiseTuple,"Field");
00105 if( m_tuple2 ) {
00106 status = m_tuple2->addItem("x", m_x2 );
00107 status = m_tuple2->addItem("y", m_y2 );
00108 status = m_tuple2->addItem("z", m_z2 );
00109 status = m_tuple2->addItem("r", m_r2 );
00110 status = m_tuple2->addItem("Bx", m_Bx2 );
00111 status = m_tuple2->addItem("By", m_By2 );
00112 status = m_tuple2->addItem("Bz", m_Bz2 );
00113 }
00114 else {
00115 msg << MSG::ERROR << " Cannot book N-tuple:" <<long(m_tuple2)<< endreq;
00116 return StatusCode::FAILURE;
00117 }
00118 }
00119
00120 NTuplePtr nt3(ntupleSvc(), "FILE1/n3");
00121 if(nt3) m_tuple3 = nt3;
00122 else {
00123 m_tuple3 = ntupleSvc()->book("FILE1/n3",CLID_ColumnWiseTuple,"Field");
00124 if( m_tuple3 ) {
00125 status = m_tuple3->addItem("x", m_x3 );
00126 status = m_tuple3->addItem("y", m_y3 );
00127 status = m_tuple3->addItem("z", m_z3 );
00128 status = m_tuple3->addItem("r", m_r3 );
00129 status = m_tuple3->addItem("phi", m_phi3 );
00130 status = m_tuple3->addItem("Bx", m_Bx3 );
00131 status = m_tuple3->addItem("By", m_By3 );
00132 status = m_tuple3->addItem("Bz", m_Bz3 );
00133 }
00134 else {
00135 msg << MSG::ERROR << " Cannot book N-tuple:" <<long(m_tuple3)<< endreq;
00136 return StatusCode::FAILURE;
00137 }
00138 }
00139
00140 NTuplePtr nt4(ntupleSvc(), "FILE1/n4");
00141 if(nt4) m_tuple4 = nt4;
00142 else {
00143 m_tuple4 = ntupleSvc()->book("FILE1/n4",CLID_ColumnWiseTuple,"Field");
00144 if( m_tuple4 ) {
00145 status = m_tuple4->addItem("time", m_time );
00146 }
00147 else {
00148 msg << MSG::ERROR << " Cannot book N-tuple:" <<long(m_tuple4)<< endreq;
00149 return StatusCode::FAILURE;
00150 }
00151 }
00152
00153 status = service("BesTimerSvc", m_timersvc);
00154 if (status.isFailure()) {
00155 msg << MSG::ERROR << name() << ": Unable to locate BesTimer Service" << endreq;
00156 return StatusCode::FAILURE;
00157 }
00158 m_timer = m_timersvc->addItem("Read field Time");
00159
00160 msg << MSG::INFO << "MagFieldReader initialized" << endreq;
00161 return StatusCode::SUCCESS;
00162 }
00163
00164
00165
00166
00167 StatusCode MagFieldReader::execute() {
00168
00169 MsgStream msg( msgSvc(), name() );
00170
00171 msg << MSG::DEBUG << "==> Execute MagFieldReader" << endreq;
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280 m_timer->start();
00281 for(int i = 0; i < 20000; i++)
00282 {
00283 double px,py,pz,bx,by,bz;
00284 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;
00285 px = min_x + (max_x - min_x)*RandFlat::shoot();
00286 py = min_y + (max_y - min_y)*RandFlat::shoot();
00287 pz = min_z + (max_z - min_z)*RandFlat::shoot();
00288 HepPoint3D r(px,py,pz);
00289 HepVector3D b;
00290 m_pIMF->fieldVector(r,b);
00291 bx = b.x()/tesla;
00292 by = b.y()/tesla;
00293 bz = b.z()/tesla;
00294 }
00295 m_timer->stop();
00296 m_time = m_timer->elapsed();
00297 m_tuple4->write();
00298
00299
00300
00301
00302
00303
00304
00305
00306 msg << MSG::INFO << "MagFieldReader executed" << endreq;
00307 return StatusCode::SUCCESS;
00308 }
00309
00310
00311
00312
00313 StatusCode MagFieldReader::finalize() {
00314
00315 MsgStream msg(msgSvc(), name());
00316 msg << MSG::DEBUG << "==> Finalize MagFieldReader" << endreq;
00317
00318 return StatusCode::SUCCESS;
00319 }
00320
00321
00322
00323
00324
00325
00326 StatusCode MagFieldReader::readField() {
00327
00328 MsgStream msg( msgSvc(), name() );
00329 msg << MSG::DEBUG << "m_pIMF = " << m_pIMF << endreq;
00330
00331 double referfield = m_pIMF->getReferField()/tesla;
00332 msg << MSG::DEBUG << "The reference field is " << referfield << " tesla" <<endreq;
00333 bool ifrealfield = m_pIMF->ifRealField();
00334 if(ifrealfield) msg << MSG::DEBUG << "Use the real field"<<endreq;
00335 else msg << MSG::DEBUG << "Use the fake field"<<endreq;
00336
00337 HepVector3D B(0.0,0.0,0.0);
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364 const int n1=241;
00365 double x[n1],y[n1],y1[n1],y2[n1],y3[n1],y4[n1];
00366
00367 const int n2=121;
00368 double x1[n2],bz1[n2],bz2[n2],bz3[n2],bz4[n2],bz5[n2],bz6[n2];
00369 double br1[n2],br2[n2],br3[n2],br4[n2],br5[n2],br6[n2];
00370 double bp1[n2],bp2[n2],bp3[n2],bp4[n2],bp5[n2],bp6[n2];
00371 double bphi1[n2],bphi2[n2],bphi3[n2],bphi4[n2],bphi5[n2],bphi6[n2];
00372
00373 const int n3=161;
00374 double x2[n3],bx1[n3],bx2[n3],bx3[n3],bx4[n3],bx5[n3],bx6[n3];
00375 double by1[n3],by2[n3],by3[n3],by4[n3],by5[n3],by6[n3];
00376
00377 const int n4=521;
00378 double globle_x[n4],globle_bz[n4];
00379 int i=0;
00380 double px,py,pz;
00381 HepPoint3D pos(px,py,pz);
00382 for( px = -2600; px <= 2600; px += 10 ) {
00383 pos[0] = px; pos[1] = 0; pos[2] = 0;
00384 m_pIMF->fieldVector( pos, B );
00385 globle_x[i]=px/m;
00386 globle_bz[i]=B.z()/tesla;
00387 i++;
00388 }
00389 i=0;
00390 for ( pz = m_zMin; pz <= m_zMax; pz += 10 ) {
00391 pos[0] = 0; pos[1] = 0; pos[2] = pz;
00392 m_pIMF->fieldVector( pos, B );
00393 x[i]=pz/m;
00394 y[i]=B.z()/tesla;
00395
00396 pos[0] = 200; pos[1] = 0; pos[2] = pz;
00397 m_pIMF->fieldVector( pos, B );
00398 y1[i]=B.z()/tesla;
00399
00400 pos[0] = 400; pos[1] = 0; pos[2] = pz;
00401 m_pIMF->fieldVector( pos, B );
00402 y2[i]=B.z()/tesla;
00403
00404 pos[0] = 600; pos[1] = 0; pos[2] = pz;
00405 m_pIMF->fieldVector( pos, B );
00406 y3[i]=B.z()/tesla;
00407
00408 pos[0] = 800; pos[1] = 0; pos[2] = pz;
00409 m_pIMF->fieldVector( pos, B );
00410 y4[i]=B.z()/tesla;
00411 i++;
00412 }
00413 int j = 0;
00414 double tbx,tby,tbz,tbr,tbp,tbphi;
00415 for ( pz = 0; pz <= m_zMax; pz += 10 ) {
00416 pos[0] = 50; pos[1] = 0; pos[2] = pz;
00417 m_pIMF->fieldVector( pos, B );
00418 x1[j]=pz/m;
00419 tbx=B.x()/tesla;
00420 tby=B.y()/tesla;
00421 tbz=B.z()/tesla;
00422 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]);
00423 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]);
00424 tbp=tbz-tbr*pz/sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
00425 bz1[j] = tbz;
00426 br1[j] = tbr*10000;
00427 bp1[j] = tbp;
00428 bphi1[j] = tbphi*10000;
00429
00430 pos[0] = 100; pos[1] = 0; pos[2] = pz;
00431 m_pIMF->fieldVector( pos, B );
00432 tbx=B.x()/tesla;
00433 tby=B.y()/tesla;
00434 tbz=B.z()/tesla;
00435 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]);
00436 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]);
00437 tbp=tbz-tbr*pz/sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
00438 bz2[j] = tbz;
00439 br2[j] = tbr*10000;
00440 bp2[j] = tbp;
00441 bphi2[j] = tbphi*10000;
00442
00443 pos[0] = 200; pos[1] = 0; pos[2] = pz;
00444 m_pIMF->fieldVector( pos, B );
00445 tbx=B.x()/tesla;
00446 tby=B.y()/tesla;
00447 tbz=B.z()/tesla;
00448 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]);
00449 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]);
00450 tbp=tbz-tbr*pz/sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
00451 bz3[j] = tbz;
00452 br3[j] = tbr*10000;
00453 bp3[j] = tbp;
00454 bphi3[j] = tbphi*10000;
00455
00456 pos[0] = 400; pos[1] = 0; pos[2] = pz;
00457 m_pIMF->fieldVector( pos, B );
00458 tbx=B.x()/tesla;
00459 tby=B.y()/tesla;
00460 tbz=B.z()/tesla;
00461 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]);
00462 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]);
00463 tbp=tbz-tbr*pz/sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
00464 bz4[j] = tbz;
00465 br4[j] = tbr*10000;
00466 bp4[j] = tbp;
00467 bphi4[j] = tbphi*10000;
00468
00469 pos[0] = 600; pos[1] = 0; pos[2] = pz;
00470 m_pIMF->fieldVector( pos, B );
00471 tbx=B.x()/tesla;
00472 tby=B.y()/tesla;
00473 tbz=B.z()/tesla;
00474 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]);
00475 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]);
00476 tbp=tbz-tbr*pz/sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
00477 bz5[j] = tbz;
00478 br5[j] = tbr*10000;
00479 bp5[j] = tbp;
00480 bphi5[j] = tbphi*10000;
00481
00482 pos[0] = 800; pos[1] = 0; pos[2] = pz;
00483 m_pIMF->fieldVector( pos, B );
00484 tbx=B.x()/tesla;
00485 tby=B.y()/tesla;
00486 tbz=B.z()/tesla;
00487 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]);
00488 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]);
00489 tbp=tbz-tbr*pz/sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
00490 bz6[j] = tbz;
00491 br6[j] = tbr*10000;
00492 bp6[j] = tbp;
00493 bphi6[j] = tbphi*10000;
00494 j++;
00495 }
00496 int l = 0;
00497 for(py = m_yMin; py<= m_yMax; py +=10){
00498 pos[0] = 0; pos[1] = py; pos[2] = 0;
00499 m_pIMF->fieldVector( pos, B );
00500 tbx=B.x()/tesla*10000;
00501 tby=B.y()/tesla*10000;
00502 tbz=B.z()/tesla*10000;
00503 x2[l] = py/m;
00504 bx1[l] = tbx;
00505 by1[l] = tby;
00506
00507 pos[0] = 100; pos[1] = py; pos[2] = 100;
00508 m_pIMF->fieldVector( pos, B );
00509 tbx=B.x()/tesla*10000;
00510 tby=B.y()/tesla*10000;
00511 tbz=B.z()/tesla*10000;
00512 bx2[l] = tbx;
00513 by2[l] = tby;
00514
00515 pos[0] = 400; pos[1] = py; pos[2] = 400;
00516 m_pIMF->fieldVector( pos, B );
00517 tbx=B.x()/tesla*10000;
00518 tby=B.y()/tesla*10000;
00519 tbz=B.z()/tesla*10000;
00520 bx3[l] = tbx;
00521 by3[l] = tby;
00522 l++;
00523 }
00524
00525 TGraph2D *dt = new TGraph2D();
00526 int k = 0;
00527 for ( pz = -3000; pz <= 3000; pz += 50 ) {
00528 for( py = -2700; py <= 2700; py += 50){
00529 pos[0] = 0; pos[1] = py; pos[2] = pz;
00530 m_pIMF->fieldVector( pos, B );
00531 tbx=B.x()/tesla;
00532 tby=B.y()/tesla;
00533 tbz=B.z()/tesla;
00534 dt->SetPoint(k,pz/m,py/m,tbz);
00535 k++;
00536 }
00537 }
00538 TGraph2D *dt1 = new TGraph2D();
00539 k = 0;
00540 for ( pz = -3000; pz <= 3000; pz += 50 ) {
00541 for( py = -3000; py <= 3000; py += 50){
00542 pos[0] = 0; pos[1] = py; pos[2] = pz;
00543 m_pIMF->fieldVector( pos, B );
00544 tbx=B.x()/tesla;
00545 tby=B.y()/tesla;
00546 tbz=B.z()/tesla;
00547 double btot = sqrt(tbx*tbx+tby*tby+tbz*tbz);
00548 dt1->SetPoint(k,pz/m,py/m,btot);
00549 k++;
00550 }
00551 }
00552 TGraph2D *dt2 = new TGraph2D();
00553 k = 0;
00554 for ( pz = m_zMin; pz <= m_zMax; pz += 50 ) {
00555 for( py = m_yMin; py <= m_yMax; py += 50){
00556 pos[0] = 0; pos[1] = py; pos[2] = pz;
00557 m_pIMF->fieldVector( pos, B );
00558 tbx=B.x()/tesla;
00559 tby=B.y()/tesla;
00560 tbz=B.z()/tesla;
00561
00562 dt2->SetPoint(k,pz/m,py/m,tbz);
00563 k++;
00564 }
00565 }
00566 gStyle->SetPalette(1);
00567 gStyle->SetOptTitle(1);
00568 gStyle->SetOptStat(0);
00569
00570 TFile* f1=new TFile("graph.root","RECREATE");
00571 TGraph *gr1 = new TGraph(n1,x,y);
00572 TGraph *gr2 = new TGraph(n1,x,y1);
00573 TGraph *gr3 = new TGraph(n1,x,y2);
00574 TGraph *gr4 = new TGraph(n1,x,y3);
00575 TGraph *gr5 = new TGraph(n1,x,y4);
00576
00577 TGraph *g1 = new TGraph(n2,x1,bz1);
00578 TGraph *g2 = new TGraph(n2,x1,bz2);
00579 TGraph *g3 = new TGraph(n2,x1,bz3);
00580 TGraph *g4 = new TGraph(n2,x1,bz4);
00581 TGraph *g5 = new TGraph(n2,x1,bz5);
00582 TGraph *g6 = new TGraph(n2,x1,bz6);
00583
00584 TGraph *g7 = new TGraph(n2,x1,br1);
00585 TGraph *g8 = new TGraph(n2,x1,br2);
00586 TGraph *g9 = new TGraph(n2,x1,br3);
00587 TGraph *g10 = new TGraph(n2,x1,br4);
00588 TGraph *g11 = new TGraph(n2,x1,br5);
00589 TGraph *g12 = new TGraph(n2,x1,br6);
00590
00591 TGraph *g13 = new TGraph(n2,x1,bp1);
00592 TGraph *g14 = new TGraph(n2,x1,bp2);
00593 TGraph *g15 = new TGraph(n2,x1,bp3);
00594 TGraph *g16 = new TGraph(n2,x1,bp4);
00595 TGraph *g17 = new TGraph(n2,x1,bp5);
00596 TGraph *g18 = new TGraph(n2,x1,bp6);
00597
00598 TGraph *g19 = new TGraph(n2,x1,bphi1);
00599 TGraph *g20 = new TGraph(n2,x1,bphi2);
00600 TGraph *g21 = new TGraph(n2,x1,bphi3);
00601 TGraph *g22 = new TGraph(n2,x1,bphi4);
00602 TGraph *g23 = new TGraph(n2,x1,bphi5);
00603 TGraph *g24 = new TGraph(n2,x1,bphi6);
00604
00605 TGraph *g25 = new TGraph(n3,x2,bx1);
00606 TGraph *g26 = new TGraph(n3,x2,bx2);
00607 TGraph *g27 = new TGraph(n3,x2,bx3);
00608
00609 TGraph *g28 = new TGraph(n3,x2,by1);
00610 TGraph *g29 = new TGraph(n3,x2,by2);
00611 TGraph *g30 = new TGraph(n3,x2,by3);
00612
00613 TGraph *g31 = new TGraph(n4,globle_x,globle_bz);
00614
00615 TCanvas *c1 = new TCanvas("c1","Two Graphs",200,10,600,400);
00616 TCanvas *c2 = new TCanvas("c2","Two Graphs",200,10,600,400);
00617 c1->Divide(2,2);
00618 c2->Divide(2,2);
00619 c1->cd(1);
00620 gr1->SetLineColor(2);
00621 gr1->SetLineWidth(2);
00622 gr1->Draw("ALP");
00623 gr1->SetTitle("bz vs z (y=0,x=0,200,400,600,800mm)");
00624 gr1->GetXaxis()->SetTitle("m");
00625 gr1->GetYaxis()->SetTitle("tesla");
00626 gr2->SetLineWidth(2);
00627 gr2->SetLineColor(3);
00628 gr2->Draw("LP");
00629 gr3->SetLineWidth(2);
00630 gr3->SetLineColor(4);
00631 gr3->Draw("LP");
00632 gr4->SetLineWidth(2);
00633 gr4->SetLineColor(5);
00634 gr4->Draw("LP");
00635 gr5->SetLineWidth(2);
00636 gr5->SetLineColor(6);
00637 gr5->Draw("LP");
00638
00639 c1->cd(2);
00640 g1->SetLineColor(2);
00641 g1->SetLineWidth(2);
00642 g1->Draw("ALP");
00643 g1->SetTitle("bz(red),bendp vs z (y=0,x=50,100,200,400,600,800mm)");
00644 g1->GetXaxis()->SetTitle("m");
00645 g1->GetYaxis()->SetTitle("tesla");
00646 g1->GetYaxis()->SetRangeUser(0.2,2);
00647 g2->SetLineWidth(2);
00648 g2->SetLineColor(2);
00649 g2->Draw("LP");
00650 g3->SetLineWidth(2);
00651 g3->SetLineColor(2);
00652 g3->Draw("LP");
00653 g4->SetLineWidth(2);
00654 g4->SetLineColor(2);
00655 g4->Draw("LP");
00656 g5->SetLineWidth(2);
00657 g5->SetLineColor(2);
00658 g5->Draw("LP");
00659 g6->SetLineColor(2);
00660 g6->SetLineWidth(2);
00661 g6->Draw("LP");
00662
00663 g13->SetLineWidth(2);
00664 g13->SetLineColor(4);
00665 g13->Draw("LP");
00666 g14->SetLineWidth(2);
00667 g14->SetLineColor(4);
00668 g14->Draw("LP");
00669 g15->SetLineWidth(2);
00670 g15->SetLineColor(4);
00671 g15->Draw("LP");
00672 g16->SetLineWidth(2);
00673 g16->SetLineColor(4);
00674 g16->Draw("LP");
00675 g17->SetLineWidth(2);
00676 g17->SetLineColor(4);
00677 g17->Draw("LP");
00678 g18->SetLineWidth(2);
00679 g18->SetLineColor(4);
00680 g18->Draw("LP");
00681
00682 c1->cd(3);
00683 g7->SetLineWidth(2);
00684 g7->SetLineColor(2);
00685 g7->Draw("ALP");
00686 g7->SetTitle("br vs z (y=0,x=50,100,200,400,600,800mm)");
00687 g7->GetXaxis()->SetTitle("m");
00688 g7->GetYaxis()->SetTitle("gauss");
00689 g7->GetYaxis()->SetRangeUser(-1100,1000);
00690 g8->SetLineWidth(2);
00691 g8->SetLineColor(3);
00692 g8->Draw("LP");
00693 g9->SetLineWidth(2);
00694 g9->SetLineColor(4);
00695 g9->Draw("LP");
00696 g10->SetLineWidth(2);
00697 g10->SetLineColor(5);
00698 g10->Draw("LP");
00699 g11->SetLineWidth(2);
00700 g11->SetLineColor(6);
00701 g11->Draw("LP");
00702 g12->SetLineWidth(2);
00703 g12->SetLineColor(7);
00704 g12->Draw("LP");
00705
00706 c1->cd(4);
00707 g19->SetLineWidth(2);
00708 g19->SetLineColor(2);
00709 g19->Draw("ALP");
00710 g19->SetTitle("bphi vs z (y=0,x=50,100,200,400,600,800mm)");
00711 g19->GetXaxis()->SetTitle("m");
00712 g19->GetYaxis()->SetTitle("gauss");
00713 g19->GetYaxis()->SetRangeUser(-1000,200);
00714 g20->SetLineWidth(2);
00715 g20->SetLineColor(3);
00716 g20->Draw("LP");
00717 g21->SetLineWidth(2);
00718 g21->SetLineColor(4);
00719 g21->Draw("LP");
00720 g22->SetLineWidth(2);
00721 g22->SetLineColor(5);
00722 g22->Draw("LP");
00723 g23->SetLineWidth(2);
00724 g23->SetLineColor(6);
00725 g23->Draw("LP");
00726 g24->SetLineWidth(2);
00727 g24->SetLineColor(7);
00728 g24->Draw("LP");
00729
00730 TCanvas *c3 = new TCanvas("c3","Two Graphs",200,10,600,400);
00731 c3->cd(1);
00732
00733
00734
00735 dt->Draw("z sinusoidal");
00736 dt->SetTitle("bz vs y,z (x=0)");
00737
00738 TCanvas *c4 = new TCanvas("c4","Two Graphs",200,10,600,400);
00739 c4->cd(1);
00740 dt1->Draw("z sinusoidal");
00741 dt1->SetTitle("btot vs y,z (x=0)");
00742
00743 TCanvas *c5 = new TCanvas("c5","Two Graphs",200,10,600,400);
00744 c5->cd(1);
00745 dt2->Draw("z sinusoidal");
00746 dt2->SetTitle("bz vs y,z (x=0)");
00747
00748 c2->cd(2);
00749 g25->SetLineWidth(2);
00750 g25->SetLineColor(2);
00751 g25->Draw("ALP");
00752 g25->SetTitle("bx(red),by vs y (x,z=0,100,400mm)");
00753 g25->GetXaxis()->SetTitle("m");
00754 g25->GetYaxis()->SetTitle("gauss");
00755 g25->GetYaxis()->SetRangeUser(-200,300);
00756 g26->SetLineWidth(2);
00757 g26->SetLineColor(2);
00758 g26->Draw("LP");
00759 g27->SetLineWidth(2);
00760 g27->SetLineColor(2);
00761 g27->Draw("LP");
00762
00763 g28->SetLineWidth(2);
00764 g28->SetLineColor(3);
00765 g28->Draw("LP");
00766 g29->SetLineWidth(2);
00767 g29->SetLineColor(3);
00768 g29->Draw("LP");
00769 g30->SetLineWidth(2);
00770 g30->SetLineColor(3);
00771 g30->Draw("LP");
00772
00773 c2->cd(3);
00774 g31->SetLineWidth(2);
00775 g31->SetLineColor(2);
00776 g31->Draw("ALP");
00777 g31->SetTitle("bz vs x (y,z=0)");
00778 g31->GetXaxis()->SetTitle("m");
00779 g31->GetYaxis()->SetTitle("tesla");
00780
00781 c1->Write();
00782 c2->Write();
00783 c3->Write();
00784 c4->Write();
00785 c5->Write();
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807 return StatusCode::SUCCESS;
00808 }