lcsim/src/org/lcsim/util/step
diff -u -r1.13 -r1.14
--- TrackStepper.java 14 Jul 2006 15:54:21 -0000 1.13
+++ TrackStepper.java 14 Jul 2006 19:14:55 -0000 1.14
@@ -19,7 +19,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.13 2006/07/14 15:54:21 caroline Exp $
+ * @version $Id: TrackStepper.java,v 1.14 2006/07/14 19:14:55 caroline Exp $
*NEXT To BE DONE: get subdet/layer info diretly
*get materiel info from a util/step/stepperCondition.class
*
@@ -110,7 +110,7 @@
*/
/**
- Gets coilSlices from class SteConditions
+ Gets coilSlices from class StepConditions
*/
public Vector<CoilSubLayer> getCoilSlices(Detector det)
{
@@ -118,7 +118,7 @@
return cond.getCoilSlices();
}
//default
- public void tkSteps(double r,double zmax,double[]stpcd)
+ private void Stepping(double r,double zmax,double[]stpcd)
{
if(debug) System.out.println("Stepper.tkSteps() called: rstep: r="+r+", zmax="+zmax);
@@ -133,7 +133,8 @@
double pabsEntry =partPabs(rpn) ;
a_r = r;
a_zmax=zmax;
- dThick = stpcd[0]; BField_i = stpcd[1];
+ dThick = stpcd[0];
+ BField_i = stpcd[1];
sumdL =0. ;
sumdpFM=0. ;
@@ -208,165 +209,29 @@
{
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(det);
-
- double pabsEntry =partPabs(rpn) ;
- a_r = r;
- a_zmax=zmax;
-
+ cond = new StepConditions(det);
dThick = stpcd[0] = cond.getSliceThickness(sliceN);
double []b = new double[3];
b = cond.getField(rpn);
BField_i = stpcd[1] =b[2] ;
stpcd[2]= cond.getSliceDeDx(sliceN);
- 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);
-
+ Stepping(r, zmax, stpcd);
}
//for the material
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);
+ 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(rpn);
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);
+ Stepping(r, zmax, stpcd);
}
//for the subdetector
@@ -374,159 +239,16 @@
{
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();
+ rpn=trakPs.getRp(); // get phaseSpace at entry
+ cond = new StepConditions(r, rpn, det, subdetName, layerN);
+ dThick = stpcd[0] = cond.getLayerThickness();
double []b = new double[3];
b = cond.getField(rpn);
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]));
- 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);
-
- }
- 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);
-
+ Stepping(r, zmax, stpcd);
}
+
/**
* make a step along this track
* @return position after step