---

Medium.cxx

Go to the documentation of this file.

// $Id: Medium.cxx,v 1.4 2001/09/04 19:26:01 atwood Exp $
//
//
// This file is part of Gismo 2
//
//

#ifdef __GNUG__
#pragma implementation
#endif

#include "gismo/Medium.h"

#include "gismo/Detector.h"
#include "gismo/Material.h"

#include "gismo/Field.h"
#include "gismo/MCParticle.h"
#include "facilities/error.h"

#include "geometry/Box.h"
#include "geometry/Ray.h"

#include "geomrep/BoxRep.h"
#include "geomrep/TubeRep.h"

#include <typeinfo>
#include <cassert>

//              constructors

Medium::Medium(Medium * prnt, float size)
   : _volume(new Box(size,size,size))
   , _material(0)
   , _detector(0)
   , _field(0)
   , _parent(prnt)
   , _title(0)
{
     set_defaults();
}
Medium::Medium(Medium* parent, Shape* vol, const char* matName, Detector* det)
   :  _volume(vol)
   ,  _material(Material::hatch(matName))
   ,  _detector(det)
   ,  _field(0)
   ,  _parent(parent)
   ,  _title(0)
{
        set_defaults();
}

Medium::Medium(const Medium & old)
: _volume(&old._volume->copy())
, _material(old._material)
, _detector(new Detector(*old._detector))
, _field(old._field)
, _parent(old._parent)
, _title(0)
{
    set_defaults();
}
void Medium::set_defaults()
{

   if(_parent) {
     // parent was specified: set appropiate attributes to be identical, then add
      _keCutOff = _parent->kECutOff();
      _maxStep = _parent->maxStep();
      if(!_material) _material = &_parent->material();
      if(!_field) _field = &(_parent->field());
      _parent->addMedium(this);
   }
   else {
        // standard defaults for tracking attributes
      _keCutOff   = 0.1f;
      _maxStep   = 100.;
      if(!_material) _material = Material::hatch("vacuum");
   }
   s_count ++;
}


Medium& Medium::setTitle(const char* newTitle)
{
   if( _title) delete [] _title;
   strcpy(_title = new char[strlen(newTitle)+1], newTitle);
   return *this;
}
Medium::~Medium()
{
   if(_title) delete [] _title;
   if(_volume) delete _volume;
   if(_detector) delete _detector;
    s_count--;
}

//               set tracking attributes: for this and all inner media

Medium& Medium::setKECutOff(float keCut) {
  _keCutOff = keCut;
  return *this;
}

Medium& Medium::setMaxStep(float mxStep) {
  _maxStep = mxStep;
  return *this;
}

Medium& Medium::setField(Field* fld)
{
  if(_parent)
     if(&_parent->field() != _field) delete _field;
  _field = fld;
  return *this;
}

const char* Medium::nameOf() const {
  return "Medium";
}

int Medium::isComposite() const {
  return 0;
}

//                manage inner media

Medium& Medium::addMedium ( Medium* nextMedium)
{
    FATAL("Non-composite Medium tried to add medium");
    return *nextMedium;
}

Medium& Medium::removeMedium (Medium* oldMedium)
{
    FATAL("Non-composite Medium tried to remove medium");
    return *oldMedium;
}

//                       coordinate transformation
GeomObject&
Medium::transform(const CoordTransform& T)
{

    if( _volume )
        _volume->transform(T); // could happen during initialzation
   if( _detector )
     _detector->transform(T);  // note that detector may ignore this!
   return *this;
}

//              Detector  processing methods
void Medium::clear(){
                                if(_detector) _detector->clear();}
void Medium::generateResponse(){
                                if(_detector) _detector->generateResponse();}
void Medium::accept(DetectorVisitor& a){
                                if(_detector) _detector->accept(a);}
void Medium::writeData(std::ostream & os){
                                if(_detector) _detector->writeData(os);}
void Medium::readData(std::istream & is){
                                if(_detector) _detector->readData(is); }
