lcsim/src/org/lcsim/util/step
diff -u -r1.2 -r1.3
--- CamKalmanStepper2.java 30 Nov 2007 17:05:00 -0000 1.2
+++ CamKalmanStepper2.java 14 Nov 2008 23:35:27 -0000 1.3
@@ -16,482 +16,824 @@
import org.lcsim.geometry.Layered;
import org.lcsim.geometry.LayeredSubdetector;
import org.lcsim.material.*;
- /**
- * CamKalmanStepper2- A Kalman Stepper
- * trasnsported from JAS2
- *JINST - published October 2006
- *arXiv: physics/0604197
- *@Author-C. Milstene <[log in to unmask]>
- * Code developped with Fedor Ignatov
- * and Hans Wenzel consulting
- * More comments and references to come
- * Hard-coded material Iron for the dE/dx
- * (Testing program-TestKalm.java contains
- * a main- and is the (Driver/MainProject))
- */
+/**
+ * CamKalmanStepper2- A Kalman Stepper
+ * trasnsported from JAS2- A fukll version in
+ * JINST - published October 2006 following address
+ * http://www.iop.org/EJ/abstract/1748-0221/1/10/P10003
+ * Partial versions in 2 arXiv:physics
+ * arXiv: physics/0604197; physics/0605015
+ *@Author Caroline Milstene <[log in to unmask]>
+ * Code developped with Fedor Ignatov'shelp
+ * and modified consulting Hans Wenzel and
+ * Rob Kutschke
+ * (Testing program-TestCamKalm.java contains
+ * a main- and is the (Driver/MainProject))
+ */
public class CamKalmanStepper2{
- private boolean Debug=false;
- private final double kAlmost1=0.999;
- private final double kAlmost0=1E-33;
- private final double kdistsmall=0.001;
- private final double kVeryBig =(1./kAlmost0);
- private final double kB2C=0.299979245*1.E-03;
- private final double kAlmost0Field =1.E-33;
- private final double kMostProbablePt=0.35;
- private final double cmTomm=10.;
- private final double mevTogev=1./1000.;
- private final double gevTomev=1000.;
- private double [] zRhoZoverA;
- private double kQ;
- private double kB=5.; // default- reset in propagateTo
- private double kMass;
- private String MaterialName;
- private boolean kRungeKutta;
- public boolean kupDate = false;
- private static boolean ktransport=true;
- private static boolean stopTkELow=false;
- private double [] kRp = new double [6];
- private double [] kRPlus ;
- private double [] measState;
- private double [][]measCov=new double[3][3];
- private double [] nwRp = new double[8];
- private Matrix PkMinus = new Matrix(6,6);
- private Matrix PkPlus = new Matrix(6,6) ;
- private Matrix Hk ;
- private double[][] kCovariance = new double [6][6];// construct a 6 by 6 array of zeroes
- private Matrix fi = new Matrix(6,6); // construct
- private AIDA aida = AIDA.defaultInstance();
- private CalHitMapMgr evtptr = null;
- private boolean StopTkElow=false;
- private Random rndcm = new Random();
- private Detector det ;
- private static int ncount1=0;
- /** default constructor*/
- public CamKalmanStepper2(){
- evtptr = CalHitMapMgr.getInstance();
- kQ = -1. ;
- kMass = 0.105;
- kRungeKutta = true;
- for(int i=0; i<6; i++)kRp[i]=0.;
- for(int i=0; i<6; i++){
- for(int j=0; j<6; j++)kCovariance[i][j]=0.; // or any starting value
- }
- }
- /** Copy Constructor-
- *@param rp an Array with the particle position and Momnetum
- * (x,y,z,px,py,pz) - a phase-space point
- */
- public CamKalmanStepper2(double [] rp){
- evtptr = CalHitMapMgr.getInstance();
- kQ = rp[6] ;
- kMass = rp[7] ;
- kRungeKutta = true;
- for(int i=0; i<6; i++)this.kRp[i]=rp[i];
- for(int i=0; i<6; i++){
- for(int j=0; j<6; j++)this.kCovariance[i][j]=kCovariance[i][j];
- }
- }
- /** Copy Constructor-
- *@param rp an array with the particle position and momentum(x,y,z,px,py,pz)
- *@param cova a covariant matrix
+ private static boolean stopTkELow = false;
+ private static boolean Debug = false;
+ private static boolean ktransport = true;
+ private static boolean materialByName = false;
+ private final double kAlmost1 = 0.999;
+ private final double kAlmost0 = 1E-33;
+ private final double kdistsmall = 0.000001;
+ private final double kVeryBig = (1./kAlmost0);
+ private final double kB2C = 2.99792458*1.E-04;
+ private final double kAlmost0Field =1.E-33;
+ private final double kMostProbablePt=0.35;
+ private final double cmTomm = 10.;
+ private final double mevTogev = 1./1000.;
+ private final double gevTomev = 1000.;
+ private double [] zRhoZoverA;
+ private double kQ;
+ private double kB = 5.; // default- reset in propagateTo
+ private double kMass ;
+ private double _dedx = 1.E-20, _radLength=1.E-20, _density=1.E-19;
+ private static String _materialName;
+ public boolean kupDate = false;
+ private double [] kRp = new double [6];
+ private double [][]measCov = new double[3][3];
+ private double [] nwRp = new double[8];
+ private double _pNow,_ptNow, _rNow;
+ private Matrix M6x6tmp = new Matrix(6,6);
+ private Matrix M3x6tmp = new Matrix(3,6);
+ private Matrix M3x3tmp = new Matrix(3,6);
+ private Matrix II = Matrix.identity(6,6);
+ private Matrix Zero6 = new Matrix(6,6);
+ private Matrix PkMinus = new Matrix(6,6);
+ private Matrix PkPlus = new Matrix(6,6);
+ // measured position x,y,z
+ double [][]val1 = {{0.,0.,0.,1.,0.,0.},{0.,0.,0.,0.,1.,0.},{0.,0.,0.,0.,0.,1.}};
+ private Matrix Hk =new Matrix(val1) ;
+ private double[][] kCovariance = new double [6][6];// construct a 6 by 6 array of zeroes
+ private Matrix fi ;
+ private CalHitMapMgr evtptr = null;
+ private Random rndcm = new Random();
+ private Detector det ;
+ public static AIDA aida = AIDA.defaultInstance();
+ private static int ncount1=0;
+ /** default constructor*/
+ public CamKalmanStepper2(){
+ evtptr = CalHitMapMgr.getInstance();
+ kQ = 1. ;
+ kMass = 0.107;
+ for(int i=0; i<6; i++)kRp[i]=0.;
+ setNewRp();
+ for(int i=0; i<6; i++){
+ for(int j=0; j<6; j++)kCovariance[i][j]=0.; // or any starting value
+ }
+ }
+ /** Copy Constructor-
+ *@param rp an Array with the particle position and Momnetum
+ * (x,y,z,px,py,pz) - a phase-space point
+ * kQ=rp[6] , kMass=rp[7] the particle Charge and Mass
+ */
+ public CamKalmanStepper2(double [] rp){
+ evtptr = CalHitMapMgr.getInstance();
+ kQ = rp[6] ;
+ kMass = rp[7] ;
+ for(int i=0; i<6; i++){double xyptmp=rp[i];this.kRp[i]=xyptmp;}
+ setNewRp();
+ for(int i=0; i<6; i++){
+ for(int j=0; j<6; j++){
+ double covtmp=kCovariance[i][j];
+ this.kCovariance[i][j]=covtmp;
+ }
+ }
+ }
+ /** Copy Constructor-
+ *@param rp an array with the particle position and momentum(x,y,z,px,py,pz)
+ *kQ and kMass the particle charge and Mass in rp[6] and rp[7]
+ *@param covar a covariant matrix
+ */
+ public CamKalmanStepper2(double [] parm, double [][]covar){ // Copy constructor
+ evtptr = CalHitMapMgr.getInstance();
+ kQ = parm[6] ;
+ kMass = parm[7] ;
+ setRpCov(parm,covar); // create external track parameters from given arguments
+ }
+ /** set rp and covariant Matrix
+ *@param rp an array with the particle position and momentum
+ *the particle charge and Mass in rp[6] and rp[7]
+ *@param covar the covariant matrix set to covar
+ */
+ public void setRpCov(double [] rp, double [] [] covar){
+ for(int i=0; i<6; i++){double xyptmp=rp[i]; this.kRp[i]=xyptmp;}
+ kQ=rp[6]; kMass=rp[7];
+ setNewRp();
+ for(int i=0; i<6; i++){
+ for(int j=0;j<6;j++){
+ double cov=covar[i][j]; kCovariance[i][j] = cov;
+ }
+ }
+ double zero=0., scal=covar[1][1];
+ PkPlus= PkPlus.times(zero).plus(II.times(scal));
+ }
+ /** Propagate the particle from a position r
+ * to s new position in a magnetic field b
+ *@param r the position
+ *@param b the magnetic field
+ *@param matInfo={dedx,radLength/density)
+ * energy loss by ionization in material & number radiation Length.
+ */
+ public int [] propagateInBarrelByDr(double rToGo, double b,double[] matInfo){
+ int[]numSteps=new int[2];
+ double dr = rToGo-Math.sqrt(kRp[0]*kRp[0]+kRp[1]*kRp[1]);
+ double pm = Math.sqrt(kRp[3]*kRp[3]+kRp[4]*kRp[4]+kRp[5]*kRp[5]);
+ double bStep = (pm>1.)?5.:2.; // MagField step 5mm above pm=1GeV and 2mm below
+ int nStepsToDo = (int)(Math.abs(dr/bStep)+0.5), nStepsDone=nStepsToDo;
+ numSteps = propagateInBarrelByDr(rToGo,b,nStepsToDo,matInfo);
+ return (numSteps);
+ }
+ /**
*/
- public CamKalmanStepper2(double [] parm, double [][]covar){ // Copy constructor
- evtptr = CalHitMapMgr.getInstance();
- kQ = parm[6] ;
- kMass = parm[7] ;
- kRungeKutta = true;
- set(parm,covar); // create external track parameters from given arguments
- if(Debug)
- System.out.println("First value of kRp[0]"+ kRp[0]+" kRp[1]"+kRp[1]+" kRp[2]"+ kRp[2]);
- }
-/** set rp and covariant Matrix
- *@param rp an array with the particle position and momentum
- *@param covar the covariant matrix
- */
- public void set(double [] rp, double [] [] covar){
- for(int i=0; i<6; i++)this.kRp[i]=rp[i];
- for(int i=0; i<6; i++){
- for(int j=0;j<6;j++){
- kCovariance[i][j] = covar[i][j];
- }
- }
- PkMinus= new Matrix(kCovariance,6,6);
- System.out.println(" PkMinus "); // NOW
- PkMinus.print(6,6); // NOW
- }
- /** Propagate the particle from a position x
- * to s new position in a magnetic field b
- *@param x the position
+ public int [] propagateInBarrelByDr(double rToGo, double b,int nStepsToDo, double[] matInfo){
+ int[]numSteps=new int[2];
+ setBField1Dim(b); _dedx=matInfo[0]; _radLength=matInfo[1];
+ double dr = rToGo-Math.sqrt(kRp[0]*kRp[0]+kRp[1]*kRp[1]);
+ double drStep = dr/nStepsToDo; int nStepsDone=0;
+ for (int i=1;i<=nStepsToDo;i++){
+ stepByDr(drStep);
+ nStepsDone=i;
+ if(getStopTkELow()) break;
+ }
+ numSteps[0]=nStepsToDo;numSteps[1]=nStepsDone;
+ return (numSteps);
+ }
+ /** Propagate the particle to a position r
+ * dr away from its position in a magnetic field b
+ *@param r the position
+ *@param b the magnetic field
+ *@param dedx energy loss by ionization
+ *@param materialName the name of material the
+ * particle goes through
+ */
+ public int [] propagateInBarrelByDr(double rToGo, double b,String materialName){
+ int[]numSteps = new int[2];
+ double dr = rToGo-Math.sqrt(kRp[0]*kRp[0]+kRp[1]*kRp[1]);
+ double pm = Math.sqrt(kRp[3]*kRp[3]+kRp[4]*kRp[4]+kRp[5]*kRp[5]);
+ double bStep = (pm>1.)?5.:2.; // steps5mm above pm=1GeV and 2mm below
+ int nStepsToDo= (int)(Math.abs(dr/bStep)+0.5), nStepsDone=nStepsToDo;
+ numSteps = propagateInBarrelByDr(rToGo, b, nStepsToDo,materialName);
+ return (numSteps);
+ }
+ /** Propagate the particle to a position r
+ *dr away from its position in a magnetic field b
+ *@param r the position
+ *@param b the magnetic field
+ *@param nStepsToDo the number of steps to do
+ *@param materialName the name of material the
+ * particle goes through
+ */
+ public int [] propagateInBarrelByDr(double rToGo, double b,int nStepsToDo,String materialName){
+ int[]numSteps = new int[2]; double pm ;
+ double dr = rToGo-Math.sqrt(kRp[0]*kRp[0]+kRp[1]*kRp[1]);
+ int nStepsDone = 0; double drStep = dr/nStepsToDo;
+ //setMatNam(materialName); materialByName=true; REMOVE
+ for (int i=1;i<=nStepsToDo;i++){
+ pm = Math.sqrt(kRp[3]*kRp[3]+kRp[4]*kRp[4]+kRp[5]*kRp[5]);
+ setBField1Dim(b); setDeDx(pm, kMass, materialName);setRadLength(materialName);
+ stepByDr(drStep);
+ nStepsDone=i;
+ if(getStopTkELow()) break;
+ }
+ numSteps[0]=nStepsToDo;numSteps[1]=nStepsDone;
+ return (numSteps);
+ }
+ /** Propagate the particle from a position (x,y,z)
+ * s to a new position s+ds in a magnetic field b
+ *@param ds change in path
+ *@param b the magnetic field
+ *@param nStepsToDo number of steps to do to cover ds
+ *@param materialName the name of material the particle
+ * goes through
+ * return a 2dim array with number of steps to do and
+ * number of steps done in case the track
+ * does not have enough energy to continue
+ */
+ public int[] propagateInBarrelByDs(double ds, double b, int nStepsToDo, String materialName){
+
+ int[]numSteps = new int[2];
+ double dsSteps = ds/nStepsToDo; int nStepsDone = 0;
+ for (int i=1;i<=nStepsToDo;i++) {
+ // setMatNam(materialName); materialByName=true; REMOVE
+ double pm = Math.sqrt(kRp[3]*kRp[3]+kRp[4]*kRp[4]+kRp[5]*kRp[5]);
+ setBField1Dim(b); setDeDx(pm, kMass, materialName); setRadLength(materialName);
+ stepByDs(dsSteps);
+ nStepsDone=i;
+ if(getStopTkELow()) break;
+ }
+ numSteps[0]=nStepsToDo; numSteps[1]=nStepsDone;
+ return (numSteps);
+ }
+ /** Propagate the particle from a position s
+ * to a new position s+ds in a magnetic field b
+ *@param ds change in path
+ *@param b the magnetic field
+ * assumes a step of 5mm for particles>1 GeV in magnetic field b
+ * assumes a step of 2mm for particles<=1 GeV in magnetic field b
+ *@param materialName the material encountered
+ * calls the overloaded propagateInBarrel(ds,b,nStepsToDo,materialName)
+ * return a 2dim array with number of steps to do and
+ * number of steps done in case the track
+ * does not have enough energy to continue
+ */
+ public int [] propagateInBarrelByDs(double ds, double b, String materialName){
+ int[]numSteps = new int[2];
+ double pm = Math.sqrt(kRp[3]*kRp[3]+kRp[4]*kRp[4]+kRp[5]*kRp[5]);
+ double bStep = (pm>1.)?5.:2.; // due to MagField step 5mm above pm=1GeV and 2mm below
+ int nStepsToDo= (int)(Math.abs(ds/bStep)+0.5);
+ numSteps = propagateInBarrelByDs(ds, b, nStepsToDo,materialName);
+ return (numSteps);
+ }
+ /** Propagate the particle from a position (x,y,z)
+ * s to a new position s+ds in a magnetic field b
+ *@param ds change in path
+ *@param b the magnetic field
+ * assumes a step of 5mm for particles>1 GeV in magnetic field b
+ * assumes a step of 2mm for particles<=1 GeV in magnetic field b
+ *@param dedx=materialInfo[0] energy loss by ionization
+ *@param radlength=materialInfo[1] radiation length
+ * calls the overloaded propagateInBarrel(ds,b,nStepsToDo,materialInfo);
+ * return a 2dim array with number of steps to do and
+ * number of steps done in case the track
+ * does not have enough energy to continue
+ */
+ public int[] propagateInBarrelByDs(double ds, double b, double[] materialInfo){
+ int[]numSteps = new int[2];
+ double pm = Math.sqrt(kRp[3]*kRp[3]+kRp[4]*kRp[4]+kRp[5]*kRp[5]);
+ double bStep = (pm>1.)?5.:2.; // due to MagField step 5mm above pm=1GeV and 2mm below
+ int nStepsToDo = (int)(Math.abs(ds/bStep)+0.5), nStepsDone=nStepsToDo;
+ numSteps= propagateInBarrelByDs(ds,b,nStepsToDo,materialInfo);
+ return (numSteps);
+ }
+ /**propagate the particle from a position (x,y,z)
+ * s to a new position s+ds in a magnetic field b
+ *@param ds change in path
*@param b the magnetic field
+ *@nStepsToDo number of steps to do to cover ds
+ *@param materialInfo[] of dimension 2
+ * dedx = materialInfo[0] energy loss by ionization
+ * radLength = materialInfo[1], radiation Length (cm)
+ * return a 2dim array with number of steps to do and
+ * number of steps done in case the track
+ * does not have enough energy to continue
+ */
+ public int[] propagateInBarrelByDs(double ds, double b,int nStepsToDo, double[] materialInfo){
+ int[]numSteps = new int[2];
+ setBField1Dim(b); _dedx=materialInfo[0]; _radLength=materialInfo[1];
+ int nStepsDone = 0 ; double dsSteps = ds/nStepsToDo;
+ for (int i=1;i<=nStepsToDo;i++) {
+ stepByDs(dsSteps);
+ nStepsDone=i;
+ if(getStopTkELow()) break;
+ }
+ numSteps[0]=nStepsToDo; numSteps[1]=nStepsDone;
+ return (numSteps);
+ }
+ /** Propagate Back the particlefrom a position (x,y,z)
+ * far-away to the starting position in a magnetic field b
+ *@param ds change in path
+ *@param b the magnetic field
+ *@param klb[] with 2 elements
+ * klb[0] number of steps to do
+ * klb[1] number of steps done
+ *@param materialInfo[] with 2 elements
+ *dedx =materialInfo[0] loss by ionization
+ *radlength =materialInfo[1] length it takes to interact(cm)
+ */
+ public void propagateInBarrelBackByDs(double ds, double b,double[]materialInfo , int [] klb){
+ // ds is negative
+ setBField1Dim(b); _dedx=materialInfo[0];_radLength=materialInfo[1];
+ int nStepsToDo = klb[0], nStepsDone = klb[1]; double dsSteps = ds/nStepsToDo;
+ for (int i=nStepsDone;i>0;i--) {
+ stepByDs(dsSteps);
+ }
+ }
+/** Propagate Back the particlefrom a position (x,y,z)
+ * far-away to the starting position in a magnetic field b
+ *@param ds change in path
+ *@param b the magnetic field
+ *@param klb[] with 2 elements
+ * klb[0] number of steps to do
+ * klb[1] number of steps done
+ *@param materialName the material encountered
+ */
+ public void propagateInBarrelBackByDs(double ds, double b,String materialName , int [] klb){
+ // ds is negative
+ double pm;
+ int nStepsToDo = klb[0], nStepsDone = klb[1]; double dsSteps = ds/nStepsToDo;
+ for (int i=nStepsDone;i>0;i--) {
+ pm = Math.sqrt(kRp[3]*kRp[3]+kRp[4]*kRp[4]+kRp[5]*kRp[5]);
+ setDeDx(pm, kMass, materialName); setRadLength(materialName); setBField1Dim(b);
+ stepByDs(dsSteps);
+ }
+ }
+ /** Propagate Back the particlefrom a position (x,y,z)
+ * far-away to the starting position in a magnetic field b
+ * param r radius to go to in the detector
+ * @param b the magnetic field
+ * @param klb[] with 2 elements
+ * klb[0] number of steps to do
+ * klb[1] number of steps done
+ *@param materialInfo[] with 2 elements
+ *dedx =materialInfo[0] loss by ionization
+ *radlength =materialInfo[1] length it takes to interact(cm)
+ */
+ public void propagateInBarrelBackByDr(double r, double b,double[] materialInfo, int [] klb){
+ // ds is negative
+ setBField1Dim(b); _dedx=materialInfo[0]; _radLength=materialInfo[1];
+ double drToGo = r-Math.sqrt(kRp[0]*kRp[0]+kRp[1]*kRp[1]);
+ int nStepsToDo = klb[0], nStepsDone=klb[1];
+ double drSteps = drToGo/nStepsToDo;
+ for (int i=nStepsDone;i>0;i--) {
+ stepByDr(drSteps);
+ }
+ }
+ /** Propagate Back the particlefrom a position (x,y,z)
+ * far-away to the starting position in a magnetic field b
+ * param r radius to go to in the detector
+ * @param b the magnetic field
+ * @param klb[] with 2 elements
+ * klb[0] number of steps to do
+ * klb[1] number of steps done
+ *@param materialName
*/
- public boolean propagateTo(double r, double b){
- setBField1Dim(b);
- double dr = r-Math.sqrt(kRp[0]*kRp[0]+kRp[1]*kRp[1]);
- int nsteps = (int)(Math.abs(dr/0.005)+0.5);
- for (int i=0;i<nsteps;i++) stepByDr(dr/nsteps);
- return (ktransport=true);
-}
- /** Propagate the particle from a position
- * xyz defined either by x, y, or z
- * to s new position in a magnetic field b
- *@param xyz the position
- *@param b the magnetic field
+ public void propagateInBarrelBackByDr(double r, double b,String materialName, int [] klb){
+ // ds is negative
+ double drToGo = r-Math.sqrt(kRp[0]*kRp[0]+kRp[1]*kRp[1]); double pm;
+ int nStepsToDo = klb[0], nStepsDone=klb[1]; double drSteps = drToGo/nStepsToDo;
+ for (int i=nStepsDone;i>0;i--) {
+ pm = Math.sqrt(kRp[3]*kRp[3]+kRp[4]*kRp[4]+kRp[5]*kRp[5]);
+ setBField1Dim(b); setDeDx(pm, kMass, materialName); setRadLength(materialName);
+ stepByDr(drSteps);
+ }
+ }
+ /** Propagate the particle to a position z
+ * in a magnetic field b
+ *@param z the position to go to in the detecotr
+ *@param b the magnetic field-
+ *@param materialInfo 2 dim
+ * dedx = materialInfo[0]
+ * radLength = materialInfo[1]
+ *TO BE TESTED
+ */
+ public int[] propagateToEndcaps(double z,double b,double [] materialInfo ){
+ int[]numSteps = new int[2];
+ double dz = z-kRp[2], pm= Math.sqrt(kRp[3]*kRp[3]+kRp[4]*kRp[4]+kRp[5]*kRp[5]);
+ double bStep = (pm>1.)?5.:2.; // MagField step 5mm above pm=1GeV and 2mm below
+ int nStepsToDo = (int)(Math.abs(dz/bStep)+0.5);
+ double dStep = dz/nStepsToDo;
+ numSteps = propagateToEndcaps(z,b,nStepsToDo,materialInfo);
+ return (numSteps);
+ }
+ /** Propagate the particle to a position z in a magnetic field b
+ *@param z the position to go to in the detecotr
+ *@param b the magnetic field-
+ *@param nStepsToDo number of steps to do
+ *@param materialInfo[] of dimension 2
+ * dedx = materialInfo[0] energy loss by ionization
+ * radLength = materialInfo[1], radiation Length (cm)
+ * return a 2dim array with number of steps to do & done
+ * in case the track does not have enough energy
+ *TO BE TESTED
+ */
+ public int[] propagateToEndcaps(double z,double b,int nStepsToDo, double [] materialInfo ){
+ int[]numSteps = new int[2];
+ double dz = z-kRp[2];
+ double dStep = dz/nStepsToDo; int nStepsDone=0;
+ for (int i = 1;i<=nStepsToDo;i++){
+ _dedx = materialInfo[0]; _radLength= materialInfo[1] ; setBField1Dim(b);
+ stepByDz(dz/nStepsToDo);
+ nStepsDone=i;
+ if(getStopTkELow()) break;
+ }
+ numSteps[0] = nStepsToDo; numSteps[1] = nStepsDone;
+ return (numSteps);
+ }
+ /** Propagate the particle to a position z
+ * in a magnetic field b
+ *@param z the position to go to in the detecotr
+ *@param b the magnetic field
+ *@param materialName- to calculate dedx and radLength
+ *TO BE TESTED
+ */
+ public int[] propagateToEndcaps(double z,double b,String materialName ){
+ int[]numSteps = new int[2];
+ double dz = z-kRp[2], pm= Math.sqrt(kRp[3]*kRp[3]+kRp[4]*kRp[4]+kRp[5]*kRp[5]);
+ double bStep = (pm>1.)?5.:2.; // MagField step 5mm above pm=1GeV and 2mm below
+ int nStepsToDo = (int)(Math.abs(dz/bStep)+0.5);
+ numSteps = propagateToEndcaps(z,b,nStepsToDo,materialName);
+ return (numSteps);
+ }
+ /** Propagate the particle to a position z
+ * in a magnetic field b
+ *@param z the position to go to in the detecotr
+ *@param b the magnetic field-
+ *@param nStepsToDo number of steps to do
+ *@param materialInfo[] of dimension 2
+ * dedx = materialInfo[0] energy loss by ionization
+ * radLength = materialInfo[1], radiation Length (cm)
+ * return a 2dim array with number of steps to do and
+ *TO BE TESTED
+ */
+ public int[] propagateToEndcaps(double z,double b,int nStepsToDo, String materialName){
+ int[]numSteps = new int[2];
+ double dz = z-kRp[2]; double pm;
+ double dzSteps = dz/nStepsToDo; int nStepsDone=nStepsToDo;
+ for (int i = 1;i<=nStepsToDo;i++){
+ pm= Math.sqrt(kRp[3]*kRp[3]+kRp[4]*kRp[4]+kRp[5]*kRp[5]);
+ setBField1Dim(b); setDeDx(pm, kMass, materialName); setRadLength(materialName);
+ stepByDz(dzSteps);
+ nStepsDone=i;
+ if(getStopTkELow()) break;
+ }
+ numSteps[0] = nStepsToDo; numSteps[1] = nStepsDone;
+ return (numSteps);
+ }
+ /** Propagate Back the particlefrom a position (x,y,z)
+ * far-away to the starting position in a magnetic field b
+ * param z to go to in the detector
+ * @param b the magnetic field
+ * @param klb[] with 2 elements
+ * klb[0] number of steps to do
+ * klb[1] number of steps done
+ *@param materialInfo[] with 2 elements
+ *dedx =materialInfo[0] loss by ionization
+ *radlength =materialInfo[1] length it takes to interact(cm)
+ */
+ public void propagateFromEndcapsBack(double z, double b,double[] materialInfo, int [] klb){
+ //*** dzToGo is negative
+ setBField1Dim(b); _dedx=materialInfo[0]; _radLength=materialInfo[1];
+ double dzToGo = z-kRp[2]; // negative
+ int nStepsToDo = klb[0], nStepsDone=klb[1];
+ double dzSteps = dzToGo/nStepsToDo;
+ for (int i = nStepsDone;i>0;i--) {
+ stepByDz(dzSteps);
+ }
+ }
+ /** Propagate Back the particlefrom a position (x,y,z)
+ * far-away to the starting position in a magnetic field b
+ * param r radius to go to in the detector
+ * @param b the magnetic field
+ * @param klb[] with 2 elements
+ * klb[0] number of steps to do
+ * klb[1] number of steps done
+ *@param materialName
*/
+ public void propagateFromEndCapsBack(double z, double b, String materialName, int [] klb){
+ // dzToGo is negative
+ double dzToGo = z-kRp[2]; double pm ;
+ int nStepsToDo = klb[0], nStepsDone=klb[1]; double dzSteps = dzToGo/nStepsToDo;
+ for (int i = nStepsDone;i>0;i--) {
+ //setMatNam(materialName); materialByName=true; REMOVE
+ pm = Math.sqrt(kRp[3]*kRp[3]+kRp[4]*kRp[4]+kRp[5]*kRp[5]);
+ setBField1Dim(b); setDeDx(pm, kMass, materialName); setRadLength(materialName);
+ stepByDz(dzSteps);
+ }
+ }
-public boolean propagateTo(int ixyz,double xyz, double b){
- ncount1++;
- setBField1Dim(b);
- double dxyz = xyz-kRp[ixyz];
- if(ixyz==0&& Debug)
- System.out.println("x="+xyz+" kRp[0]="+kRp[ixyz]+" ixyz="+ixyz);
- if((ixyz==0)&& (ncount1<5)&&Debug)
- System.out.println(" new x="+(kRp[0]+dxyz));
- // int nsteps = (int)(Math.abs(dxyz/0.01)+0.5); debug
- int nsteps = (int)(Math.abs(dxyz/0.2)+0.5);
- if(Debug)
- System.out.println("In propagateTo dxyz="+dxyz+" nsteps="+nsteps);
- for (int i=0;i<nsteps;i++){
- if(getStopTkELow()){ System.out.println("not enough energy left");
- break;}
- stepByDx(ixyz,dxyz/nsteps);
- if(ixyz==0&&Debug) System.out.println("!!!! In Propagate:after stepBy new x ="+kRp[0]+" step number ="+i);
- }
- return (ktransport=true);
-}
- /** stepby ds on the particle pathLength
- *@param ds one step length
- */
+ /** stepby ds on the particle pathLength-Transport Matrix
+ *@param ds one step length
+ *One of the 2 MOST IMPORTANT modules of this class ( UpDate is the second)
+ *Calculate the transport Matrix-
+ *Including Loss by Ionization DeDx
+ * Multiple Scattering
+ */
+ public boolean stepByDs(double ds){ // one step by ds
+
+ if(Debug&&ncount1<100)System.out.println("ds at entry of stepByDs"+ds);
+ double[]wk = new double[3];
+ double pIn = Math.sqrt(kRp[3]*kRp[3]+ kRp[4]*kRp[4]+kRp[5]*kRp[5]) ;
+ double E = Math.sqrt(pIn*pIn+kMass*kMass) ;
+ kB = getBField1Dim();
+ double k = kB2C*kQ*kB/pIn ;
+ double dpdx=ds*getDeDx()*(E/pIn) ;
+//*** Multiple scattering- _radLength set to 1. E-20
+ wk=vwk((getRadLength())*Math.abs(ds));
+ if(Debug&&(ncount1<100))
+ System.out.println("wk[0]="+wk[0]+" wk[1]="+wk[1]+" wk[2]="+wk[2]);
+ double[] [] valsq={{0.,0.,0.,0.,0.,0.},{0.,0.,0.,0.,0.,0.},{0.,0.,0.,0.,0.,0.},
+ {0.,0.,0.,(wk[0]*wk[0]),0.,0.},{0.,0.,0.,0.,(wk[1]*wk[1]),0.},{0.,0.,0.,0.,0.,(wk[2]*wk[2])}};
+ //*** MSk is the multiple scattering
+ Matrix MSk = new Matrix(valsq);
+ double[][]fifi= new double[6][6];
+ k*=ds;
+ if(Debug&&ncount1<100)
+ System.out.println(" k after product with ds "+k);
+ double pxx = kRp[3]; double pyy = kRp[4]; double pzz = kRp[5];
+ double xxx = kRp[0]; double yyy = kRp[1]; double zzz = kRp[2];
+ double[] pra = new double[6] ; double[] prb = new double[6];
+ pra[3]=xxx; pra[4]=yyy; pra[5]=zzz; pra[0]=pxx; pra[1]=pyy; pra[2]=pzz;
+// *** Equations of motion
+ double npxx = pxx+(k/(1.+0.25*k*k))*pyy -((0.5*k*k)/(1.0+0.25*k*k))*pxx;
+ double npyy = pyy+(-k/(1.+0.25*k*k))*pxx-((0.5*k*k)/(1.0+0.25*k*k))*pyy;
+ double npzz = pzz;
+ double mpxx = npxx - dpdx*npxx/pIn ;
+ double mpyy = npyy - dpdx*npyy/pIn ;
+ double mpzz = npzz - dpdx*npzz/pIn ;
+ double pOut = Math.sqrt(mpxx*mpxx+mpyy*mpyy+mpzz*mpzz);
+ double pMean = (pIn+pOut)/2;
+ if(dpdx > pMean){
+ System.out.println("!!!! StepByDs !!!!! : dpdx> p -track stopped not enough energy dpdx="+dpdx+" pMean="+pMean);
+ return stopTkELow=true ;
+ }
+ double twoPmean=(pMean*2);
+ kRp[3] = mpxx; kRp[4] = mpyy; kRp[5] = mpzz;
+ kRp[0] = xxx + (ds/twoPmean)*(kRp[3]+pxx);
+ kRp[1] = yyy + (ds/twoPmean)*(kRp[4]+pyy);
+ kRp[2] = zzz + (ds/twoPmean)*(kRp[5]+pzz ) ;
+// *** Equations of motions in a matriciel format = transport Matrix
+ fifi[0][0] = -((0.5*k*k)/(1+0.25*k*k))*(1.-(dpdx/pIn)) ;
+ fifi[0][1] = k/(1.+0.25*k*k)*(1.-(dpdx/pIn)) ;
+ fifi[1][0] = (-k/(1.+0.25*k*k))*(1.-(dpdx/pIn)) ;
+ fifi[1][1] = -((0.5*k*k)/(1+0.25*k*k))*(1.-(dpdx/pIn)) ;
+ fifi[2][2] = 0.;
+ fifi[3][0] = (ds/twoPmean)*(2.- ((((0.5*k*k)/(1.0+0.25*k*k)))*(1.-(dpdx/pIn)))-(dpdx/pIn));
+ fifi[3][1] = (k*ds/((1.+0.25*k*k)*twoPmean))*(1.-(dpdx/pIn)) ;
+ fifi[4][0] = (-k*ds/((1.+0.25*k*k)*twoPmean))*(1.-(dpdx/pIn)) ;
+ fifi[4][1] = (ds /twoPmean)*(2.- ((((0.5*k*k)/(1.0+0.25*k*k)))*(1.-(dpdx/pIn)))-(dpdx/pIn));
+ fifi[5][2] = (ds/twoPmean)*(2.-dpdx/pIn);
+ fi = new Matrix(fifi,6,6) ;
+ ncount1++;
+ double[] [] oldValue={{-dpdx/pIn,0.,0.,0.,0.,0.},{0.,-dpdx/pIn,0.,0.,0.,0.},{0.,0.,-dpdx/pIn,0.,0.,0.},
+ {0.,0.,0.,0.,0.,0.},{0.,0.,0.,0.,0.,0.},{0.,0.,0.,0.,0.,0.}};
+ Matrix IIdpdx =new Matrix(oldValue) ;
+ fi = (fi.plus(II)).plus(IIdpdx);
+ prb = fi.times(pra);
+//*** Check:Mtx against equations of motion
+ if(Debug&&ncount1%10==0){ //Inverted location of positions and momenta in vectors prb and in kRp
+ System.out.println("prb[3]="+prb[3]+"prb[4]="+prb[4]+"prb[5]="+prb[5]); // momentum prb[0,1,2], position prb[3,4,5]
+ System.out.println("kRp[0]="+kRp[0]+"kRp[1]="+kRp[1]+"kRp[2]="+kRp[2]); // momentum kRp[3,4,5], position kRp[0,1,2]
+ System.out.println("prb[0]="+prb[0]+"prb[1]="+prb[1]+"prb[2]="+prb[2]);
+ System.out.println("kRp[3]="+kRp[3]+"kRp[4]="+kRp[4]+"kRp[5]="+kRp[5] );
+ }
+ setNewRp(); // nwRp set to new values
+ // PkPlus has three types of value- first : set with the covariant matrix values,
+ // second : the covariant matrix transforming through phi,
+ // third : propagated covariant matrix this step
+ M6x6tmp = (fi.times(PkPlus));
+ PkMinus = (M6x6tmp.times(fi.transpose())).plus(MSk); //Uncomment to include multiple scattering
+ double ccc = 1.;
+ PkPlus = PkMinus.times(ccc);
+ kCovariance = PkPlus.getArrayCopy();
+ return (ktransport=true);
+ }
+ /** stepby dz on the z direction
+ *@param dz one step in z
+ */
+ public boolean stepByDz(double dz){ // one step by dx
+ double P = Math.sqrt(kRp[3]*kRp[3]+ kRp[4]*kRp[4]+kRp[5]*kRp[5]);
+ double dstoo = Math.abs(dz*P/kRp[(5)]);
+ return stepByDs(dstoo);
+ }
+ /** stepby dr along the radius direction
+ *@param dr one step in r
+ */
+ public boolean stepByDr(double dr){ // one step by dr
+ double P = Math.sqrt(kRp[3]*kRp[3]+kRp[4]*kRp[4]+kRp[5]*kRp[5]);
+ double ds = dsFromDr(dr);
+ if(dr<0.)System.out.println( " Negative ds stepByDr0:ds="+ds );
+ return stepByDs(ds);
+ }
+ /** ds change in path due to change in radius dr
+ *@param dr one step in r
+ */
+ public double dsFromDr(double dr){
+ // Based on Geometry and kinematics
+ double ds;
+ double drr=Math.abs(dr);
+ double pIn = Math.sqrt(kRp[3]*kRp[3]+ kRp[4]*kRp[4]+kRp[5]*kRp[5]);
+ double r = Math.sqrt(kRp[0]*kRp[0]+kRp[1]*kRp[1]);
+ double a = (kRp[3]*kRp[3]+kRp[4]*kRp[4])/(pIn*pIn);
+ double bhalf= ((kRp[0]*kRp[3]+kRp[1]*kRp[4])/pIn);
+ double c =-(drr*drr+2*drr*Math.sqrt(kRp[0]*kRp[0]+kRp[1]*kRp[1]));
+ ds=(dr>=0.)? ((-bhalf+Math.sqrt(bhalf*bhalf-a*c))/a) : (-1)*((-bhalf+Math.sqrt(bhalf*bhalf-a*c))/a);
+ //System.out.println(" stepByDr: ds="+ds);
+ return (ds);
+ }
+ /** update the phase space point and covariant Matrix
+ * This is the 2nd MOST IMPORTANT MODULE
+ *Using the measured position and errors
+ *@param mstate a 3dim array with the measure position
+ *@param mCov the measurement error Matrix
+ */
+ public boolean upDate(double []mState,double[][]mCov){
+ //*** update vector with current extrapolation
+ if(getStopTkELow()){ return (kupDate = false);}
+ double [] cState = new double[6];
+ for(int i=0;i<3;i++) {double pxyztmp= kRp[i+3]; cState[i] = pxyztmp;} // px,py,pz
+ for(int i=3;i<6;i++) {double xyztmp = kRp[i-3]; cState[i] = xyztmp;} // x,y,z
+//*** measurement covariante Matrix
+ for(int i=0;i<3;i++){ double covsm = mCov[i][i]; measCov[i][i]=covsm;}
+//*** update error covariance Matrix- Kalman Gain Matrix
+ M3x6tmp = Hk.times(PkMinus);
+ M3x3tmp = M3x6tmp.times(Hk.transpose());
+ Matrix Meas = new Matrix(measCov,3,3);
+ Matrix KkPart=M3x3tmp.plus(Meas);
+ Matrix Kk = (PkMinus.times(Hk.transpose())).times(KkPart.inverse());
+ M6x6tmp = Matrix.identity(6,6).minus(Kk.times(Hk));
+ PkPlus = M6x6tmp.times(PkMinus);
+ double []rrpa = new double[3]; double []rrpb = new double[3];
+ rrpb = Hk.times(cState); // [3][6]*[6][1] cState={px,py,pz,x,y,z)
+ for(int j=0;j<3;j++){
+ rrpb[j] = (mState[j]-rrpb[j]); //[3x1
+ }
+ double[]gainKalm=new double[6];
+ gainKalm = Kk.times(rrpb); // [6]
+ for(int i=0;i<3;i++){
+ double xyz = cState[i+3]+gainKalm[i+3] ; kRp[i] = xyz ;
+ double pxyz = cState[i]+gainKalm[i] ; kRp[i+3] = pxyz ;
+ }
+ setNewRp();
+ kCovariance=PkPlus.getArrayCopy();
+ return (kupDate = true);
+ }
+ /** dedx for a given material calculated from
+ *the Bethe-Block-(Sternheimer) calculation
+ *@param pmom (double)absolute particle momentum in GeV/c
+ *@param mass (double)particle mass in GeV
+ *@param materialName a String with the material gone through
+ *@return dedx (double)the energy loss in GeV/mm
+ *Uses the MaterialManager of Jeremy McCormick
+ *
+ */
- public boolean stepByDs(double ds){ // one step by ds
- double[]wk=new double[3];
- double P = Math.sqrt(kRp[3]*kRp[3]+ kRp[4]*kRp[4]+kRp[5]*kRp[5]) ;
- double E = Math.sqrt(P*P+kMass*kMass) ;
- double k=kB2C*kQ*kB/P ;
- double dpdx = getDeDx(P, kMass,"Iron")*E/P*ds ;
- // double dpdx=0.;
- if(Debug)
- System.out.println("I'm in StepByDs"+" k="+k);
- if(dpdx>P){
- System.out.println(" dpdx> p -track stopped not enough energy");
- System.out.println(" dpdx="+dpdx+" P="+P);
- return stopTkELow=true ;
- }
- if(ncount1<6) System.out.println("I'm in StepByDs");
- if(dpdx>P) return stopTkELow=true ;
- double rdL= getRadLength("Iron"); //im mm
- if(Debug)
- System.out.println(" number 1: dpdx="+ dpdx);
- wk=vwk((getRadLength("Iron"))*Math.abs(ds));
- if(Debug)
- System.out.println("wk[0]="+wk[0]+" wk[1]="+wk[1]+" wk[2]="+wk[2]);
- //Next to be checked again
-
- double[] [] valsq={{0.,0.,0.,0.,0.,0.},{0.,0.,0.,0.,0.,0.},{0.,0.,0.,0.,0.,0.},
- {0.,0.,0.,(wk[0]*wk[0]),0.,0.},{0.,0.,0.,0.,(wk[1]*wk[1]),0.},{0.,0.,0.,0.,0.,(wk[2]*wk[2])}};
- // MSk is the multiple scattering
- Matrix MSk = new Matrix(valsq);
- double[][]fifi=new double[6][6];
- if(Debug)
- System.out.println("1st value kRp[0]" + kRp[0]+" kRp[3]= " + kRp[3]+" kRp[4]=" + kRp[4]+ " ds="+ds);
- if(Debug)
- System.out.println(" number 2: dpdx="+ dpdx);
- kRp[0]+=(kRp[3]/P)*ds ;
- kRp[1]+=(kRp[4]/P)*ds ;
- kRp[2]+=(kRp[5]/P)*ds ;
- P=Math.sqrt(kRp[3]*kRp[3]+ kRp[4]*kRp[4]+kRp[5]*kRp[5]);
- // if(ncount1<6) System.out.println("1st value kRp[0]" + kRp[0]+" kRp[3]= " + kRp[3]+" kRp[4]=" + kRp[4]);
- //if(ncount1<6) System.out.println(" 1st value P="+ P+ " ds= "+ds);
- // System.out.println("finding out certain variables kRp[3]="+ kRp[3]+ " P="+P+" ds="+ds);
- System.out.println(" !!!!!!! HELLO FEDOR !!!!");
- k*=ds;
- double pxx = kRp[3];
- double pyy = kRp[4];
- kRp[3]+=(k/(1.+0.25*k*k/2))*pyy-0.5*k*k*pxx-dpdx*pxx/P ;
- kRp[4]+=(-k/(1.+0.25*k*k/2))*pxx-0.5*k*k*pyy-dpdx*pyy/P ;
- double p2=Math.sqrt(kRp[3]*kRp[3]+ kRp[4]*kRp[4]+kRp[5]*kRp[5]);
- if(ncount1<6) System.out.print("StepByDs ncount: "+ncount1+ " kRp[0]="+kRp[0]+" kRp[3]= " + kRp[3]+" kRp[4]=" + kRp[4]);
- kRp[5]+=-dpdx*kRp[5]/P;
- // System.out.println(" p2-P="+(p2-P)+" dpdx="+dpdx+" p2="+p2);
- fifi[0][3] =-kRp[3]*kRp[3]/(P*P*P)+1./P ;
- fifi[0][4] =-kRp[3]*kRp[4]/(P*P*P) ;
- fifi[0][5] =-kRp[3]*kRp[5]/(P*P*P) ;
- fifi[1][3] =-kRp[4]*kRp[3]/(P*P*P) ;
- fifi[1][4] =-kRp[4]*kRp[4]/(P*P*P)+1./P ;
- fifi[1][5] =-kRp[4]*kRp[5]/(P*P*P) ;
- fifi[2][3] =-kRp[5]*kRp[3]/(P*P*P) ;
- fifi[2][4] =-kRp[5]*kRp[4]/(P*P*P) ;
- fifi[2][5] =-kRp[5]*kRp[3]/(P*P*P)+1./P ;
- fifi[3][3] =-k*kRp[4]*kRp[3]/(P*P) ;
- fifi[3][4] =-k*kRp[4]*kRp[4]/(P*P)+k ;
- fifi[3][5] =-k*kRp[4]*kRp[5]/(P*P) ;
- fifi[4][3] =k*kRp[3]*kRp[3]/(P*P)-k ;
- fifi[4][4] =k*kRp[3]*kRp[4]/(P*P) ;
- fifi[4][5] =k*kRp[3]*kRp[5]/(P*P) ;
- Matrix fi = new Matrix(fifi,6,6) ;
- ncount1++;
- if ((ncount1==1)&& Debug)System.out.println("Print Matrix fi- ncount1="+ncount1);
- if(ncount1==0)fi.print(6,6);
- // fi.constructWithCopy(fifi);
-// fi=fi.times(ds); Just Checking
- Matrix fiMtx = fi.copy();
- Matrix II = Matrix.identity(6,6);
- fi=fi.plus(II);
- if(Debug){
- System.out.println("Matrix (fi*ds+I(6,6))");
- fi.print(6,6);
- }
- // update PkMinus
- if(Debug){
- System.out.println("1/Print PkMinus reenter");
- PkMinus.print(6,6);
- }
- //Pkminus =((fiMtx.times(PkMinus)).times(fiMtx.transpose())).plus(MSk); done one by one below
- if(Debug)System.out.println("2/ fi*PkMinus");
- PkMinus =((fi.times(PkMinus)));
- if(Debug)System.out.println(" fi*Pkminus");
- if(Debug)PkMinus.print(6,6);
- if(Debug)System.out.println( "3/((fi.times(PkMinus)).times(fi.transpose())) ");
- PkMinus = PkMinus.times(fi.transpose());
- if(Debug)PkMinus.print(6,6);
- if(Debug)System.out.println(" MSk");
- if(Debug)MSk.print(6,6);
- if(Debug)System.out.println("PkMinus+MSk");
- PkMinus = PkMinus.plus(MSk);
- if(Debug)PkMinus.print(6,6);
- return (ktransport=true);
-}
- /** stepby dxyz on the xyz direction
- *@param dxyz one step in xyz
- */
- private boolean stepByDx(int ixyz,double dxyz){ // one step by dx
- double P = Math.sqrt(kRp[3]*kRp[3]+ kRp[4]*kRp[4]+kRp[5]*kRp[5]);
- double dstoo=Math.abs(dxyz*P/kRp[(ixyz+3)]);
- if(Math.abs(kRp[ixyz+3])<1.E-010){
- System.out.println("Zero Pxyz, can't propagate in direction"+ixyz);
- return (stopTkELow=true);
- }
- if(dstoo<=kdistsmall)return(ktransport=true);
- return stepByDs(dstoo);
- }
+ public void setDeDx(double pmom, double mass, String materialName){
+ _dedx=1.e-20;
+ if ((Math.abs(pmom/mass )>1.e-3)&& (materialName!="")){
+ double [] p={gevTomev*kRp[3],gevTomev*kRp[4],gevTomev*kRp[5]};
+ if(Debug)System.out.println("MeV px="+p[0]+" py="+p[1]+" pz="+p[2]);
+ MaterialManager manager = MaterialManager.instance();
+ try{
+ org.lcsim.material.Material mat = manager.findMaterial(materialName);
+ double oneCm=1.;
+ _dedx = MaterialCalculator.computeBetheBloch(mat, p,(1000*mass), kQ, oneCm);// dedx MeV/cm
+ if(Debug)
+ System.out.println("getDeDx !!! dedx(MeV/cm)="+_dedx);
+ _dedx*=mevTogev/cmTomm;
+ if(Debug)
+ System.out.println("getDeDx !!! dedx(GeV/mm)="+_dedx);
+ } catch (MaterialNotFoundException e){}
+ if(_dedx<0.) _dedx=1.e-20;
+ }
+ }
+ /** set the dedx to a default value or for debug or a mixture of materials.
+ */
+ public void setDeDx(double dedx){_dedx=dedx;}
+ /** Brings the value
+ */
+ public double getDeDx(){
+ if(_dedx<0.) return (1e-20);
+ return _dedx;
+ }
+/** set the radLength using the materialCalculator of Jeremy McCormick
+ *@param materialName
+ */
+ public void setRadLength(String materialName){
+ MaterialManager manager = MaterialManager.instance();
+ double radLength=1.E-20;
+ _density=1.E-19;
+ try{
+ org.lcsim.material.Material mat = manager.findMaterial(materialName);
+ double zeff = mat.getZeff();
+ double aeff = mat.getAeff();
+ _density= mat.getDensity();
+ radLength=MaterialCalculator.computeRadiationLengthTsai(aeff,zeff); // rad Length in g/cm^-2
+ _radLength= radLength/_density ; // rad Length in cm
+ } catch (MaterialNotFoundException e){}
+ }
+ /** set the radLength and density to a default value or for debug or for a mixture
+ */
+ public void setRadLength(double xxx){_radLength=xxx;}
+ /** Brings the value
+ */
+ public double getRadLength() {
+ if(_radLength<0.) return 1.e-20;
+ return(_radLength); // in cm
+ }
+ /** Return 2 relevant informations
+ * _dedx and ;_radLength/_density
+ */
+ public double[] getMaterialInfo(){
+ double[] matParm = new double[2];
+ matParm[0]=getDeDx(); matParm[1]=getRadLength();
+ return matParm;
+ }
+ public void setMatNam(String mtn){_materialName=mtn;}
+ public String getMatNam(){return _materialName;}
+ public double getX(){return kRp[0];}
+ public double getY(){return kRp[1];}
+ public double getZ(){return kRp[2];}
-/** stepby dr along the radius direction
- *@param dr one step in r
- */
- private boolean stepByDr(double dr){ // one step by dr
- double P = Math.sqrt(kRp[3]*kRp[3]+ kRp[4]*kRp[4]+kRp[5]*kRp[5]);
- double a = kRp[3]*kRp[3]+kRp[4]*kRp[4]/P/P;
- double b = kRp[0]*kRp[3]+kRp[1]*kRp[4]/P ;
- double c =-(dr*dr+2*dr*Math.sqrt(kRp[0]*kRp[0]+kRp[1]*kRp[1]));
- double ds=(-b+Math.sqrt(b*b-a*c))/a ;
- return stepByDs(ds);
- }
-/** update the phase space point and covariant Matrix
- *Using the measured position and errors
- *@param mstate a 3dim array with the measure position
- *@param mCov the measurement error Matrix
- */
- public boolean upDate(double []mState,double[][]mCov){
- // update vector with current extrapolation
- double [] cState= new double[6];
- for(int i=0;i<6;i++) {cState[i]=kRp[i];}
- // measured position x,y,z
- measState = new double [3];
- double [][]val1 = {{1.,0.,0.,0.,0.,0.},{0.,1.,0.,0.,0.,0.},
- {0.,0.,1.,0.,0.,0.}};
- Hk=new Matrix(val1);
- if((ncount1==1)&& Debug)System.out.println( "print HK Matrix");
- if((ncount1==1)&&Debug)Hk.print(3,6);
- // measurement covariante Matrix
- for(int i=0;i<3;i++) measCov[i][i]=mCov[i][i];
- /** update error covariance Matrix- Kalman Gain Matrix
- *Problem here !!!!!!
[truncated at 1000 lines; 288 more skipped]
