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