---

RCParticle.cxx

Go to the documentation of this file.

//  $Id: RCParticle.cxx,v 1.1 2001/09/22 16:41:38 burnett Exp $
//
// This file is part of Gismo 2
//
//
// Implemention of all RCParticle member functions

#include "gismo/RCParticle.h"

#include "gismo/Medium.h"
#include "gismo/Material.h"
#include "gismo/Detector.h"
#include "gismo/Units.h"
#include "gismo/PDT.h"
#include "gismo/World.h"

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

#include "facilities/error.h"
#include <cassert>

/********************************************************************
 *                RCParticle                                        *
 ********************************************************************/
#ifdef DUMP_STEP // special debug output
# include <fstream>
 static std::ofstream fout("dump_step.prn");
#endif

const Medium* RCParticle::s_world=0;

RCParticle::RCParticle(Point x0, Vector t0, float dist = 0)
        : GParticle("e-")
        , m_maxArcLen(dist)
{   
    setPosition(x0);
    setMomentum(t0);

    assert(s_world != 0);

    m_currentMedium = s_world->inside(x0);
    ini();
    if(dist > 0) { propagate();}
}

RCParticle::~RCParticle()
{ 
}

// Initializes RCParticle for propagation
void RCParticle::ini()
{
  m_arcLength = m_properTime = m_radLength = 0.;
  m_stepcount = 0;
  Medium::setLastMedium(0);
  m_lastMedium = m_currentMedium;
  while (m_stepList.size() > 0) {
      m_stepList.pop_back();
  }
}

//  RCParticle::propagate -- propagate the particle for a distance 
//              specified by maxArcLen or until it leaves the World
//
void RCParticle::propagate() {

// Set max-step here incase a continuation of previous projection
  m_maxStep = m_maxArcLen+.0001;
  setStatus( ALIVE );
// Loop until we're out.... 
  while (status()==ALIVE)
  {
    const Medium * nextMedium=m_currentMedium;  // a temporary for new medium

        // return a Ray that estimates the trajectory (or a point)

        Ray *m_segment=
            m_currentMedium->CreateRay( position(), momentum(), charge(), m_maxStep );

        // flag this segment
        m_segment->setFlag(m_stepcount);

        // ask the Medium for a distance along that ray, (unless
        // the particle is at rest),and for possible detector
        m_step = m_currentMedium->distanceToLeave(*m_segment, nextMedium, m_maxStep);
        if (m_step < m_maxStep) {  // hit a boundary
            setStatus( ALIVE );
            //NEW CODE >>>>>> Cross the multiple shared surfaces..... 
            while(m_step < Volume::Surface_EPSILON && status()==ALIVE) {
                if(m_currentMedium == nextMedium) setStatus (STUCK);
                Medium::setLastMedium(m_currentMedium);
                m_lastMedium = m_currentMedium;
                m_currentMedium = nextMedium;
                m_step = m_currentMedium->distanceToLeave(*m_segment, nextMedium, m_maxStep);
            }
        }       
        if (m_step == FLT_MAX || (m_step<=0 && m_currentMedium==nextMedium) )   {
            WARNING( "RCParticle::propagate -- we are lost or stuck!");
            setStatus( STUCK );
            break;
        }
        else if( m_step > m_maxStep) {
            // maxStep was ignored above: assume still in same medium
            m_step = m_maxStep;
            nextMedium = m_currentMedium;
        }
        
        // increment along trajectory by arclength step
        if(m_step > 0) {
            stepBy(m_step, *m_segment );
        }
    Medium::setLastMedium(m_currentMedium);
    m_lastMedium = m_currentMedium;
    delete m_segment; 

    // We've reached the boundary: update current medium and see if we've left
    if ( status()==ALIVE && (m_currentMedium=nextMedium)==0 ){ setStatus(LEFT);}
    else if( m_arcLength > m_maxArcLen)                      { setStatus(DONE);}
    else {
      m_stepcount ++;
      m_maxStep = m_maxArcLen - m_arcLength + .001;
    }
  }
}

