MagFieldReader Class Reference

#include <MagFieldReader.h>

Inheritance diagram for MagFieldReader:

Public Member Functions
	MagFieldReader (const std::string &name, ISvcLocator *pSvcLocator)
	Standard constructor.
virtual	~MagFieldReader ()
virtual StatusCode	initialize ()
	Destructor.
virtual StatusCode	execute ()
	Algorithm execution.
virtual StatusCode	finalize ()
	Algorithm finalization.

Detailed Description

Parameters

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.

Definition at line 23 of file MagFieldReader.h.

Constructor & Destructor Documentation

MagFieldReader()

MagFieldReader::MagFieldReader	(	const std::string &	name,
		ISvcLocator *	pSvcLocator )

Standard constructor.

Definition at line 44 of file MagFieldReader.cxx.

    : Algorithm( name, pSvcLocator ), m_pIMF( 0 ) {
  declareProperty( "Zmin", m_zMin = -1200.0 * mm );
  declareProperty( "Zmax", m_zMax = 1200.0 * mm );
  declareProperty( "Step", m_step = 50.0 * mm );
  declareProperty( "Xmin", m_xMin = -900.0 * mm );
  declareProperty( "Xmax", m_xMax = 900.0 * mm );
  declareProperty( "Ymin", m_yMin = -900.0 * mm );
  declareProperty( "Ymax", m_yMax = 900.0 * mm );
  declareProperty( "filename", m_filename = "rawmode3_out.dat" );
}

~MagFieldReader()

virtual MagFieldReader::~MagFieldReader ( )

inlinevirtual

Definition at line 28 of file MagFieldReader.h.

{}; ///< Destructor

Member Function Documentation

execute()

StatusCode MagFieldReader::execute ( )

virtual

Algorithm execution.

