The B field is not constant in the FTPC along the z axis. Because of this, before fitting the track to
an helix the position of the clusters must be modified to take into account the real B values.
This is done in the  StFtpcTrack method "MomentumFit"
in StRoot/StFtpcTrackMaker/StFtpcTrack 

First the momentum is calculated guessing B is constant and equal to B in the position (0,0,0)
Going cluster by cluster for one track, we calculate the local momentum for each helix segment (which is
obiously bigger and bigger including new clusters). 
We read the real local measurement of the B field and we Calculate the new Momentum. Then, we modify the cluster
position taking into account the B distorsion and, with these new values, we refit the helix.
After finishing the procedure for all clusters in one track we set final momentum Values.

void StFtpcTrack::MomentumFit(StFtpcVertex *vertex)
{
  Double_t xWeight[11], yWeight[11];
  Double_t xval[11], yval[11], zval[11];
  Double_t xhelix[11], yhelix[11], zhelix[11];
  Int_t i, j;
  
  mIterSteps = 10;
  mYCenter = 0.;
  mXCenter = 0.;

  if (vertex) {
    //additional point needed
    mVertexPointOffset = 1;
    
    //vertex used as additional point on track
    xval[0] = vertex->GetX() * centimeter;
    yval[0] = vertex->GetY() * centimeter;
    zval[0] = vertex->GetZ() * centimeter;
    
    // use vertex error as weight (if it exists)
    xWeight[0] = (vertex->GetXerr() == 0.) ? 100. : 1./(vertex->GetXerr() * centimeter);
    yWeight[0] = (vertex->GetYerr() == 0.) ? 100. : 1./(vertex->GetYerr() * centimeter);
  }

  else {
    //no additional point needed
    mVertexPointOffset = 0;
  }
  
  Int_t backw_counter = 0;
  // initialize position arrays
  // these values will later be manipulated, mPoint array will not be touched   
  for(i = 0; i < GetNumberOfPoints(); i++) {
    backw_counter = GetNumberOfPoints() - 1 - i + mVertexPointOffset;
    StFtpcPoint *point = (StFtpcPoint*)mPoints->At(i);

    xval[backw_counter] = point->GetX() * centimeter;
    yval[backw_counter] = point->GetY() * centimeter;
    zval[backw_counter] = point->GetZ() * centimeter;
    
    // hit resolution taken from hit error (or set to 100)
    xWeight[backw_counter] = (point->GetXerr() == 0.) ? 100. : 1./(point->GetXerr() * centimeter);
    yWeight[backw_counter] = (point->GetYerr() == 0.) ? 100. : 1./(point->GetYerr() * centimeter);
  }

  /////////////////////////////////////////////////////////////////////
  // calculate first guess momentum from helix fit
  /////////////////////////////////////////////////////////////////////
  
  CircleFit(xval, yval, xWeight, yWeight, GetNumberOfPoints() + mVertexPointOffset);
  LineFit(xval, yval, zval, xWeight, yWeight, GetNumberOfPoints() + mVertexPointOffset);

  // determine helix parameters
  Double_t dipAngle = fabs(atan(1/(mFitRadius*mArcSlope)));
  
  if (zval[1] < 0) {
    dipAngle *= -1;
  }
  
  Double_t startPhase = atan((yval[0]-mYCenter)/(xval[0]-mXCenter));
  
  if (xval[0]-mXCenter < 0) {
    startPhase += pi;
  }
  
  else if (yval[0]-mYCenter < 0) {
    startPhase += twopi;
  }
  
  Int_t orientation = 1;

  if (mArcSlope * zval[1] < 0) {
    orientation = -1;
  }

  // create helix
  Double_t startAngle = mArcOffset + mArcSlope * zval[0];
  Double_t startX = mXCenter + mFitRadius * cos(startAngle);
  Double_t startY = mYCenter + mFitRadius * sin(startAngle);

  StThreeVector<Double_t> startHit(startX, startY, zval[0]);
  setParameters(1/mFitRadius, dipAngle, startPhase, startHit, orientation);

  // get z-component of B-field at 0,0,0 for first momentum guess
  Float_t pos[3] = {0, 0, 0};
  Float_t centralField[3];
  StFtpcTrackingParams::Instance()->MagField()->B3DField(pos, centralField);
  centralField[0] *= kilogauss;
  centralField[1] *= kilogauss;
  centralField[2] *= kilogauss;
  mZField = (Double_t) centralField[2];
  
  // get momentum at track origin and charge
  StThreeVector<Double_t> rv(0, 0, zval[0]);
  StThreeVector<Double_t> nv(0, 0, 1);
  Double_t pl = pathLength(rv, nv);
  mHelixMomentum = momentumAt(pl, mZField);
  SetCharge(charge(mZField));

  // store helix fitted hit positions
  for(i = 0; i < GetNumberOfPoints() + mVertexPointOffset; i++) {
    StThreeVector<Double_t> rvec(0, 0, zval[i]);
    StThreeVector<Double_t> nvec(0, 0, 1);
    Double_t plength = pathLength(rvec, nvec);
    xhelix[i] = x(plength);
    yhelix[i] = y(plength);
    zhelix[i] = z(plength);
  }

  ///////////////////////////////////////////////////////////////////////
  // track helix momentum through measured field:
  ///////////////////////////////////////////////////////////////////////

  // initialize position and momentum
  StThreeVector<Double_t> currentPosition(xhelix[0 + mVertexPointOffset], 
					  yhelix[0 + mVertexPointOffset],
					  zhelix[0 + mVertexPointOffset]);
  pl = pathLength(currentPosition, nv);
  StThreeVector<Double_t> currentMomentum(momentumAt(pl, mZField));


  if (StFtpcTrackingParams::Instance()->MagFieldFactor()) {

    // iterate over points
    Double_t stepSize;
    
    for(i = 1 + mVertexPointOffset; i < GetNumberOfPoints() + mVertexPointOffset; i++) {
      stepSize = (zval[i] - zval[i-1]) / mIterSteps;
      
      // iterate between points
      for(j = 0; j < mIterSteps; j++) {
	// store momentum for position propagation
	Double_t propagateXMomentum = currentMomentum.x();
	Double_t propagateYMomentum = currentMomentum.y();
	Double_t propagateZMomentum = currentMomentum.z();
	
	// get local magnetic field
	Float_t positionArray[3] = {currentPosition.x(), 
				    currentPosition.y(), 
				    currentPosition.z() + stepSize/2};
	Float_t localField[3];
	StFtpcTrackingParams::Instance()->MagField()->B3DField(positionArray, localField);
	
	StThreeVector<Double_t> fieldVector
	  ((Double_t) localField[0] * kilogauss/tesla*c_light*nanosecond/meter, 
	   (Double_t) localField[1] * kilogauss/tesla*c_light*nanosecond/meter, 
	   (Double_t) localField[2] * kilogauss/tesla*c_light*nanosecond/meter); 
	
	// calculate new momentum as helix segment
	Double_t absMomentum = abs(currentMomentum);
	StThreeVector<Double_t> perpField = 
	  currentMomentum.cross(fieldVector) * (Double_t)mQ/absMomentum;
	Double_t twistRadius = (absMomentum/abs(perpField)) * meter/GeV;
	
	Double_t stepLength = stepSize/cos(currentMomentum.theta());
	
	Double_t newMomentumCross = absMomentum*stepLength/twistRadius;
	Double_t newMomentumParallel = 
	  sqrt(absMomentum * absMomentum - newMomentumCross * newMomentumCross);
	currentMomentum.setMagnitude(newMomentumParallel);
	StThreeVector<Double_t> momentumChange(perpField);
	momentumChange.setMagnitude(newMomentumCross);
	currentMomentum = currentMomentum+momentumChange;
	
	// propagate position
	propagateXMomentum = (propagateXMomentum + currentMomentum.x())/2;
	propagateYMomentum = (propagateYMomentum + currentMomentum.y())/2;
	propagateZMomentum = (propagateZMomentum + currentMomentum.z())/2;
	currentPosition.setX(currentPosition.x() + stepSize *
			     (propagateXMomentum / propagateZMomentum));
	currentPosition.setY(currentPosition.y() + stepSize *
			     (propagateYMomentum / propagateZMomentum));
	currentPosition.setZ(currentPosition.z() + stepSize);
      }
      
      // change position array to compensate for distortion
      StThreeVector<Double_t> rvec(0, 0, zval[i]);
      StThreeVector<Double_t> nvec(0, 0, 1);
      Double_t plength=pathLength(rvec, nvec);
      
      if (zval[1] > 0) {
	xval[i] += (x(plength) - currentPosition.x());
	yval[i] += (y(plength) - currentPosition.y());
      }
      
      else {
	xval[i] -= (x(plength) - currentPosition.x());
	yval[i] -= (y(plength) - currentPosition.y());
      }
      
      // calculate fit quality indicators only if needed
      //       Double_t distHitHelix=sqrt((xval[i]-x(plength))*(xval[i]-x(plength))+(yval[i]-y(plength))*(yval[i]-y(plength)));
      //       Double_t distHelixFit=sqrt((x(plength)-currentPosition.x())*(x(plength)-currentPosition.x())+(y(plength)-currentPosition.y())*(y(plength)-currentPosition.y()));
      
    }
    
    //////////////////////////////////////////////////////////////////////
    // refit helix
    //////////////////////////////////////////////////////////////////////
    
    CircleFit(xval, yval, xWeight, yWeight, GetNumberOfPoints()+mVertexPointOffset);
    LineFit(xval, yval, zval, xWeight, yWeight, GetNumberOfPoints()+mVertexPointOffset);
    
    // determine helix parameters
    dipAngle = fabs(atan(1/(mFitRadius*mArcSlope)));
    
    if (zval[1] < 0) {
      dipAngle*=-1;
    }
    
    startPhase = atan((yval[0]-mYCenter)/(xval[0]-mXCenter));
    
    if (xval[0] - mXCenter < 0) {
      startPhase+=pi;
    }
    
    else if (yval[0]-mYCenter<0) {
      startPhase+=twopi;
    }
    
    orientation = 1;
    
    if (mArcSlope * zval[1] < 0) {
      orientation = -1;
    }
    
    // set helix parameters to new values
    startAngle = mArcOffset + mArcSlope * zval[0];
    startX = mXCenter + mFitRadius * cos(startAngle);
    startY = mYCenter + mFitRadius * sin(startAngle);
    startHit.setX(startX);
    startHit.setY(startY);
    
    setParameters(1/mFitRadius, dipAngle, startPhase, startHit, orientation);
  }
  
  // set final momentum value
  pl = pathLength(rv, nv);
  mFullMomentum = momentumAt(pl, mZField);
}