bool RCParticle::propagate(const char * stop_med_title) {

// Set max-step - done here incase a continuation of previous projection
  m_maxStep = m_maxArcLen+.0001;
  setStatus( ALIVE );

// Loop until we're out.... 
  while (status()==ALIVE)
  {
    const Medium * nextMedium=m_currentMedium;  // a temporary for new medium

        // return a Ray that estimates the trajectory (or a point)

        Ray *m_segment= 
            m_currentMedium->CreateRay( position(), momentum(), charge(), m_maxStep );

       // flag this segment
        m_segment->setFlag(m_stepcount);

        // ask the Medium for a distance along that ray, (unless
        // the particle is at rest),and for possible detector

        m_step = m_currentMedium->distanceToLeave(*m_segment, nextMedium, 100.);
        if (m_step < m_maxStep) {  // hit a boundary
            setStatus( ALIVE );
        }       
        if (m_step == FLT_MAX || (m_step<=0 && m_currentMedium==nextMedium) )   {
            WARNING( "RCParticle::propagate -- we are lost or stuck!");
            setStatus( STUCK );
            break;
        }
        else if( m_step > m_maxStep) {
            // maxStep was ignored above: assume still in same medium
            m_step = m_maxStep;
            nextMedium = m_currentMedium;
        }
        
        // increment along trajectory by arclength step
        if(m_step > 0) {
            stepBy(m_step, *m_segment );
        }

    Medium::setLastMedium(m_currentMedium);
    m_lastMedium = m_currentMedium;
    delete m_segment;

    // At a boundary... check the current Medium is the stopping place
    if(strcmp(stop_med_title, m_currentMedium->title()) == 0 ) {
        setStatus(DONE);
        return true;
    } 

    // We've reached the boundary: update current medium and see if we've left
    if ( status()==ALIVE && (m_currentMedium=nextMedium)==0 ){ setStatus(LEFT);}
    else if( m_arcLength > m_maxArcLen)                      { setStatus(DONE);}
    else {
      m_stepcount ++;
      m_maxStep = m_maxArcLen - m_arcLength + .0001; 
    }
  }
  return false;
}

void
RCParticle::stepBy(float step, Ray & path)
{
    double paverage = momentum().mag();

    setPosition(path.position(step));
    path.setArcLength(step);

    m_arcLength  += step;
    if( paverage ) {
        m_properTime += step * mass() / paverage;
        setT( t() + step * e() / paverage);
    }

    float rad_len = m_currentMedium->material().radiationLength();
    float x0;
    if(rad_len > 1.e8 || rad_len < 1.e-4) { x0 = 0.;}
    else                                  { x0 = step / rad_len;}
    m_radLength += x0; 
    m_stepList.push_back(Step(m_currentMedium, path.position(0), path.direction(), step, x0));
}

float RCParticle::mScat_Angle(float momentum, float s) const
{// method to sum up the multiple scattering contributions to the track angle 
 // over a distance s (s must be less than maxArcLen) 
  
  float scat_angle = 0.;
  float dist = 0.;
  int num_step = m_stepList.size();
  for(int i=0; i<num_step; i++){
          float x0s =    m_stepList[i].rad_len();
          float s_dist = m_stepList[i].length();
          if(dist+s_dist > s) { // pro-rate the last step
                  if(s_dist > 0) x0s *= (s - dist)/s_dist;
          }
          if(x0s > 0) {
              float ms_Angle = .014*sqrt(x0s)*(1+0.038*log(x0s))/momentum; 
              scat_angle += ms_Angle*ms_Angle;
          }
          dist += s_dist;
          if(dist >= s ) break;
  }
  return sqrt(scat_angle);       
}

float RCParticle::mScat_Angle(float momentum, int step_no) const
{// method to sum up the multiple scattering contributions to the track angle 
 // up to and including step #step_no   
  
  float scat_angle = 0.;
  int num_step = m_stepList.size();
  if(num_step > step_no) num_step = step_no;
  for(int i=0; i<num_step; i++){
          float x0s = m_stepList[i].rad_len();
          if(x0s > 0) {
              float ms_Angle = .014*sqrt(x0s)*(1+0.038*log(x0s))/momentum; 
              scat_angle += ms_Angle*ms_Angle;
          }
  }
  return sqrt(scat_angle);       
}

