lcsim/src/org/lcsim/util/step
diff -u -r1.14 -r1.15
--- TrackStepper.java 14 Jul 2006 19:14:55 -0000 1.14
+++ TrackStepper.java 10 Oct 2006 02:01:45 -0000 1.15
@@ -19,14 +19,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.14 2006/07/14 19:14:55 caroline Exp $
+ * @version $Id: TrackStepper.java,v 1.15 2006/10/10 02:01:45 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
+ 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
{
@@ -85,8 +85,8 @@
private double sumdL ;
private double sumdpFM ;
public StepConditions cond, condMat;
-
-
+
+
/** Default Constructor */
public TrackStepper()
{
@@ -108,7 +108,7 @@
* @param stpcd
* an array ={xNumberOfSteps, BField, Mean_dE/dx_this layer}
*/
-
+
/**
Gets coilSlices from class StepConditions
*/
@@ -117,93 +117,7 @@
cond = new StepConditions(det);
return cond.getCoilSlices();
}
- //default
- 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;
- dThick = stpcd[0];
- BField_i = stpcd[1];
-
- 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
- // Check if distance cm or distance mm !!!!
-
- 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);
-
- }
//for the coil
public void tkSteps(double r,double zmax,Detector det, int sliceN)
{
@@ -215,40 +129,137 @@
double []b = new double[3];
b = cond.getField(rpn);
BField_i = stpcd[1] =b[2] ;
+ System.out.println(" sliceN Field="+b[2]);
stpcd[2]= cond.getSliceDeDx(sliceN);
Stepping(r, zmax, stpcd);
}
//for the material
- public void tkSteps(double r,double zmax,Detector det, String materialname)
+ public void tkSteps(double ToGoPos,double limitingVar,String BarrEndc,double magfild, String materialname)
{
-
+ System.out.print("stpr:I am Here1");
double[]stpcd = new double[3];
rpn=trakPs.getRp(); // get phaseSpace at entry
- condMat = new StepConditions(rpn,det, materialname, r);
-
+ System.out.println("rpn-x="+rpn[0]+" y="+rpn[1]+" z="+rpn[2]);
+ condMat = new StepConditions(rpn,BarrEndc, materialname, ToGoPos);
+ System.out.println("after stpcd:rpn-x="+rpn[0]+" y="+rpn[1]+" z="+rpn[2]);
+
dThick = stpcd[0] = condMat.getMaterialThickness();
+ System.out.println("dThick="+condMat.getMaterialThickness());
+ System.out.println("tkSteps:HELLO");
double []b = new double[3];
- b = condMat.getField(rpn);
- BField_i = stpcd[1] =b[2] ;
+ //b = condMat.getField(rpn);
+ System.out.println("Material Field="+b[2]);
+ //BField_i = stpcd[1] =b[2] ;
+ BField_i = stpcd[1]=magfild;
stpcd[2]= condMat.getMatterDeDx();
- Stepping(r, zmax, stpcd);
-
+ System.out.println("dThick="+dThick+" Field="+b[2]+"materialdEdx="+ stpcd[2]);
+ Stepping(ToGoPos, limitingVar, stpcd);
+
}
//for the subdetector
- public void tkSteps(double r,double zmax, Detector det,String subdetName, int layerN)
+ public void tkSteps(double ToGoPos,double limitingVar, Detector det,String subdetName, int layerN)
{
double[]stpcd = new double[3];
rpn=trakPs.getRp(); // get phaseSpace at entry
- cond = new StepConditions(r, rpn, det, subdetName, layerN);
- dThick = stpcd[0] = cond.getLayerThickness();
+ cond = new StepConditions(ToGoPos, rpn, det, subdetName, layerN);
+ dThick = stpcd[0] = cond.getLayerThickness();
double []b = new double[3];
b = cond.getField(rpn);
BField_i = stpcd[1] =b[2] ;
+ System.out.println(" LayerN Field="+b[2]);
stpcd[2] = cond.getlayerDeDx();
- Stepping(r, zmax, stpcd);
+ Stepping(ToGoPos, limitingVar, stpcd);
}
-
+
+ //default
+ private void Stepping(double r,double zmax,double[]stpcd)
+ {
+
+ if(r>=1270)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
+ if(r>=1270)System.out.print("rpn: x="+rpn[0]+" y="+rpn[1]+" z="+rpn[2] );
+ setVelocity(rpn) ; // particle velocity
+ double pabsEntry =partPabs(rpn) ;
+ a_r = r;
+ a_zmax=zmax;
+ dThick = stpcd[0];
+ BField_i = stpcd[1];
+
+ 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
+ // Check if distance cm or distance mm !!!!
+
+ 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
* @return position after step
