lcsim/src/org/lcsim/recon/tracking/trfcyl
diff -u -r1.2 -r1.3
--- CylEloss.java 7 Mar 2005 06:09:48 -0000 1.2
+++ CylEloss.java 11 Sep 2007 04:57:01 -0000 1.3
@@ -4,12 +4,24 @@
import org.lcsim.recon.tracking.trfbase.Interactor;
import org.lcsim.recon.tracking.trfbase.ETrack;
+import org.lcsim.recon.tracking.trfbase.PropDir;
+import org.lcsim.recon.tracking.trfbase.Propagator;
import org.lcsim.recon.tracking.trfbase.Surface;
import org.lcsim.recon.tracking.trfbase.TrackError;
import org.lcsim.recon.tracking.trfbase.TrackVector;
import org.lcsim.recon.tracking.trfeloss.DeDx;
+import static org.lcsim.recon.tracking.trfcyl.SurfCylinder.IALF;
+import static org.lcsim.recon.tracking.trfcyl.SurfCylinder.IPHI;
+import static org.lcsim.recon.tracking.trfcyl.SurfCylinder.IQPT;
+import static org.lcsim.recon.tracking.trfcyl.SurfCylinder.ITLM;
+import static org.lcsim.recon.tracking.trfcyl.SurfCylinder.IZ;
+
+import static java.lang.Math.abs;
+import static java.lang.Math.sqrt;
+import static java.lang.Math.cos;
+
/**
* Class for modifying the covariance matrix of a track which
* has been propagated to an InteractingLayer containing
@@ -119,7 +131,86 @@
}
- //
+ public void interact_dir(ETrack theTrack, PropDir direction )
+ {
+ // This can only be used with cylinders... check that we have one..
+
+ Surface srf = theTrack.surface();
+ Assert.assertTrue( srf instanceof SurfCylinder );
+
+ // Reduced direction must be forward or backward.
+
+ Propagator.reduceDirection(direction);
+ Assert.assertTrue(direction == PropDir.FORWARD || direction == PropDir.BACKWARD);
+
+ TrackError cleanError = theTrack.error();
+ TrackError newError = new TrackError(cleanError);
+
+ TrackVector theVec = theTrack.vector();
+ TrackVector newVector = new TrackVector(theVec);
+
+ double pionMass = 0.13957; // GeV
+ double ptmax = 10000.; // GeV
+
+ double tanlm = theVec.get(ITLM);
+ double coslm = 1. / sqrt(1. + tanlm*tanlm);
+ double pinv = abs(theVec.get(IQPT)*coslm);
+
+ // make sure pinv is greater than a threshold (1/ptmax)
+ // in this case assume q = 1, otherwise q = q/pt/abs(q/pt)
+
+ int sign = 1;
+ if ( pinv < 1./ptmax )
+ pinv = 1./ptmax;
+ else
+ sign = (int) (theVec.get(IQPT)/abs(theVec.get(IQPT)));
+
+ // Evaluate the initial energy assuming the particle is a pion.
+
+ double trackEnergy = sqrt(1./pinv/pinv+pionMass*pionMass);
+
+ double trueLength = _thickness/abs(cos(theVec.get(IALF)))/coslm;
+
+ // assume the energy loss distribution to be Gaussian
+ double stdEnergy = _dedx.sigmaEnergy(trackEnergy, trueLength);
+
+ double stdMomentum = stdEnergy;
+ if(direction == PropDir.BACKWARD)
+ trueLength = -trueLength;
+ _dedx.loseEnergy(trackEnergy, trueLength);
+
+ double newMomentum = trackEnergy>pionMass ?
+ sqrt(trackEnergy*trackEnergy-
+ pionMass*pionMass): 1./pinv;
+
+ // Only vec(IQPT) changes due to E loss.
+ newVector.set(IQPT, 1./newMomentum/coslm*sign);
+
+ // Also only error(IQPT,IQPT) changes.
+ double stdVec = theVec.get(IQPT)*theVec.get(IQPT)*stdMomentum*coslm;
+ stdVec *= stdVec;
+ newError.set(IQPT, IQPT, newError.get(IQPT, IQPT) + stdVec);
+
+// TODO introduce axial definition to ETrack
+// // Axial track?
+// if(theTrack.is_axial()) {
+// newError(IPHI,IZ) = 0.;
+// newError(IZ,IZ) = 0.;
+// newError(IALF,IZ) = 0.;
+// newError(ITLM,IZ) = 0.;
+// newError(IQPT,IZ) = 0.;
+// newError(IPHI,ITLM) = 0.;
+// newError(IALF,ITLM) = 0.;
+// newError(ITLM,ITLM) = 0.;
+// newError(IQPT,ITLM) = 0.;
+// }
+
+ theTrack.setVector(newVector);
+ theTrack.setError( newError );
+
+ }
+
+//
/**
*Return the thickness of material in the cylindrical shell.
@@ -131,7 +222,7 @@
return _thickness;
}
- //
+//
/**
*Return the energy loss model used in this Interactor.
@@ -143,7 +234,7 @@
return _dedx; //cng shallow copy!
}
- //
+//
/**
*Make a clone of this object.