float RCParticle::mScat_Dist(float momentum, float s) const
{// method to sum up the multiple scattering contributions to the track 
 // displacement over a distance s (s must be less than maxArcLen) 
  
  float scat_dist = 0.;
  float dist = 0.;
  int num_step = m_stepList.size();
  for(int i=0; i<num_step; i++){
          float x0s = m_stepList[i].rad_len();
          float s_dist =m_stepList[i].length();
          float s_distp = s_dist;
          if(dist+s_dist > s) { // pro-rate the last step: s_distp
                  if(s_dist > 0) x0s *= (s - dist)/s_dist;
                  s_distp = s - dist; 
          }
          if(x0s != 0) {
              float ms_Angle = .014*sqrt(x0s)*(1+0.038*log(x0s))/momentum;
              float ms_Dist = (s - dist - s_distp)*ms_Angle;  // Disp. over remaining traj
                  float ms_sDst = s_distp*ms_Angle/1.7320508;     // Disp. within step
              scat_dist += (ms_Dist*ms_Dist) + (ms_sDst)*(ms_sDst);
          }
          dist += s_dist;
          if(dist >= s ) break;
  }
  return sqrt(scat_dist);        
}

float RCParticle::mScat_Dist(float momentum, int step_no) const
{// method to sum up the multiple scattering contributions to the track angle up to 
 // and including step #step_no  
  
  float scat_dist = 0.;
  float dist = 0.;
  int num_step = m_stepList.size();
  if(num_step > step_no) num_step = step_no;
  float total_dist = 0;
  for(int j=0; j<num_step; j++){ // find total lever arm
          total_dist += m_stepList[j].length();
  }
  for(int i=0; i<num_step; i++){
          float x0s = m_stepList[i].rad_len();
          float s_dist =m_stepList[i].length();
          if(x0s > 0) {
              float ms_Angle = .014*sqrt(x0s)*(1+0.038*log(x0s))/momentum;
              float ms_Dist = (total_dist - dist - s_dist)*ms_Angle; // Disp. over remaining traj
                  float ms_sDst = s_dist*ms_Angle/1.7320508;             // Disp within step
              scat_dist += (ms_Dist*ms_Dist) + (ms_sDst)*(ms_sDst);
          }
          dist += s_dist;
  }
  return sqrt(scat_dist);        
}

int RCParticle::findMedium(const char * name) const
{// Returns the index of the Step with the given name.
 // Returns -1 if name not found
  
  int step_no = -1;
  int num_step = m_stepList.size();
  for(int i=0; i<num_step; i++){
          step_no++;
          if(strcmp(name, m_stepList[i].place()->title()) == 0) break;
  }
  return step_no;        
}

void
RCParticle::printOn(std::ostream& str )const
{
  GParticle::printOn(str);
  str << " status: ";
  switch( m_status ) {
     case ALIVE:      str << "alive";    break;
     case LEFT:       str << "left";     break;
     case LOST:       str << "lost";     break;
     case STUCK:      str << "stuck";    break;
         case DONE:       str << "done";     break;
     default:         str << "unknown";  break;
        }
  str <<'\n';
  str <<" Total Step length: "<<m_arcLength<<"  Proper Time: "<<m_properTime<<'\n';
  str <<" Step Count: "<<m_stepcount<<'\n';

  int num_step = m_stepList.size();
  for(int i=0; i<num_step; i++){
          m_stepList[i].printOn(str);
  }
}

// Class methods for internal classStep
RCParticle::Step::Step(const Medium * med, Point x, Vector t, float s, float x0)
      :  m_region(med)
      ,  m_pos(x)
      ,  m_dir(t)
      ,  m_step_length(s)
      ,  m_rad_length(x0)
{}

void RCParticle::Step::printOn(std::ostream& str )const
{
        str<<" Region: "<<m_region->title()<<" Location:"<<m_pos<<" Direction:"<<m_dir<<'\n';
        str<<"         Step Length:"<<m_step_length<<" Radiation Length: "<<m_rad_length<<'\n';
}

---