Definition at line 171 of file MagFieldReader.cxx.

                                   {
 
  MsgStream msg( msgSvc(), name() );
 
  msg << MSG::DEBUG << "==> Execute MagFieldReader" << endmsg;
  /*
    //calculate the error of Bfield
    std::ifstream infile( m_filename.c_str() );
    if(!infile.good()) msg << MSG::FATAL << "Unable to read the file: " << m_filename <<
    endmsg; int nLine = 0; char line[ 255 ]; while( infile ) {
      // parse each line of the file,
      // comment lines begin with '#' in the cdf file
      infile.getline( line, 255 );
      if ( line[0] == '#' ) continue;
      std::string  gridx, gridy, gridz, sFx, sFy, sFz;
      char* token = strtok( line, " " );
      if ( token ) { gridx = token; token = strtok( NULL, " " );} else continue;
      if ( token ) { gridy = token; token = strtok( NULL, " " );} else continue;
      if ( token ) { gridz = token; token = strtok( NULL, " " );} else continue;
      if ( token ) { sFx = token; token = strtok( NULL, " " );} else continue;
      if ( token ) { sFy = token; token = strtok( NULL, " " );} else continue;
      if ( token ) { sFz = token; token = strtok( NULL, " " );} else continue;
      if ( token != NULL ) continue;
      nLine++;
      double x,y,z,bx,by,bz;
      //Grid position
      x = atof( gridx.c_str() );
      y = atof( gridy.c_str() );
      z = atof( gridz.c_str() );
      // Field values are given in tesla in CDF file. Convert to CLHEP units
      bx = atof( sFx.c_str() );
      by = atof( sFy.c_str() );
      bz = atof( sFz.c_str() );
      //std::cout<<"x,y,z: "<<x<<" "<<y<<" "<<z<<" "<<"bx,by,bz: "<<bx<<" "<<by<<"
    "<<bz<<std::endl; x = -x*m; y = y*m; z = -z*m;
 
      HepPoint3D r(x,y,z);
      HepVector3D b;
      m_pIMF->fieldVector(r,b);
      m_Bx = b.x()/tesla;
      m_By = b.y()/tesla;
      m_Bz = b.z()/tesla;
      m_sigma_bx = (b.x()/tesla + bx)*10000.;
      m_sigma_by = (b.y()/tesla - by)*10000.;
      m_sigma_bz = (b.z()/tesla + bz)*10000.;
      m_r = std::sqrt(x/m*x/m+y/m*y/m);
      m_x = x/m;
      m_y = y/m;
      m_z = z/m;
      m_tuple1->write();
    }
    infile.close();
    std::cout<<"Totally "<<nLine<<" in file "<<m_filename<<std::endl;
 
    for(int k = 0; k < 5; k++) {
      double rr;
      if(k==0) rr = 0;
      if(k==1) rr = 200;
      if(k==2) rr = 400;
      if(k==3) rr = 600;
      if(k==4) rr = 800;
      for(int zz = -1195; zz <1200; zz+=50) {
        double bxt = 0.,byt= 0.,bzt = 0.;
        for(int i = 0; i < 100000; i++) {
          double phi = CLHEP::twopi*CLHEP::RandFlat::shoot();
          double xx = rr*std::cos(phi);
          double yy = rr*std::sin(phi);
          HepPoint3D r(xx,yy,double(zz));
          HepVector3D b;
          m_pIMF->fieldVector(r,b);
          bxt+=b.x()/tesla;
          byt+=b.y()/tesla;
          bzt+=b.z()/tesla;
        }
        m_z2 = zz;
        m_r2 = rr;
        m_Bx2 = bxt/100000;
        m_By2 = byt/100000;
        m_Bz2 = bzt/100000;
        m_tuple2->write();
      }
    }
 
    for(int k = 0; k < 3; k++) {
      double zz;
      if(k==0) zz = 0;
      if(k==1) zz = -800;
      if(k==2) zz = 800;
      for(int rr = 200; rr <=600; rr+=400) {
        for(double phi = 0; phi <= 2*3.14159; phi+=0.1745) {
          //double phi = CLHEP::twopi*CLHEP::RandFlat::shoot();
          double xx = rr*std::cos(phi);
          double yy = rr*std::sin(phi);
          HepPoint3D r(xx,yy,double(zz));
          HepVector3D b;
          m_pIMF->fieldVector(r,b);
          m_Bx3 = b.x()/tesla;
          m_By3 = b.y()/tesla;
          m_Bz3 = b.z()/tesla;
          m_phi3 = phi*180/CLHEP::pi;
          m_x3 = xx;
          m_y3 = yy;
          m_z3 = zz;
          m_r3 = rr;
          m_tuple3->write();
        }
      }
    }
  */
  m_timer->start();
  for ( int i = 0; i < 20000; i++ )
  {
    double px, py, pz, bx, by, bz;
    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;
    px = min_x + ( max_x - min_x ) * RandFlat::shoot();
    py = min_y + ( max_y - min_y ) * RandFlat::shoot();
    pz = min_z + ( max_z - min_z ) * RandFlat::shoot();
    HepPoint3D  r( px, py, pz );
    HepVector3D b;
    m_pIMF->fieldVector( r, b );
    bx = b.x() / tesla;
    by = b.y() / tesla;
    bz = b.z() / tesla;
  }
  m_timer->stop();
  m_time = m_timer->elapsed();
  m_tuple4->write();
  /*
    StatusCode sc = this->readField();
    if( sc.isFailure() ) {
      msg << MSG::FATAL << "Unable to execute MagFieldReader"
          << endmsg;
      return sc;
    }
  */
  msg << MSG::INFO << "MagFieldReader executed" << endmsg;
  return StatusCode::SUCCESS;
}

finalize()

StatusCode MagFieldReader::finalize ( )

virtual

Algorithm finalization.

Definition at line 314 of file MagFieldReader.cxx.

                                    {
 
  MsgStream msg( msgSvc(), name() );
  msg << MSG::DEBUG << "==> Finalize MagFieldReader" << endmsg;
 
  return StatusCode::SUCCESS;
}

initialize()

StatusCode MagFieldReader::initialize ( )

virtual

Destructor.

Algorithm initialization

