SwimStepOnceAction.cxx

Go to the documentation of this file.

// $Id: SwimStepOnceAction.cxx,v 1.16 2006/08/03 18:10:41 musser Exp $
//
// seun@huhepl.harvard.edu
#include "Swimmer/SwimStepOnceAction.h"
#include "Swimmer/SwimParticle.h"
#include "Swimmer/SwimStepData.h"
#include "Swimmer/SwimStepper.h"
#include "Swimmer/SwimGeo.h"
#include <cmath>

#include "MessageService/MsgService.h" // For MSG

CVSID("$Id: SwimStepOnceAction.cxx,v 1.16 2006/08/03 18:10:41 musser Exp $");

void SwimStepOnceAction::Perform(SwimParticle& particle, 
                                 SwimStepData* stepData)
{
  //====================================================================
  // Set the start point to be the previous end point and swim through 
  // again
  //====================================================================
  SwimGeo::SwimMaterial_t material = stepData->GetSwimMaterial();
  if (particle.GetMomentumModulus()!=0.0) {
    if (material==SwimGeo::kAir || 
        material==SwimGeo::kScint)
    {
      // Used to calculate path, range traveled in this step
      double density = SwimGeo::GetSwimMaterialDensity(material);
      TVector3 startpos, path;
      startpos = particle.GetPosition();
      
      // If material is Air or Scint
      // straight line approximation

      TVector3 position;
      double   xNext[3];
      double   stepSize = stepData->GetStepSize();

      if (stepData->GetIsForward()) {
        xNext[0] =
          particle.GetPosition().X()+stepSize*particle.GetDirection().X();
        xNext[1] = 
          particle.GetPosition().Y()+stepSize*particle.GetDirection().Y();
        xNext[2] = 
          particle.GetPosition().Z()+stepSize*particle.GetDirection().Z();
      }
      else {
        xNext[0] =
          particle.GetPosition().X()-stepSize*particle.GetDirection().X();
        xNext[1] = 
          particle.GetPosition().Y()-stepSize*particle.GetDirection().Y();
        xNext[2] =
          particle.GetPosition().Z()-stepSize*particle.GetDirection().Z();
      }

      //===============================================================
      // Set position, momentum of the particle
      //===============================================================
      position.SetXYZ(xNext[0],xNext[1],xNext[2]);
      particle.SetPosition(position);
  
      // Calculate path, range in this step, add to particle totals
      path = particle.GetPosition();
      double endz = path[2];
      path -= startpos;
      particle.AddS(path.Mag());
      particle.AddRange(density*path.Mag()/(Munits::g/Munits::cm2));
      MSG("Swim",Msg::kDebug)  << " start/end z " << startpos[2] <<  "/" << endz << " density " << density/1000. << " dzds " << fabs(path[2]/path.Mag()) << " dRange  " << density*path.Mag() << endl;
    }
    else {
      MSG("Swim",Msg::kDebug) << "StepOnce\n";

      bool stepOk = stepData->GetStepper()->StepOnce(particle, stepData);
      MSG("Swim",Msg::kDebug) << "StepOnce complete\n";

      if (stepOk == false) {
        // Didn't get the step we asked for
        MSG("Swim",Msg::kWarning) << "Step error\n";
      }
    }
  }
}

---