lcsim/src/org/lcsim/util/step
diff -u -r1.2 -r1.3
--- TrackStepper.java 27 Jun 2005 21:32:14 -0000 1.2
+++ TrackStepper.java 31 Oct 2005 22:33:42 -0000 1.3
@@ -14,7 +14,7 @@
* Particle assumed is a Muon, but can be changed by setting rp[6] at input
* Units are Tesla, cm, GeV/c
* @author C. Milstene.
- * @version $Id: TrackStepper.java,v 1.2 2005/06/27 21:32:14 lima Exp $
+ * @version $Id: TrackStepper.java,v 1.3 2005/10/31 22:33:42 caroline Exp $
*/
public class TrackStepper
{
@@ -22,41 +22,43 @@
private static final double clight = 2.9979e08; // m/s
private static final double halfpi = 0.5 * Math.PI;
private static final double m_default = 0.106 ;
+ private static final double BField_default=5.0;
// Instance field
private TrkParams tpar = new TrkParams(0.,0.,0.,0.,0.);
private TrkPhaseSpace trakPs = new TrkPhaseSpace();
private DecimalFormat df = new DecimalFormat();
- private double [] rpn= new double[8];
- private double [] r_ints = new double[3];
- private double [] p_ints = new double[3];
+ protected double [] rpn= new double[8];
+ protected double [] r_ints = new double[3];
+ protected double [] p_ints = new double[3];
private double [] dpField = new double[3];
private double [] dpFordL = new double[3];
private double [] vxyz = new double[3];
private double dTof;
- private double dL ;
+ private double dL ;
double []xy = new double[3];
private double []dpF = new double[3];
- private double []dpM = new double[3];
+ private double []dpM = new double[3];
private double []vstep = new double[3];
-
-
-
-
+
// Status Flags
private boolean stopTkElow = false;
private int stepFlags = 0 ;
private int atZmax = 0 ;
private int atR = 0 ;
// Material Information
- private int a_numSteps;
- private double a_BField;
- private double a_materDeDx ;
+ private int numSteps_i;
+ protected double BField_i =-10. ;
+ private double materDeDx_i ;
+
// private Detector detNow ;
// Multiple Local use
private double tmpxyz ;
private double tmpp ;
- private double a_mass ;
+ private double mass_i ;
+ private double particleR;
+ private double particlePt;
+ private double particlePabs;
/** Default Constructor */
public TrackStepper()
@@ -69,7 +71,7 @@
* both cases return with old rp values. Otherwise continue
* calling stepBy to step until the total number of steps
* calculated to reach the r boundary is exhausted. Allow a
- * Delta_r = 0.01 for round_off differences in r.
+ * Delta_r = 0.001mm for round_off differences in r.
*
* @param r the radial boundary to aim for- generally r next layer
* @param zmax the z boundary of the detector
@@ -81,30 +83,40 @@
// for(int i=0; i<6; ++i) {
// System.out.print("["+i+"] = )");
// }
-
+ // put all values used in comparison to calculations in cm
+ // put all values returned in mm
// put values of x,y,z,px,py,pz,q,amass in rp
// Check if at boundary at entry, if so return with old values
+
+ rpn=trakPs.getRp(); // get phaseSpace at entry
+
double a_r = r;
- //double [] rpn= new double[8];
- rpn=trakPs.getRp(); // get phaseSpace at entry
- a_mass = (rpn[6] ==-1.)?m_default:rpn[6];
- atZmax = (Math.abs(rpn[2])>zmax)? 1 :0;
- atR = (Math.sqrt(rpn[0]*rpn[0]+rpn[1]*rpn[1])>(a_r+0.01))?1:0;
+ double a_zmax=zmax;
+ mass_i = (rpn[6] ==-1.)?m_default:rpn[6];
+ atZmax = (Math.abs(rpn[2])>a_zmax)? 1 :0;
+ atR = (Math.sqrt(rpn[0]*rpn[0]+rpn[1]*rpn[1])>(a_r+0.001))?1:0; // add 0.001 in mm now-Sept05
+
if(atZmax==1 || stopTkElow) return;
+
// access values of dE/dx & number of steps from stpcd
- a_numSteps =(int) stpcd[0];
- a_BField = stpcd[1];
+ numSteps_i =(int) stpcd[0]; BField_i = stpcd[1];
double d=Math.abs(Math.sqrt(rpn[0]*rpn[0]+rpn[1]*rpn[1])-a_r) ;
- a_materDeDx= stpcd[2]/d ;
+ // Check if distance cm or distance mm !!!!
+ materDeDx_i= stpcd[2]/d ;
+
// Call stepBy to do the stepping && update array[] rp
- double a_dc = d/a_numSteps ;
- for( int i=0; (i<a_numSteps);i++)
+ double a_dc = d/numSteps_i ;
+ System.out.println("Before "+numSteps_i+" steps particle Radius="+partR()+" Pabs="+ partPabs()); ;
+
+ for( int i=0; (i<numSteps_i);i++)
{
- rpn = stepBy(a_dc,a_materDeDx) ; // rpn each step
- System.out.println("rpn[0..2]=( "+rpn[0]+", "+rpn[1]+", "+rpn[2]+" )");
+ rpn = stepBy(a_dc,materDeDx_i) ; // rpn each step
+ // System.out.println("Stpr:rpn[0..2]=( "+rpn[0]+", "+rpn[1]+", "+rpn[2]+" )");
+ // System.out.println("stpr:rpn[3..5]=( "+rpn[3]+", "+rpn[4]+", "+rpn[5]+" )");
+
// Check boundaries after each step
- atZmax = (Math.abs(rpn[2])>zmax)? 1 :0; //each step
- atR = (Math.sqrt(rpn[0]*rpn[0]+rpn[1]*rpn[1])>(a_r+0.01))?1:0; //each step
+ atZmax = (Math.abs(rpn[2])>a_zmax)? 1 :0; //each step
+ atR = (Math.sqrt(rpn[0]*rpn[0]+rpn[1]*rpn[1])>(a_r+0.001))?1:0; //each step
if(atZmax ==1 || stopTkElow )
{
if(atZmax==1)System.out.println(" break For z>zMax of the Detector ");
@@ -112,6 +124,7 @@
break;
}
}
+ System.out.println(" After "+numSteps_i+" steps particleRadius="+partR()+" Pabs="+partPabs());
}
/**
* make a step along this track
@@ -119,8 +132,7 @@
*/
public double [] stepBy(double dc,double matDeDx)
{
- rpn=trakPs.getRp();
- System.out.println("stepBy: rpn[0..2]=( "+rpn[0]+", "+rpn[1]+", "+rpn[2]+" )");
+ //System.out.println("stepBy: rpn[0..2]=( "+rpn[0]+", "+rpn[1]+", "+rpn[2]+" )");
setVelocity();
tmpxyz=0.; tmpp=0.;
setDTOF(dc);
@@ -129,6 +141,9 @@
dL = getPathLength();
dpF = dpFromField(dc);
dpM = dpFromMatter(dc,matDeDx);
+
+ //System.out.println(" !!!=> px,dpF0,dpM0="+rpn[3]+" "+dpF[0]+" "+dpM[0]);
+ //System.out.println("!!!=> py,dpF1,dpM1="+rpn[4]+" "+dpF[1]+" "+dpM[1]);
vstep = getVelocity();
double dpFromMat2 = dpM[0]*dpM[0]+dpM[1]*dpM[1]+dpM[2]*dpM[2];
double pAvalaible2= rpn[3]*rpn[3]+rpn[4]*rpn[4]+rpn[5]*rpn[5];
@@ -138,11 +153,9 @@
tmpxyz = rpn[i] ; tmpp = rpn[i+3];
rpn[i] = tmpxyz + vstep[i]*dTof;
if(!stopTkElow)
- rpn[i+3] = tmpp + dpF[i] - dpM[i];
- else
- rpn[i+3] = 0.;
+ rpn[i+3] = tmpp + dpF[i]; // - dpM[i];
}
- System.out.println("stepBy end: rpn[0..2]=( "+rpn[0]+", "+rpn[1]+", "+rpn[2]+" )");
+ //System.out.println("stepBy end: rpn[0..2]=( "+rpn[0]+", "+rpn[1]+", "+rpn[2]+" )");
// Histogram.find("y vs x").fill(rpn[0],rpn[1]);
// Histogram.find("py vs px").fill(rpn[3],rpn[4]);
return trakPs.upDateRp(rpn);
@@ -153,19 +166,25 @@
* system. (Inverted sign of curvatures.) by changing the sign of the
* magnetic field to -Bz
*/
+
private double[] dpFromField(double dc)
{
-
- rpn = trakPs.getRp();
+ double alpha; double delta; double adjust; double mainTerm;
double pabs2 = rpn[3]*rpn[3]+rpn[4]*rpn[4]+rpn[5]*rpn[5];
- a_mass = rpn[6];
- double a_En = Math.sqrt(pabs2 + a_mass*a_mass);
+ mass_i = rpn[6];
+ double a_En = Math.sqrt(pabs2 + mass_i*mass_i);
double a_q = rpn[7];
- //dpField[0]=a_q*0.3*(rpn[4]/a_En) * clight*(-a_BField)*getDTOF(dc);
- //dpField[1]=-a_q*0.3*(rpn[3]/a_En)* clight*(-a_BField)*getDTOF(dc);
- dpField[0]=a_q*0.3*(rpn[4]/a_En) * clight*(-a_BField)*dTof;
- dpField[1]=-a_q*0.3*(rpn[3]/a_En)* clight*(-a_BField)*dTof;
+ //dpField[0]=a_q*0.3*(rpn[4]/a_En) * clight*(-BField_i)*dTof;
+ //dpField[1]=-a_q*0.3*(rpn[3]/a_En)* clight*(-BField_i)*dTof;
+
+ assert BField_i > 4.0 : "BField too small!";
+
+ alpha = a_q*0.3*(1./a_En)*(clight*dTof);
+ delta =1.+0.25*alpha*alpha*BField_i*BField_i; mainTerm=alpha/delta;
+ adjust = 0.5*(1./delta)*alpha*alpha*BField_i*BField_i;
+ dpField[0]=mainTerm*rpn[4]*(-BField_i) - adjust*rpn[3];
+ dpField[1]=-mainTerm*rpn[3]*(-BField_i) - adjust*rpn[4];
dpField[2]=0.;
return dpField;
}
@@ -174,17 +193,16 @@
private double[] dpFromMatter(double dc,double materMeanDeDx)
{
- rpn=trakPs.getRp();
double pabs2= rpn[3]*rpn[3]+rpn[4]*rpn[4]+rpn[5]*rpn[5];
- a_materDeDx = materMeanDeDx;
+ materDeDx_i = materMeanDeDx;
double tempdpoverdl;
- double a_mass = rpn[6];
- double a_En = Math.sqrt(pabs2 + a_mass*a_mass);
+ double mass_i = rpn[6];
+ double a_En = Math.sqrt(pabs2 + mass_i*mass_i);
for (int i =0; i<3;i++)
{
// Change in Momentum due to material
- tempdpoverdl = (a_materDeDx*a_En*(rpn[i+3]/pabs2));
+ tempdpoverdl = (materDeDx_i*a_En*(rpn[i+3]/pabs2));
dpFordL[i] = tempdpoverdl*dL;
}
@@ -192,7 +210,7 @@
}
/**
- * Return vector rp(x,y,z,px,py,pz,a_mass=MuonMass,charge,BField)
+ * Return vector rp(x,y,z,px,py,pz,mass_i=MuonMass,charge,BField)
* from the information at the Interaction Point given by TrkParams.
*
* @param tkprm TrkParams at Interaction Point (IP)
@@ -203,12 +221,13 @@
this.tpar =tkprm ;
df.setMaximumFractionDigits(5);
double inM=m_default;
- trakPs.trkRpFromParams(tpar,inM, a_BField);
+ BField_i=(BField_i==-10.)?BField_default:BField_i;
+ trakPs.trkRpFromParams(tpar,inM, BField_i);
return trakPs.getRp();
}
/**
- * Return vector rp(x,y,z,px,py,pz,a_mass=inMass,charge,BField)
+ * Return vector rp(x,y,z,px,py,pz,mass_i=inMass,charge,BField)
* from the information at the Interaction Point given by TrkParams.
* @param tkprm TrkParams at Interaction Point(IP
* @param inMass
@@ -218,7 +237,9 @@
// detNow= Detector.instance();
this.tpar =tkprm ;
df.setMaximumFractionDigits(5);
- trakPs.trkRpFromParams(tpar,inMass, a_BField);
+ BField_i=(BField_i==-10.)?BField_default:BField_i;
+ System.out.println("stpr:rpFromTrk:BField="+BField_i);
+ trakPs.trkRpFromParams(tpar,inMass, BField_i);
return trakPs.getRp();
}
@@ -227,20 +248,19 @@
public void setVelocity()
{
double tempV ;
- rpn = trakPs.getRp();
- a_mass=(rpn[6]==-1.)?m_default:rpn[6] ; // Default Muon mass
- double En = Math.sqrt(rpn[3]*rpn[3]+rpn[4]*rpn[4]+rpn[5]*rpn[5]+a_mass*a_mass);
+ //rpn = trakPs.getRp();
+ mass_i=(rpn[6]==-1.)?m_default:rpn[6] ; // Default Muon mass
+ double En = Math.sqrt(rpn[3]*rpn[3]+rpn[4]*rpn[4]+rpn[5]*rpn[5]+mass_i*mass_i);
for(int i=0;i<3;i++)
{
- tempV = (rpn[i+3]/En)*clight*100. ;
+ tempV = (rpn[i+3]/En)*clight*1000. ;
vxyz[i]=tempV ;
}
-
- }
+ }
/**
* Return velocity vx,vy,vz of the particle this step
*/
- public double []getVelocity()
+ public double[] getVelocity()
{
return vxyz;
}
@@ -254,8 +274,7 @@
double a=0.;double b=0.; double c=0.;
double a_dc = dc ;
- //double dTof;
- xy = trakPs.getPosition();
+ xy[0]=rpn[0]; xy[1]=rpn[1];xy[2]=rpn[2];
vstep=getVelocity();
for(int i=0;i<2;i++)
{
@@ -293,21 +312,48 @@
*/
public double[] getRAfterStep()
{
- //double [] r_ints = new double[3];
+
return r_ints =trakPs.getPosition();
}
/** @return Momentum (px,py,pz) after step
*/
public double[] getPAfterStep()
- {
- //double [] p_ints = new double[3];
- return p_ints=trakPs.getMomentum();
+ {
+ return p_ints=trakPs.getMomentum();
}
+ /** @return phase-space point after step
+ */
public double [] getNewRp()
{
return rpn= trakPs.getRp();
}
+ /** @return Radius at which the particle is
+ */
+ public double partR()
+ {
+ particleR= rpn[0]*rpn[0]+rpn[1]*rpn[1];
+ if(particleR!=0.)particleR = Math.sqrt(particleR);
+ return particleR;
+ }
+ /** @return The Pt of the particle
+ */
+ public double partPt()
+ {
+
+ particlePt= rpn[3]*rpn[3]+rpn[4]*rpn[4];
+ if(particlePt!=0.) particlePt=Math.sqrt(particlePt);
+ return particlePt;
+ }
+ /** @return the absolute momentum pabs
+ */
+ public double partPabs()
+ {
+
+ particlePabs= rpn[3]*rpn[3]+rpn[4]*rpn[4]+rpn[5]*rpn[5];
+ if(particlePabs!=0.) particlePabs= Math.sqrt(particlePabs);
+ return particlePabs;
+ }
/** Reset the Stepping flags by calling setStopTkELow()
setAtZmax(), setCurlBack, and setting stepFlags to zero
*/
@@ -315,7 +361,7 @@
{
resetStopTkELow();
resetAtZmax() ;
-
+
int stepFlags = 0;
}
@@ -335,7 +381,7 @@
/** Set Boolean Indicator If Energy Left
*/
public void setStopTkELow(boolean b)
- {boolean stopTkElow=b;}
+ { stopTkElow=b;}
/** Reset Boolean Indicator Energy Left for the step
*/
public void resetStopTkELow()
@@ -347,7 +393,7 @@
/** Set integer to 1(track out of detector),0 otherwise
*/
public void setAtZmax(int iOneZero)
- { int atZmax=iOneZero;}
+ { atZmax=iOneZero;}
/** Reset integer to 0, hitZmax=1 when track out of the detector.
*/
public void resetAtZmax()
@@ -359,7 +405,7 @@
/* Set integer to 1(track at R Boundary),0 otherwise
*/
public void setAtR( int inOneZero)
- {int atR=inOneZero;}
+ {atR=inOneZero;}
/** Reset integer to 0, hitR=1 when track reach Boundary in z
*/
public void resetAtR()
@@ -380,6 +426,8 @@
*/
public void clearTrackStepper()
{
+ trakPs.clear();
+ BField_i = -10.;
setDTOF(0.);
setPathLength(0.);
dpFromField(0.);