Definition at line 59 of file MagFieldReader.cxx.

                                      {
 
  MsgStream msg( msgSvc(), name() );
 
  msg << MSG::INFO << "FieldReader intialize() has been called" << endmsg;
  StatusCode status = service( "MagneticFieldSvc", m_pIMF, true );
  if ( !status.isSuccess() )
  {
    msg << MSG::FATAL << "Unable to open Magnetic field service" << endmsg;
    return StatusCode::FAILURE;
  }
 
  msg << MSG::DEBUG << " Book ntuple" << endmsg;
 
  NTuplePtr nt1( ntupleSvc(), "FILE1/n1" );
  if ( nt1 ) m_tuple1 = nt1;
  else
  {
    m_tuple1 = ntupleSvc()->book( "FILE1/n1", CLID_ColumnWiseTuple, "Field" );
    if ( m_tuple1 )
    {
      status = m_tuple1->addItem( "x", m_x );
      status = m_tuple1->addItem( "y", m_y );
      status = m_tuple1->addItem( "z", m_z );
      status = m_tuple1->addItem( "r", m_r );
      status = m_tuple1->addItem( "Bx", m_Bx );
      status = m_tuple1->addItem( "By", m_By );
      status = m_tuple1->addItem( "Bz", m_Bz );
      status = m_tuple1->addItem( "SigmaBx", m_sigma_bx );
      status = m_tuple1->addItem( "SigmaBy", m_sigma_by );
      status = m_tuple1->addItem( "SigmaBz", m_sigma_bz );
    }
    else
    {
      msg << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple1 ) << endmsg;
      return StatusCode::FAILURE;
    }
  }
 
  NTuplePtr nt2( ntupleSvc(), "FILE1/n2" );
  if ( nt2 ) m_tuple2 = nt2;
  else
  {
    m_tuple2 = ntupleSvc()->book( "FILE1/n2", CLID_ColumnWiseTuple, "Field" );
    if ( m_tuple2 )
    {
      status = m_tuple2->addItem( "x", m_x2 );
      status = m_tuple2->addItem( "y", m_y2 );
      status = m_tuple2->addItem( "z", m_z2 );
      status = m_tuple2->addItem( "r", m_r2 );
      status = m_tuple2->addItem( "Bx", m_Bx2 );
      status = m_tuple2->addItem( "By", m_By2 );
      status = m_tuple2->addItem( "Bz", m_Bz2 );
    }
    else
    {
      msg << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple2 ) << endmsg;
      return StatusCode::FAILURE;
    }
  }
 
  NTuplePtr nt3( ntupleSvc(), "FILE1/n3" );
  if ( nt3 ) m_tuple3 = nt3;
  else
  {
    m_tuple3 = ntupleSvc()->book( "FILE1/n3", CLID_ColumnWiseTuple, "Field" );
    if ( m_tuple3 )
    {
      status = m_tuple3->addItem( "x", m_x3 );
      status = m_tuple3->addItem( "y", m_y3 );
      status = m_tuple3->addItem( "z", m_z3 );
      status = m_tuple3->addItem( "r", m_r3 );
      status = m_tuple3->addItem( "phi", m_phi3 );
      status = m_tuple3->addItem( "Bx", m_Bx3 );
      status = m_tuple3->addItem( "By", m_By3 );
      status = m_tuple3->addItem( "Bz", m_Bz3 );
    }
    else
    {
      msg << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple3 ) << endmsg;
      return StatusCode::FAILURE;
    }
  }
 
  NTuplePtr nt4( ntupleSvc(), "FILE1/n4" );
  if ( nt4 ) m_tuple4 = nt4;
  else
  {
    m_tuple4 = ntupleSvc()->book( "FILE1/n4", CLID_ColumnWiseTuple, "Field" );
    if ( m_tuple4 ) { status = m_tuple4->addItem( "time", m_time ); }
    else
    {
      msg << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple4 ) << endmsg;
      return StatusCode::FAILURE;
    }
  }
 
  status = service( "BesTimerSvc", m_timersvc );
  if ( status.isFailure() )
  {
    msg << MSG::ERROR << name() << ": Unable to locate BesTimer Service" << endmsg;
    return StatusCode::FAILURE;
  }
  m_timer = m_timersvc->addItem( "Read field Time" );
 
  msg << MSG::INFO << "MagFieldReader initialized" << endmsg;
  return StatusCode::SUCCESS;
}

---

The documentation for this class was generated from the following files:

• MagFieldReader.h
• MagFieldReader.cxx