---

MCParticle.cxx

Go to the documentation of this file.

//  $Id: MCParticle.cxx,v 1.8 2001/09/21 18:38:54 atwood Exp $
//
// This file is part of Gismo 2
//
//
// Implemention of all MCParticle member functions

#include "gismo/MCParticle.h"

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

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

#include "facilities/error.h"
#include "geomrep/TrackRep.h"
#include "Random.h"


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

MCParticle::MCParticle(const PData* data, float mass)
        : GParticle(data,mass)
        , _track(0)
        , _status(ALIVE)
        {}
MCParticle::MCParticle(long id, float mass)
        : GParticle(id,mass)
        , _track(0)
        , _status(ALIVE)
        {}
MCParticle::MCParticle(const char* name, float mass)
        : GParticle(name,mass)
        , _track(0)
        , _status(ALIVE)
        {}
MCParticle::MCParticle(const MCParticle& p)
        : GParticle(p)
        , _track(0)
        , _status(ALIVE)
        {}

MCParticle::~MCParticle()
{
        if(_track)
                delete _track;
}
GParticle*
MCParticle::addChild(const PData* childData, float mass)
{
    MCParticle* child = new MCParticle(childData, mass);
    GParticle::addChild(child);
    return child;
}
//
MCParticle*
MCParticle::child(int n)const
{
    return (MCParticle*)(GParticle::child(n));
}
const MCParticle*
MCParticle::parent()const
{
    return (const MCParticle*)(GParticle::parent());
}


const float MCParticle::MAX_TIME=1000;      //cm...
const float MCParticle::MAX_BEND=31;       //radians .. 5 loops...
const int MCParticle::MAX_GENERATION=500;
const int MCParticle::MAX_STEPS=10000;

//  MCParticle::propagate -- propagate the particle,and all its progeny,
//   out of the medium
//
void MCParticle::propagate(const Medium * medium) {

  currentMedium = medium;
  arcLength = properTime = turnAngle = 0.;
  timeOfDecay = GenerateTimeOfDecay();
  startingEnergy = e();
  if( ! parent() ) generation =0;
  stepcount = 0;
  Medium::setLastMedium(0);
  lastMedium = medium;
  Vector starting_momentum = momentum(); // save to restore when we are done

  while (status()==ALIVE)
  {
    const Medium * nextMedium=currentMedium;  // a temporary for new medium

    //  find maximum step, reflecting decay or interaction
    float maxStep =  checkStep( );

    if( maxStep==0)    {
        step = 0;  // means that it stopped: skip the next stuff
        s_segment = 0;
    }
    else {
        // return a Ray that estimates the trajectory (or a point)

        s_segment=
            currentMedium->CreateRay( position(), momentum(), charge(), maxStep );

        // flag the segment for use by geometry
        s_segment->setFlag(stepcount); 

        // ask the Medium for a distance along that ray, (unless
        // the particle is at rest),and for possible detector
        step = currentMedium->distanceToLeave(*s_segment, nextMedium, maxStep);

        if (step < maxStep) {  // hit a boundary: cancel decay or interaction
            setStatus( ALIVE );
        }       
        if (step == FLT_MAX || (step<=0 && currentMedium==nextMedium) ) {
            WARNING( "MCParticle::propagate -- we are lost or stuck!");
            setStatus( STUCK );
            break;
        }
        else if( step > maxStep) {
            // maxStep was ignored above: assume still in same medium
            step = maxStep;
            nextMedium = currentMedium;
        }
        
        // perform the  step, taking all physics into account, and
        if(step > 0) {
            stepBy(step, *s_segment );
        }

    }
    // peform decays or interactions at end of step
    endStep();


#ifdef DUMP_STEP
    fout << z() << '\t' << step << '\t' << currentMedium->title() << '\t'
        << currentMedium->material().name()<< '\t'
        << energyLoss()*1e6 << '\n';
#endif

    // notify the associated detector, if any
    Detector* d =currentMedium->detector();
    if( d )
      d->score(this);

    // add the segment for this step to the track, or delete it
    if(s_segment) {
       if( saveTrack && step > 0) {
          if( !_track )
             _track = new Track(s_segment, charge()!=0);
          else
             _track->addSegment(s_segment);
        }
       else
          delete s_segment;
    }

    startingEnergy = e();
    Medium::setLastMedium(currentMedium);
    lastMedium = currentMedium;

    // if particle didn't stop, decay, or interact, it must have
    // reached the boundary: update current medium and see if left

#if 0  // special debug code to enable by hand to study geometry problems
    if( status()==ALIVE && nextMedium != 0 ) {
        const Medium* temp;
        float d = nextMedium->distanceToLeave( Ray(position(),direction()),temp, 1000.);
        if( d>= 1000.)
            cerr << "MCParticle::propagate--misstep from "<< currentMedium->title()
                << " to " << nextMedium->title() << '\n';
        
    }
#endif

    if ( status()==ALIVE && (currentMedium=nextMedium)==0 )
        setStatus(LEFT);
    else
      stepcount ++;

  }
  // Finished stepping this particle. Restore initial momentum (and energy) for user access
  setMomentum(starting_momentum);

  // Note that a "Stopping" particle may have children! This can be used by
  // score() to ignore children and stop propagating
  if(status() == STOPPED){
      return;
  }


  // Now  recurse: make sure that all of the children
  // of a particle's  decays or interactions propagate out of
  // the medium in which they were created
  const Medium* parentMedium = currentMedium;
  generation ++;
  for( int i=0;i<numChildren(); i++ ) {
     child(i)->propagate(parentMedium);
  }
  generation --;

  // Finally, if not displaying, get rid of children here
  if( generation > saveTrack) removeChildren();
}