void Medium::printResponse(std::ostream& os)const{
                                if(_detector) _detector->printOn(os);}
void Medium::createResponseView(gui::DisplayRep& v) {
                                if( _detector) _detector->draw(v);}

void Medium::notify() {
                                if( _detector) _detector->notify();}

//                    Methods used in Tracking

const Medium * Medium::inside(const Point& r)const
{
     const Medium* which = this;
     while (which &&  which->_volume && ! which->volume().inside(r))
        which = which->parent();
     return which;

}

Ray * Medium::CreateRay(const Point& r, const Vector& p, float q,float maxStep )const
{
 /*
  * create a ray to extrapolate from point r, for a particle with
  * momentum p and charge q
  *
  * Last arg is maximum step anticipated, to allow eventual Runge-Kutta field
  * if zero, the particle is at rest,and the ray is only to record the
  * position
  */
   if( maxStep ==0 ) return new Ray(r,Vector(0,0,1));
   if( _field )
      return _field->CreateRay(r,p,q);
   else
        return new Ray(r, p);
}


double Medium::distanceToLeave( const Ray& r,
                               const Medium* & newStuff,
                                         double maxS ) const
{
    assert(maxS>=0);

    //  Get the distance to leave our  Shape.
    //  If there is no intersection (looper?) the distance will be a large number

    if( !_volume ) {
        // no associated volume--means that this is the world?
        newStuff = this;
        return maxS;
    }

    double distance = volume().distanceToLeave( r,   maxS );
    newStuff = distance>=maxS? this : parent();

    if( distance == FLT_MAX) {
        // No solution: either trapped inside this Medium (looper) or not in Medium
        // so, check for geometry error
        const Medium *testStuff = inside(r.position());
        if( this != testStuff)  {
            std::ostrstream cerr;
            cerr << "Medium::distanceToLeave: geometry error at " << r.position()
              << "\n\ttrying to leave \"" << title() << "\""; ;
            if( testStuff )
                cerr << "\tbut actually inside \"" <<  testStuff->title()  << "\"\n";
            else
                cerr << " but beyond known world?\n";
            cerr << '\0'; WARNING(cerr.str());
            return distance;  // force lost here
        }
    }

    // Check solution against max allowed step size for this Medium
    if(distance > maxStep())
    {  distance = maxStep();
       newStuff = this;   // if cut back, must be still here
    }
    return distance;
}

double Medium::distanceToEnter( const Ray& r,
                               const Medium* & newStuff,
                               double maxDist ) const
{
    // called by a parent CompoundMedium
    //  If there is no intersection the distance will return FLT_MAX
    //  If the Ray origin is already inside the Medium, then returns zero

    newStuff = this;
    double distance = volume().distanceToEnter( r,  maxDist );

    if(distance == FLT_MAX) {
        if( this != lastMedium &&  volume().inside(r.position(Volume::Surface_EPSILON)))
            distance = 0.;
        else
            newStuff = parent();
                
    }
    return distance;
}

//                 Print and display
//
void Medium::printOn( std::ostream& os ) const
{
    os << "  Medium subClass:" << nameOf() <<" - title: "<< title() <<"\n"
       << "             maxStep = " << maxStep()
       << " keCutOff = "<< kECutOff() << " \n";
    volume().printOn(os);
    if(_field)
            _field->printOn(os);
    os << " Material Name: "
       << _material->name() << '\n';

}

void Medium::createDetectorView(gui::DisplayRep& v) {
    const type_info& t = typeid(volume());

    if( t==typeid(Box) )    v.append(BoxRep(volume() ));
    else if (t==typeid(Tube) )v.append(TubeRep(volume() ));
    else {
#if 0 // disable if want other volumes (like LayeredMedium::Slice) to be ignored
        WARNING_MACRO("no display representation was set for volume " << t.name() );
#endif
    }
}

//-----------------------------------------------------------------
//                  gobal statics
const Medium * Medium::lastMedium = 0;
unsigned Medium::s_count=0;

---