SwimDefStepper.cxx

Go to the documentation of this file.

// SwimDefStepper
//
#include "Swimmer/SwimDefStepper.h"
#include "Swimmer/SwimParticle.h"
#include "Swimmer/SwimStepData.h"
#include "BField/BField.h"
#include "Conventions/Munits.h"
#include "Conventions/Mphysical.h"
#include <iostream>
#include "TMath.h"
#include "MessageService/MsgService.h"

CVSID("$Id: SwimDefStepper.cxx,v 1.20 2008/02/20 11:40:34 bckhouse Exp $");

//......................................................................

SwimDefStepper::SwimDefStepper(BField* f) : fBfield(f) 
{
  this->SetUpIntegrator();
}

//......................................................................

SwimDefStepper::SwimDefStepper(BField* f, double min, double acc) : 
  fMinStepSize(min),
  fAcc(acc),
  fBfield(f)
{
  this->SetUpIntegrator();
}

//......................................................................

void SwimDefStepper::SetUpIntegrator() 
{
//======================================================================
// Purpose: Setup the integrator object
//======================================================================
  double s[6] = {1.0,1.0,1.0,  // Typical sizes of x,y,z and
                 1.0,1.0,1.0}; // px,py,pz
  // Setup the integration object
  fIntegrator.SetNvar(6);
  fIntegrator.SetScales(s);
  fIntegrator.SetAccuracy(fAcc);
  fIntegrator.SetMinStepSize(fMinStepSize);
  fIntegrator.SetDerivFunc(fDerivFunc);
}

//......................................................................

bool SwimDefStepper::StepOnce(SwimParticle& particle, 
                              SwimStepData* stepData)
{
  // Setup step
  fX[0] = particle.GetPosition().X();
  fX[1] = particle.GetPosition().Y();
  fX[2] = particle.GetPosition().Z();
  fX[3] = particle.GetMomentum().X();
  fX[4] = particle.GetMomentum().Y();
  fX[5] = particle.GetMomentum().Z();
  double stepSize = stepData->GetStepSize();
  // Used to calculate path, range traveled in this step
  double density = SwimGeo::GetSwimMaterialDensity(stepData->GetSwimMaterial());
  TVector3 startpos, path, startmom, deltamom;
  startpos = particle.GetPosition();
  startmom = particle.GetMomentum();

  fDerivFunc.SetChargeBField(particle.GetCharge(),fBfield);
  // TVector3 Bxyz = fBfield->GetBField(startpos);
  // cout << " pos " << fX[0] << " " << fX[1] << " z  " << fX[2] << " " << fX[4]/fX[5] << " "  << fX[5] << endl ;
  // cout << " density " << density << endl;
  //  cout << " bfield " << Bxyz.X() << " " <<      Bxyz.Y() << " " <<  Bxyz.Z() << " " <<   endl;
  // Integrate across the step
  fIntegrator.SetInitialCondition(fX);
  fIntegrator.SetStartPoint(0.0);
  if (stepData->GetIsForward())
    fIntegrator.SetEndPoint(stepSize);
  else
    fIntegrator.SetEndPoint(-stepSize);

  // 20 Feb 2008 - Avoid assert(stepSize>0) in fIntegrator.
  if(stepSize > 0){
    fIntegrator.SetStepSize(stepSize);
  }
  else{
    MSG("Swim", Msg::kError) << "Unexpected stepSize of " << stepSize << " bailing...\n";
    return false;
  }

  bool aok = fIntegrator.Integrate(fX);

  // Shuffle to output
  TVector3 position;
  TVector3 momentum;
  double   pModulus2 = fX[3]*fX[3]+fX[4]*fX[4]+fX[5]*fX[5];
  position.SetXYZ(fX[0],fX[1],fX[2]);
  particle.SetPosition(position);
  if (pModulus2!=0.0)
    momentum.SetXYZ(fX[3],fX[4],fX[5]);
  else
    momentum.SetXYZ(0.0,0.0,0.0);
  particle.SetMomentum(momentum);

  // Calculate path, range in this step, add to particle totals
  path = particle.GetPosition();
  
  double endz =path[2];
  deltamom = particle.GetMomentum();
  path -= startpos;
  deltamom -= startmom;
  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;

  particle.AddVxB(deltamom.Mag());
  
  return aok;
}

//......................................................................

SwimDefStepper::DerivFunc::DerivFunc(double q, BField* b) 
{
  // Initialize force coeff. and BField
  this->SetChargeBField(q,b);
}

//......................................................................

void SwimDefStepper::DerivFunc::SetChargeBField(double q, BField* b)
{
  // The units conversion is c[m/s]*[eV]/[GeV]
  static const double unitsC = 
    (Mphysical::c_light/(Munits::m/Munits::s))/1.0e9;
  fForceCoef = unitsC*q;
  fBfield    = b;
}

//......................................................................

void SwimDefStepper::DerivFunc::operator()(double /* x */, 
                                           const double* y, 
                                           double* dydx) 
{
  //const double oneOverSqrt2 = 1./sqrt(2.);
  // y[0] - x coordinate (base length units)
  // y[1] - y coordinate (base length units)
  // y[2] - z coordinate (base length units)
  // y[3] - x momentum   (base momentum units)
  // y[4] - y momentum   (base momentum units)
  // y[5] - z momentum   (base momentum units)
  double ptotSqr = y[3]*y[3]+y[4]*y[4]+y[5]*y[5];
  if (ptotSqr != 0.0) {
    double oneOverPtot = 1.0/TMath::Sqrt(ptotSqr);      
    dydx[0] = y[3]*oneOverPtot; // (d/ds)x = Px/Ptot
    dydx[1] = y[4]*oneOverPtot; // (d/ds)y = Py/Ptot
    dydx[2] = y[5]*oneOverPtot; // (d/ds)z = Pz/Ptot
    
    //the point is passed in in xyz, even though the measurements are in uvz - 
    //the rotation is done in the fitter.
    TVector3 xyz(y[0],y[1],y[2]);

    //since the Bfield is also in xyz, we are all happy.
    TVector3 B = fBfield->GetBField(xyz);

    // Calculate derivatives
    double qOverE = fForceCoef*oneOverPtot;
    dydx[3] = (y[4]*B[2]-y[5]*B[1])*qOverE; // dP/ds = q*(p x B)/E/v
    dydx[4] = (y[5]*B[0]-y[3]*B[2])*qOverE; // dP/ds = q*(p x B)/E/v
    dydx[5] = (y[3]*B[1]-y[4]*B[0])*qOverE; // dP/ds = q*(p x B)/E/v
  }
  else {
    // Particle is stopped
    dydx[0] = 0.0;
    dydx[1] = 0.0;
    dydx[2] = 0.0;
    dydx[3] = 0.0;
    dydx[4] = 0.0;
    dydx[5] = 0.0;
  }
}

---