lcsim/src/org/lcsim/util/step
diff -u -r1.7 -r1.8
--- TrackStepper.java 12 Dec 2005 05:46:55 -0000 1.7
+++ TrackStepper.java 18 Dec 2005 16:40:10 -0000 1.8
@@ -17,17 +17,19 @@
* 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.7 2005/12/12 05:46:55 lima Exp $
+ * @version $Id: TrackStepper.java,v 1.8 2005/12/18 16:40:10 caroline Exp $
*/
public class TrackStepper
{
+ double mm = 0.001;
+ boolean debug = false;
private static AIDA aida = AIDA.defaultInstance();
private static CalHitMapMgr evtptr = null;
// constants
- private static final double clight = 2.9979e08; // m/s
+ private static final double clight = 2.99792458e08; // 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;
+ private static final double m_default = 0.105658369; // mu mass in GeV
+ private static final double BField_default=5.0; // in Tesla
// Instance field
private TrkParams tpar = new TrkParams(0.,0.,0.,0.,0.);
private TrkPhaseSpace trakPs = new TrkPhaseSpace();
@@ -37,7 +39,7 @@
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 [] vxyz = new double[3]; // velocity in mm/s
private double dTof;
private double dL ;
@@ -85,7 +87,7 @@
*/
public void tkSteps(double r,double zmax,double[]stpcd)
{
- System.out.println("Stepper.tkSteps() called: rstep: r="+r+", zmax="+zmax);
+ if(debug) System.out.println("Stepper.tkSteps() called: rstep: r="+r+", zmax="+zmax);
// for(int i=0; i<6; ++i) {
// System.out.print("["+i+"] = )");
// }
@@ -101,7 +103,6 @@
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
@@ -111,26 +112,29 @@
materDeDx_i= stpcd[2]/d ;
// Call stepBy to do the stepping && update array[] rp
- double a_dc = d/numSteps_i ;
- System.out.println("Before "+numSteps_i+" steps particle Radius="+partR()+" Pabs="+ partPabs()); ;
+ double steplen = d/numSteps_i ;
+ if(debug) System.out.println("Before "+numSteps_i+" steps particle Radius="+partR()+" Pabs="+ partPabs()); ;
- for( int i=0; (i<numSteps_i);i++)
+ for( int i=0; i<numSteps_i; i++)
{
- rpn = stepBy(a_dc,materDeDx_i) ; // rpn each step
+ rpn = stepBy(steplen,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])>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
+ double aux = Math.sqrt(rpn[0]*rpn[0]+rpn[1]*rpn[1]);
+ atR = ( Math.abs(aux-a_r) < 1.e-3 ) ? 1 : 0; //each step
if(atZmax ==1 || stopTkElow )
{
+ if(debug) {
if(atZmax==1)System.out.println(" break For z>zMax of the Detector ");
if(stopTkElow)System.out.println("The Track is left with too little energy to continue");
- break;
+ }
+ break;
}
}
- System.out.println(" After "+numSteps_i+" steps, particle Radius="+partR()+" Pabs="+partPabs());
+ if(debug) System.out.println(" After "+numSteps_i+" steps, particle Radius="+partR()+" Pabs="+partPabs());
}
/**
* make a step along this track
@@ -154,17 +158,18 @@
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];
stopTkElow = (pAvalaible2 < dpFromMat2);
- for(int i=0;i<3;i++)
- {
- tmpxyz = rpn[i] ; tmpp = rpn[i+3];
- rpn[i] = tmpxyz + vstep[i]*dTof;
- if(!stopTkElow)
- rpn[i+3] = tmpp + dpF[i]; // - dpM[i];
- }
+ if(!stopTkElow) {
+ for(int i=0;i<3;i++)
+ {
+ tmpxyz = rpn[i] ; tmpp = rpn[i+3];
+ rpn[i] = tmpxyz + vstep[i]*dTof;
+ rpn[i+3] = tmpp + dpF[i] - dpM[i];
+ }
+ }
//System.out.println("stepBy end: rpn[0..2]=( "+rpn[0]+", "+rpn[1]+", "+rpn[2]+" )");
- String evtnum = Integer.toString( evtptr.getEvent().getEventNumber() );
+// String evtnum = Integer.toString( evtptr.getEvent().getEventNumber() );
// System.out.println("filling cloud2D, rpn: "+rpn[0]+" "+rpn[1]);
- aida.cloud2D("y vs x, evt "+evtnum).fill(rpn[0],rpn[1]);
+// aida.cloud2D("y vs x, evt "+evtnum).fill(rpn[0],rpn[1]);
// Histogram.find("py vs px").fill(rpn[3],rpn[4]);
return trakPs.upDateRp(rpn);
}
@@ -183,12 +188,12 @@
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*(-BField_i)*dTof;
- //dpField[1]=-a_q*0.3*(rpn[3]/a_En)* clight*(-BField_i)*dTof;
+ //dpField[0]=a_q*(clight*1e-9)*(rpn[4]/a_En) * clight*(-BField_i)*dTof;
+ //dpField[1]=-a_q*(clight*1e-9)*(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);
+ alpha = a_q*(clight*1e-9)*(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];
@@ -249,7 +254,7 @@
this.tpar =tkprm ;
df.setMaximumFractionDigits(5);
BField_i=(BField_i==-10.)?BField_default:BField_i;
- System.out.println("stpr:rpFromTrk:BField="+BField_i);
+ if(debug) System.out.println("stpr:rpFromTrk:BField="+BField_i);
trakPs.trkRpFromParams(tpar,inMass, BField_i);
return trakPs.getRp();
}
