lcsim/src/org/lcsim/util/step
diff -u -r1.12 -r1.13
--- TrackStepper.java 13 Jul 2006 23:35:11 -0000 1.12
+++ TrackStepper.java 14 Jul 2006 15:54:21 -0000 1.13
@@ -3,6 +3,7 @@
import java.util.*;
import java.text.*;
import org.lcsim.util.aida.*;
+import org.lcsim.recon.muon.CoilSubLayer;
import org.lcsim.recon.cluster.util.CalHitMapMgr;
import org.lcsim.geometry.Detector;
@@ -18,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.12 2006/07/13 23:35:11 caroline Exp $
+ * @version $Id: TrackStepper.java,v 1.13 2006/07/14 15:54:21 caroline Exp $
*NEXT To BE DONE: get subdet/layer info diretly
*get materiel info from a util/step/stepperCondition.class
*
@@ -85,6 +86,7 @@
private double sumdpFM ;
public StepConditions cond, condMat;
+
/** Default Constructor */
public TrackStepper()
{
@@ -107,8 +109,16 @@
* an array ={xNumberOfSteps, BField, Mean_dE/dx_this layer}
*/
-
- public void tkSteps(double r,double zmax,double[]stpcd)
+ /**
+ Gets coilSlices from class SteConditions
+ */
+ public Vector<CoilSubLayer> getCoilSlices(Detector det)
+ {
+ cond = new StepConditions(det);
+ return cond.getCoilSlices();
+ }
+ //default
+ public void tkSteps(double r,double zmax,double[]stpcd)
{
if(debug) System.out.println("Stepper.tkSteps() called: rstep: r="+r+", zmax="+zmax);
@@ -193,8 +203,90 @@
setDTOF(rpn,drLength);
}
-
-
+ //for the coil
+ public void tkSteps(double r,double zmax,Detector det, int sliceN)
+ {
+
+ 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;
+
+ 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);
+
+ }
+ //for the material
public void tkSteps(double r,double zmax,Detector det, String materialname)
{
@@ -212,7 +304,7 @@
dThick = stpcd[0] = condMat.getMaterialThickness();
double []b = new double[3];
- b = condMat.getField();
+ b = condMat.getField(rpn);
BField_i = stpcd[1] =b[2] ;
stpcd[2]= condMat.getMatterDeDx();
sumdL =0. ;
@@ -277,7 +369,7 @@
setDTOF(rpn,drLength);
}
-
+ //for the subdetector
public void tkSteps(double r,double zmax, Detector det,String subdetName, int layerN)
{
@@ -294,7 +386,7 @@
a_zmax=zmax;
dThick = stpcd[0] = cond.getLayerThickness();
double []b = new double[3];
- b = cond.getField();
+ b = cond.getField(rpn);
BField_i = stpcd[1] =b[2] ;
stpcd[2] = cond.getlayerDeDx();