lcsim/src/org/lcsim/util/step
diff -u -r1.11 -r1.12
--- TrackStepper.java 16 Jun 2006 20:45:36 -0000 1.11
+++ TrackStepper.java 13 Jul 2006 23:35:11 -0000 1.12
@@ -4,6 +4,7 @@
import java.text.*;
import org.lcsim.util.aida.*;
import org.lcsim.recon.cluster.util.CalHitMapMgr;
+import org.lcsim.geometry.Detector;
/**Given a particle of initial position momentum and charge
* (x,y,z,px,py,pz,and q),calculates its trajectory in a uniform
@@ -17,9 +18,14 @@
* 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.11 2006/06/16 20:45:36 caroline Exp $
+ * @version $Id: TrackStepper.java,v 1.12 2006/07/13 23:35:11 caroline Exp $
*NEXT To BE DONE: get subdet/layer info diretly
*get materiel info from a util/step/stepperCondition.class
+ *
+ New modifications by Judith Odili and Caroline Milstene - 07/13/06 -
+ Modified so the Stepper information get read
+ directly from StepConditions class.
+
*/
public class TrackStepper
{
@@ -30,7 +36,7 @@
// constants
//private static final double clight = 2.99792458e08; // m/s
private static final double clight = 3.0e08; // m/s
- private static final double halfpi = 0.5 * Math.PI;
+ // private static final double halfpi = 0.5 * Math.PI;
private static final double m_default = 0.105658369; // mu mass in GeV
private static final double BField_default=5.0; // in Tesla
private static double dT0=1.0E-13;
@@ -77,6 +83,7 @@
private double particlePabs ;
private double sumdL ;
private double sumdpFM ;
+ public StepConditions cond, condMat;
/** Default Constructor */
public TrackStepper()
@@ -99,7 +106,9 @@
* @param stpcd
* an array ={xNumberOfSteps, BField, Mean_dE/dx_this layer}
*/
- public void tkSteps(double r,double zmax,double[]stpcd)
+
+
+ public void tkSteps(double r,double zmax,double[]stpcd)
{
if(debug) System.out.println("Stepper.tkSteps() called: rstep: r="+r+", zmax="+zmax);
@@ -123,9 +132,7 @@
drLength =(dThick>=1200.)?6.:((dThick/10)*pFactor);
materDeDx_i=(stpcd[2]/stpcd[0]) ; //mean per mmm multiplied by dL-in dPFromMatter
numSteps_i =(dThick>=1200.)?210:(int)(dThick/drLength);
- System.out.println("dThick="+dThick+" rpn-012="+rpn[0]+" "+rpn[1]+" "+rpn[2]);
- System.out.print("drLength="+drLength+"numSteps_i="+numSteps_i);
- System.out.println("materDeDx_i"+materDeDx_i);
+
setDTOF(rpn,drLength);
vv=getVelocity();
rv=((rpn[0]*vv[0]+rpn[1]*vv[1]))/(partR(rpn)*Math.sqrt(vv[0]*vv[0]+vv[1]*vv[1]));
@@ -136,24 +143,185 @@
if(atZmax==1 || stopTkElow) return;
// access values of dE/dx & number of steps from stpcd
// Check if distance cm or distance mm !!!!
- // Call stepBy to do the stepping && update array[] rp
- System.out.println("Before "+numSteps_i+" steps:Px Py x y Radius (r.v)" );
- System.out.println(rpn[3]+" "+rpn[4]+" "+rpn[0]+" "+rpn[1]+" "+partR(rpn)+" "+rv);
+
if(partR(rpn)<2.)dT0=getDTOF();
- System.out.println("Before bigloop DTOF="+getDTOF());
+
int nLoop=0;
atR=0;
- System.out.print("Loop Entry:atZmax=" +atZmax+" stopTkElow="+stopTkElow);
- System.out.println(" curlBack="+curlBack+" atR="+atR);
- System.out.println(" RLIMIT="+a_r);
+
bigLoop:while(((atZmax!=1)&&!stopTkElow)&&((nLoop<=10)||((a_r>1000)&&(nLoop<=1270)))){
nLoop++;
- if(nLoop==1)System.out.println("HERE: numSteps_i"+numSteps_i+"drLength="+drLength+"R="+partR(rpn));
+
if(partR(rpn)==0)break;
- System.out.println(" ENTRY OF BIGLOOP nLoop="+nLoop+" drLength="+drLength+" dt="+getDTOF()+" dTof="+dTof);
- if(nLoop==1)System.out.println(" ENTRY OF FOR LOOP numsteps_i="+numSteps_i);
- rpn = stepBy(drLength,materDeDx_i) ; // rpn each step
+ rpn = stepBy(drLength,materDeDx_i) ; // rpn each step
+ pFactor=(partPabs(rpn)>=1.)?1.:(5*partPabs(rpn));
+ setVelocity(rpn); vv=getVelocity();
+ rv = (rpn[0]*vv[0]+rpn[1]*vv[1])/(partR(rpn)*Math.sqrt(vv[0]*vv[0]+vv[1]*vv[1]));
+ setDTOF(rpn,drLength);
+ aida.cloud2D("behavior of rv in loop").fill(getDTOF(),rv);
+ pFactor=(partPabs(rpn)>=1.)?1.:(5*partPabs(rpn));
+ setCurlBack(rv<=0.35);
+ if((rv<0.35) && (Math.abs(rv)>0.35)&&(curlBackFlag==1)){
+ curlBackFlag=2; break;
+ }
+ if(( (curlBackFlag==0) && (rv<=0.35) )||(getDTOF()>1.5*dT0))
+ {
+ if(rv<=0.35)curlBackFlag=1 ;
+ drLength =(dThick>1200.)?6.:((dThick/100)*pFactor);
+ numSteps_i =(dThick>1200.)?210:(int)(dThick/drLength);
+ }
+ // Check boundaries after each step
+ atZmax = (Math.abs(rpn[2])>(a_zmax+(drLength/10)))? 1 :0; //each step
+ atR = ( Math.abs(partR(rpn)-a_r) < drLength) ? 1 : 0; //each step
+ if((atZmax==1)||stopTkElow||(atR==1)) break;
+ }
+ if(atZmax ==1 || stopTkElow ||curlBack)
+ {
+ if((atZmax==1)||stopTkElow){
+ 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");
+ }else if(curlBack){
+ System.out.println(" The track curls back go to a smaller radius");
+ }
+ }
+ double pabsExit = partPabs(rpn);
+ double dpFromDedxLay =sumdpFM;
+
+ rv=((rpn[0]*vv[0]+rpn[1]*vv[1]))/(partR(rpn)*Math.sqrt(vv[0]*vv[0]+vv[1]*vv[1]));
+ aida.cloud1D("Evolution of scalar productrv ").fill(rv);
+ setDTOF(rpn,drLength);
+
+ }
+
+
+ public void tkSteps(double r,double zmax,Detector det, String materialname)
+ {
+
+ double[]stpcd = new double[3];
+ double[]vv=new double[3];
+ double rv; double drLength;
+ rpn=trakPs.getRp(); // get phaseSpace at entry
+ setVelocity(rpn) ; // particle velocity
+
+ condMat = new StepConditions(rpn,det, materialname, r);
+
+ double pabsEntry =partPabs(rpn) ;
+ a_r = r;
+ a_zmax=zmax;
+
+ dThick = stpcd[0] = condMat.getMaterialThickness();
+ double []b = new double[3];
+ b = condMat.getField();
+ BField_i = stpcd[1] =b[2] ;
+ stpcd[2]= condMat.getMatterDeDx();
+ sumdL =0. ;
+ sumdpFM=0. ;
+ mass_i = (rpn[6] ==-1.)?m_default:rpn[6];
+ pFactor=(partPabs(rpn)>=1.)?1.:(5*partPabs(rpn));
+ drLength =(dThick>=1200.)?6.:((dThick/10)*pFactor);
+ materDeDx_i=(stpcd[2]/stpcd[0]) ; //mean per mmm multiplied by dL-in dPFromMatter
+ numSteps_i =(dThick>=1200.)?210:(int)(dThick/drLength);
+
+ setDTOF(rpn,drLength);
+ vv=getVelocity();
+ rv=((rpn[0]*vv[0]+rpn[1]*vv[1]))/(partR(rpn)*Math.sqrt(vv[0]*vv[0]+vv[1]*vv[1]));
+ double vv1=Math.sqrt(vv[0]*vv[0]+vv[1]*vv[1]);
+ aida.cloud1D("Evolution vradial ").fill(vv1);
+ atZmax = (Math.abs(rpn[2])>(a_zmax+(drLength/10)) )? 1 :0;
+ atR = ((Math.sqrt(rpn[0]*rpn[0]+rpn[1]*rpn[1])-(a_r))<(drLength))?1:0; // add 0.001 in mm now-Sept05
+ if(atZmax==1 || stopTkElow) return;
+
+ if(partR(rpn)<2.)dT0=getDTOF();
+ int nLoop=0;
+ atR=0;
+ bigLoop:while(((atZmax!=1)&&!stopTkElow)&&((nLoop<=10)||((a_r>1000)&&(nLoop<=1270)))){
+ nLoop++;
+ if(partR(rpn)==0)break;
+ rpn = stepBy(drLength,materDeDx_i) ; // rpn each step
+ pFactor=(partPabs(rpn)>=1.)?1.:(5*partPabs(rpn));
+ setVelocity(rpn); vv=getVelocity();
+ rv = (rpn[0]*vv[0]+rpn[1]*vv[1])/(partR(rpn)*Math.sqrt(vv[0]*vv[0]+vv[1]*vv[1]));
+ setDTOF(rpn,drLength);
+ aida.cloud2D("behavior of rv in loop").fill(getDTOF(),rv);
+ pFactor=(partPabs(rpn)>=1.)?1.:(5*partPabs(rpn));
+ setCurlBack(rv<=0.35);
+ if((rv<0.35) && (Math.abs(rv)>0.35)&&(curlBackFlag==1)){
+ curlBackFlag=2; break;
+ }
+ if(( (curlBackFlag==0) && (rv<=0.35) )||(getDTOF()>1.5*dT0))
+ {
+ if(rv<=0.35)curlBackFlag=1 ;
+ drLength =(dThick>1200.)?6.:((dThick/100)*pFactor);
+ numSteps_i =(dThick>1200.)?210:(int)(dThick/drLength);
+ }
+ // Check boundaries after each step
+ atZmax = (Math.abs(rpn[2])>(a_zmax+(drLength/10)))? 1 :0; //each step
+ atR = ( Math.abs(partR(rpn)-a_r) < drLength) ? 1 : 0; //each step
+ if((atZmax==1)||stopTkElow||(atR==1)) break;
+ }
+ if(atZmax ==1 || stopTkElow ||curlBack)
+ {
+ if((atZmax==1)||stopTkElow){
+ 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");
+ }else if(curlBack){
+ System.out.println(" The track curls back go to a smaller radius");
+ }
+ }
+ double pabsExit = partPabs(rpn);
+ double dpFromDedxLay =sumdpFM;
+
+ rv=((rpn[0]*vv[0]+rpn[1]*vv[1]))/(partR(rpn)*Math.sqrt(vv[0]*vv[0]+vv[1]*vv[1]));
+ aida.cloud1D("Evolution of scalar productrv ").fill(rv);
+ setDTOF(rpn,drLength);
+
+ }
+
+ public void tkSteps(double r,double zmax, Detector det,String subdetName, int layerN)
+ {
+
+ double[]stpcd = new double[3];
+ double[]vv=new double[3];
+ double rv; double drLength;
+ rpn=trakPs.getRp(); // get phaseSpace at entry
+ setVelocity(rpn) ; // particle velocity
+
+ cond = new StepConditions(r, rpn, det, subdetName, layerN);
+
+ double pabsEntry =partPabs(rpn) ;
+ a_r = r;
+ a_zmax=zmax;
+ dThick = stpcd[0] = cond.getLayerThickness();
+ double []b = new double[3];
+ b = cond.getField();
+ BField_i = stpcd[1] =b[2] ;
+ stpcd[2] = cond.getlayerDeDx();
+
+ sumdL =0. ;
+ sumdpFM=0. ;
+ mass_i = (rpn[6] ==-1.)?m_default:rpn[6];
+ pFactor=(partPabs(rpn)>=1.)?1.:(5*partPabs(rpn));
+ drLength =(dThick>=1200.)?6.:((dThick/10)*pFactor);
+ materDeDx_i=(stpcd[2]/stpcd[0]) ; //mean per mmm multiplied by dL-in dPFromMatter
+ numSteps_i =(dThick>=1200.)?210:(int)(dThick/drLength);
+
+ setDTOF(rpn,drLength);
+ vv=getVelocity();
+ rv=((rpn[0]*vv[0]+rpn[1]*vv[1]))/(partR(rpn)*Math.sqrt(vv[0]*vv[0]+vv[1]*vv[1]));
+ double vv1=Math.sqrt(vv[0]*vv[0]+vv[1]*vv[1]);
+ aida.cloud1D("Evolution vradial ").fill(vv1);
+ atZmax = (Math.abs(rpn[2])>(a_zmax+(drLength/10)) )? 1 :0;
+ atR = ((Math.sqrt(rpn[0]*rpn[0]+rpn[1]*rpn[1])-(a_r))<(drLength))?1:0; // add 0.001 in mm now-Sept05
+ if(atZmax==1 || stopTkElow) return;
+
+ if(partR(rpn)<2.)dT0=getDTOF();
+ int nLoop=0;
+ atR=0;
+ bigLoop:while(((atZmax!=1)&&!stopTkElow)&&((nLoop<=10)||((a_r>1000)&&(nLoop<=1270)))){
+ nLoop++;
+ if(partR(rpn)==0)break;
+ rpn = stepBy(drLength,materDeDx_i) ; // rpn each step
pFactor=(partPabs(rpn)>=1.)?1.:(5*partPabs(rpn));
setVelocity(rpn); vv=getVelocity();
rv = (rpn[0]*vv[0]+rpn[1]*vv[1])/(partR(rpn)*Math.sqrt(vv[0]*vv[0]+vv[1]*vv[1]));
@@ -186,16 +354,86 @@
}
double pabsExit = partPabs(rpn);
double dpFromDedxLay =sumdpFM;
-// {
- //System.out.print("meanDeDx="+stpcd[2]+" dpFromDedxLay="+dpFromDedxLay+" sumdL="+sumdL);
- //System.out.println(" pEntry"+pabsEntry+" pExit="+pabsExit+" pExit-PEntry="+(pabsExit-pabsEntry));
-// }
- System.out.println("After "+numSteps_i+"Px Py x y En Radius r.v DTOF" );
+
rv=((rpn[0]*vv[0]+rpn[1]*vv[1]))/(partR(rpn)*Math.sqrt(vv[0]*vv[0]+vv[1]*vv[1]));
aida.cloud1D("Evolution of scalar productrv ").fill(rv);
setDTOF(rpn,drLength);
- System.out.print(rpn[3]+" "+rpn[4]+" "+rpn[0]+" "+rpn[1]+" "+partR(rpn)+" "+rv);
- System.out.println(getDTOF());
+
+ }
+ private void Stepping (double r, double zmax, double stpcd[])
+ {
+ if(debug) System.out.println("Stepper.tkSteps() called: rstep: r="+r+", zmax="+zmax);
+ // 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
+ double[]vv=new double[3];
+ double rv; double drLength;
+ //rpn=trakPs.getRp(); // get phaseSpace at entry
+ setVelocity(rpn) ; // particle velocity
+ double pabsEntry =partPabs(rpn) ;
+ a_r = r;
+ a_zmax=zmax;
+ sumdL =0. ;
+ sumdpFM=0. ;
+ mass_i = (rpn[6] ==-1.)?m_default:rpn[6];
+ pFactor=(partPabs(rpn)>=1.)?1.:(5*partPabs(rpn));
+ drLength =(dThick>=1200.)?6.:((dThick/10)*pFactor);
+ materDeDx_i=(stpcd[2]/stpcd[0]) ; //mean per mmm multiplied by dL-in dPFromMatter
+ numSteps_i =(dThick>=1200.)?210:(int)(dThick/drLength);
+ setDTOF(rpn,drLength);
+ vv=getVelocity();
+ rv=((rpn[0]*vv[0]+rpn[1]*vv[1]))/(partR(rpn)*Math.sqrt(vv[0]*vv[0]+vv[1]*vv[1]));
+ double vv1=Math.sqrt(vv[0]*vv[0]+vv[1]*vv[1]);
+ aida.cloud1D("Evolution vradial ").fill(vv1);
+ atZmax = (Math.abs(rpn[2])>(a_zmax+(drLength/10)) )? 1 :0;
+ atR = ((Math.sqrt(rpn[0]*rpn[0]+rpn[1]*rpn[1])-(a_r))<(drLength))?1:0; // add 0.001 in mm now-Sept05
+ if(atZmax==1 || stopTkElow) return;
+ // access values of dE/dx & number of steps from stpcd
+ // Call stepBy to do the stepping && update array[] rp
+ if(partR(rpn)<2.)dT0=getDTOF();
+ int nLoop=0;
+ atR=0;
+ bigLoop:while(((atZmax!=1)&&!stopTkElow)&&((nLoop<=10)||((a_r>1000)&&(nLoop<=1270)))){
+ nLoop++;
+ if(partR(rpn)==0)break;
+ rpn = stepBy(drLength,materDeDx_i) ; // rpn each step
+ pFactor=(partPabs(rpn)>=1.)?1.:(5*partPabs(rpn));
+ setVelocity(rpn); vv=getVelocity();
+ rv = (rpn[0]*vv[0]+rpn[1]*vv[1])/(partR(rpn)*Math.sqrt(vv[0]*vv[0]+vv[1]*vv[1]));
+ setDTOF(rpn,drLength);
+ aida.cloud2D("behavior of rv in loop").fill(getDTOF(),rv);
+ pFactor=(partPabs(rpn)>=1.)?1.:(5*partPabs(rpn));
+ setCurlBack(rv<=0.35);
+ if((rv<0.35) && (Math.abs(rv)>0.35)&&(curlBackFlag==1)){
+ curlBackFlag=2; break;
+ }
+ if(( (curlBackFlag==0) && (rv<=0.35) )||(getDTOF()>1.5*dT0))
+ {
+ if(rv<=0.35)curlBackFlag=1 ;
+ drLength =(dThick>1200.)?6.:((dThick/100)*pFactor);
+ numSteps_i =(dThick>1200.)?210:(int)(dThick/drLength);
+ }
+ // Check boundaries after each step
+ atZmax = (Math.abs(rpn[2])>(a_zmax+(drLength/10)))? 1 :0; //each step
+ atR = ( Math.abs(partR(rpn)-a_r) < drLength) ? 1 : 0; //each step
+ if((atZmax==1)||stopTkElow||(atR==1)) break;
+ }
+ if(atZmax ==1 || stopTkElow ||curlBack)
+ {
+ if((atZmax==1)||stopTkElow){
+ 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");
+ }else if(curlBack){
+ System.out.println(" The track curls back go to a smaller radius");
+ }
+ }
+ double pabsExit = partPabs(rpn);
+ double dpFromDedxLay =sumdpFM;
+ rv=((rpn[0]*vv[0]+rpn[1]*vv[1]))/(partR(rpn)*Math.sqrt(vv[0]*vv[0]+vv[1]*vv[1]));
+ aida.cloud1D("Evolution of scalar productrv ").fill(rv);
+ setDTOF(rpn,drLength);
+
}
/**
* make a step along this track