float MCParticle::checkStep()

//  calculate maximum step, allowing for decay, or interaction

{
    float pmagnitude =  momentum().mag();
    float distance = FLT_MAX;

    // get minimum energy allowed by medium for this part.
    float eCutMin = currentMedium->kECutOff() + mass();
    ecut = currentMedium->material().isVacuum()? 0: interactor()->ecut(this,currentMedium);
    if( ecut>0 &&  ecut < eCutMin) ecut = eCutMin;
        
    if( ( e() < ecut )
         || generation >= MAX_GENERATION
        || stepcount > MAX_STEPS)
    {
        if( generation >= MAX_GENERATION )
            WARNING("MCParticle::propagate - stopped particle because too many generations");
        if( stepcount > MAX_STEPS )
            WARNING("MCParticle::propagate - too many steps, stopping");
        

        setStatus( STOPPED );
        if( charge() && pmagnitude > 0. )
            setMomentum(Vector(0.,0.,0.));

        distance = 0.;
    }

    // see if decayed before taking step

    if (mass()>0. && timeOfDecay < MAX_TIME )  {
        float distToDecay =
                 (timeOfDecay-properTime)*eta();
        if (distToDecay <= distance)  // will catch decay at rest
        {   setStatus( DECAYED);
            distance = (distToDecay > 0.) ? distToDecay : 0.;
        }
    }

    // if in material, check range and distance to interact, compare w/ decay

    if (!currentMedium->material().isVacuum())    {

        if (charge() && distance > 0. ) {

            // compare with maximum from Material or Medium
            float maxstep = interactor()->maxStepSize(this,currentMedium);
            if( maxstep < distance ) {
                distance = maxstep;
                setStatus(ALIVE);  // cutting back: cancel decay if any
            }
         }

         double intlen =interactor()->interactionLength(this,currentMedium),
                distToInteract = distance;
         if( intlen !=FLT_MAX && intlen > 0.
                && intlen < 1E5    ) // avoid problem with exponential
              distToInteract = Random::exponential( intlen );

         if (distToInteract < distance) {
             setStatus( INTERACTED);
             distance = distToInteract;
         }
    }
    else if (distance != 0 && status()!= DECAYED )
        distance = currentMedium->maxStep();
    else if ( status()==DECAYED && distance >currentMedium->maxStep() ){
        setStatus(ALIVE); // cancel decay if too far
        distance = currentMedium->maxStep();
    }
    return  distance;
}

void
MCParticle::stepBy(float step, Ray & path)
{

    double paverage = momentum().mag();
    if (paverage==0 || step==0 ) return; // shouldn't happen

    if( charge() ) {
        Vector dir = path.direction(step);
        setMomentum( paverage * dir);
        //  apply energy loss, if appropriate

        interactor()->afterStep(step, this,currentMedium);
        float pfinal   = momentum().mag();
        if (pfinal==0.)
        {
            // special case: at rest, so decay it or declare it stopped
            if (timeOfDecay > 0 && timeOfDecay < MAX_TIME)
            {
                setT( t() + timeOfDecay );
                setStatus( DECAYED );
            }else
                setStatus( STOPPED );
        }
        paverage = ( paverage + pfinal)/2.;
        turnAngle  += step*path.curvature();
        if ( fabs(turnAngle) > MAX_BEND )
                        setStatus( LOOPING );
    }

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

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

}

void
MCParticle::endStep()
{
    if( status() ==STOPPED )
    {

        // if at rest, it might decay or annihilate

        if (timeOfDecay > 0 && timeOfDecay < MAX_TIME)
        {  setT(t() + timeOfDecay);   // is this right????
                        setStatus( DECAYED );
        }
        else if( charge()
                         && interactor()->interactionLength(this,currentMedium)==0 )
        {
                // this case for a stopped positron or anti-proton
                setStatus( INTERACTED );
        }

     }
     // Now actually perform interaction or decay

    if  (status()==DECAYED)
                decayOnce();

    else if  (status()==INTERACTED) {
         // status has been set to interacted, but if only E-loss,turn if off
        if( !interactor()->interact(this,currentMedium))
                setStatus( ALIVE );
    }
    return;

}

void
MCParticle::printOn(std::ostream& str )const
{
  GParticle::printOn(str);
  str << " status: ";
  switch( _status ) {
     case ALIVE:      str << "alive";    break;
     case STOPPED:    str << "stopped";  break;
     case DECAYED:    str << "decayed";  break;
     case INTERACTED: str << "interacted";break;
     case LOOPING:    str << "looping";  break;
     case LEFT:       str << "left";     break;
     case STUCK:      str << "stuck";    break;
     case PRIMARY:    str << "primary";  break;
     case SHOWER:     str << "shower";  break;
     default:         str << "unknown";  break;
        }
}

void MCParticle::addInteractor(Interactor * i)
{
   thePDT()->addInteractor(i);
}


// allocate static variables

Ray * MCParticle::s_segment;               // Trajectory of last step

const Medium* MCParticle::currentMedium; // the current Medium for propagation
float MCParticle::arcLength;             // total distance swum so far
float MCParticle::properTime;            // total proper time since creation.
float MCParticle::turnAngle;             // total angular deflection so far
float MCParticle::timeOfDecay;           // actual proper time of decay
double MCParticle::startingEnergy;        // keep track of last reported energy
                                         // for eloss calculation
float MCParticle::step;                  // current step length
float MCParticle::ecut;                  // minimum energy allowed

int MCParticle::stepcount;               // step counter
int MCParticle::generation;              // generation counter

int MCParticle::saveTrack=0;

---