lcsim/src/org/lcsim/recon/muon
diff -u -r1.8 -r1.9
--- BarrelCalSegmentFinder.java 6 Dec 2005 16:11:21 -0000 1.8
+++ BarrelCalSegmentFinder.java 7 Dec 2005 17:27:14 -0000 1.9
@@ -13,7 +13,7 @@
import org.lcsim.geometry.segmentation.ProjectiveCylinder;
import org.lcsim.geometry.util.IDDescriptor;
import org.lcsim.geometry.layer.LayerSlice;
-import org.lcsim.digisim.CalHitMapMgr;
+import org.lcsim.recon.cluster.util.CalHitMapMgr;
import org.lcsim.recon.cluster.density.MyTools;
// Muon Segment Finder for barrel calorimeter
@@ -25,7 +25,7 @@
{
public BarrelCalSegmentFinder(Detector aDet, String subdetectorName)
{
-
+
if(dedx==null) dedx = DeDx.instance();
subdetName = subdetectorName;
@@ -33,8 +33,8 @@
det = aDet;
CylindricalBarrelCalorimeter calsub = (CylindricalBarrelCalorimeter)det.getSubdetectors().get(subdetName);
- /**
- * Get materials (of 1st layer only) from LayeredCalorimeter class.
+ /**
+ * Get materials (of 1st layer only) from LayeredCalorimeter class.
*
* WARNING: This approach is prone to error, because the slices are
* not guaranteed to be the same thickness in every layer.
@@ -63,7 +63,7 @@
dataMgr = CalHitMapMgr.getInstance();
segm = calsub.getReadout().getSegmentation();
-
+
if( segm instanceof ProjectiveCylinder ) {
ProjectiveCylinder projCal = (ProjectiveCylinder)segm;
nPhi = projCal.getPhiBins();
@@ -116,7 +116,7 @@
double dr = (rmax - rmin)/nLayers; // layer thickness
//double [] rint = new double[3]; // intersection coords.
numLayersHit = 0;
-
+
double xNumSteps;
double r;
// C.M. Dec-03- Tkparams to stepper
@@ -138,7 +138,7 @@
System.out.println("Tracker r=" +r+" zmax="+zmax);
System.out.println("Barrel:out of Tracker");
}
-
+
if(subdetName.equals(hcalSubdetName))
{
System.out.println(" !!!=== rmin="+rmin+" r="+r);
@@ -151,9 +151,9 @@
double d_EMHD = Math.abs(r-rnow);
System.out.println(" distance between calorimeters EM-HD "+d_EMHD);
stepConditions = steprConditions((10*d_EMHD),BField,0.0);
- stpr.tkSteps(r,zmax,stepConditions);
+ stpr.tkSteps(r,zmax,stepConditions);
}
-
+
rpVect = stpr.getNewRp();
xNumSteps = 10.;
stepConditions = steprConditions(xNumSteps,BField,meanDEdxPerLay(subdetName));
@@ -161,27 +161,27 @@
// Forms a muon segment by looking for calorimeter hits that match
// cells which the track passes through
// Begin by finding intersection of track with innermost layer
-
+
r = rmin+0.5*dr;
-
+
stpr.tkSteps(r, zmax,stepConditions);
rpVect = stpr.getNewRp();
System.out.println(subdetName+" !!!!!! rmin="+rmin);
// Get intersections with successive layers as long as track is
// within the detector, and create a list of the cells that
// the track passes through
-
+
while ((Math.abs(rpVect[2]) < zmax )&& (nr < nLayers))
{
// Add cell to list
-
+
int iphi = xyzToPhiBin(rpVect);
int itheta = xyzToThetaBin(rpVect);
double xx = rpVect[0];
double yy = rpVect[1];
double rho = Math.sqrt(xx*xx+yy*yy);
// System.out.println("matchHits: ("+rpVect[0]+"; "+rpVect[1]+"; "+rpVect[2]+"), rho="+rho);
-
+
// System.out.println(" CellID: layer="+nr+", iphi="+iphi+", itheta="+itheta);
encoder.setValue("layer", nr);
encoder.setValue("theta", itheta);
@@ -197,7 +197,7 @@
} // Go on to the next layer
stpr.tkSteps(r, zmax,stepConditions);
rpVect = stpr.getNewRp();
- nr++;
+ nr++;
}
// Save number of layers track passed through
numLayersHit = nr;
@@ -205,8 +205,8 @@
double rr=stpr.partR();
if((Math.abs(rpVect[2])<=zmax-0.1)|| (rr>=rmax+0.1))
stpr.resetAtZmax(); // reset for next sub-detector
-
-
+
+
// Now match the cells on this list to calorimeter hits
// System.out.println("Names: subdetName=<"+subdetName+">"
// +", EcalSubdetName="+ecalSubdetName
@@ -225,56 +225,56 @@
matchHitsFast(segm, dataMgr.getCollHitMap(ecalHitmapName) );
// HistogramFolder.setDefaultFolder(".."); // C.M.-27Jan03- Folder
}
-
+
}
-
+
public void reset()
{
- // Erase lists of hits and cells
+ // Erase lists of hits and cells
hitList.removeAllElements();
cellList.removeAllElements();
}
-
+
public void setPhiNNCut(int cut)
{
-
+
// Set the number of nearest neighbors cut in phi. Cells
// with abs(phi_track - phi_cell) <= cut will be added
// to the list
-
+
phiNNCut = cut;
}
-
+
public void setThetaNNCut(int cut)
{
-
+
// Set the number of nearest neighbors cut in theta. Cells
// with abs(theta_track - theta_cell) <= cut will be added
// to the list
-
+
thetaNNCut = cut;
}
-
+
public double[] getNNCuts()
{
// Return the number of nearest neighbors cuts in phi and theta
-
+
return(new double[]
{phiNNCut, thetaNNCut});
}
-
+
public int getNLayersHit()
{
// Return the number of layers the track passed through
-
+
return(numLayersHit);
}
-
-
+
+
public int xyzToPhiBin(double[] r)
{
// Determine phi bin corresponding to this Cartesian coordinate
-
+
double phi;
int bin;
nc++;
@@ -284,7 +284,7 @@
bin = (int)(phi/phiBinSize);
return(bin);
}
-
+
public int xyzToThetaBin(double[] r)
{
// Determine theta bin corresponding to this Cartesian coordinate
@@ -297,15 +297,15 @@
bin = (int)(theta/thetaBinSize);
return(bin);
}
-
-
+
+
// The following two methods implement the CalorimeterHits interface
-
+
public Enumeration getHits()
{
return(hitList.elements());
}
-
+
public int getNHits()
{
return(hitList.size());
@@ -340,9 +340,9 @@
return matterDEdx ;
}
-
+
//------------------------------------------------------
-
+
protected void matchHitsFast(Segmentation cell,
Map<Long,SimCalorimeterHit> calHits)
{
lcsim/src/org/lcsim/util/step
diff -u -r1.5 -r1.6
--- TrackStepper.java 6 Dec 2005 16:33:54 -0000 1.5
+++ TrackStepper.java 7 Dec 2005 17:27:15 -0000 1.6
@@ -3,7 +3,7 @@
import java.util.*;
import java.text.*;
import org.lcsim.util.aida.*;
-import org.lcsim.digisim.CalHitMapMgr;
+import org.lcsim.recon.cluster.util.CalHitMapMgr;
/**Given a particle of initial position momentum and charge
* (x,y,z,px,py,pz,and q),calculates its trajectory in a uniform
@@ -17,7 +17,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.5 2005/12/06 16:33:54 caroline Exp $
+ * @version $Id: TrackStepper.java,v 1.6 2005/12/07 17:27:15 lima Exp $
*/
public class TrackStepper
{
@@ -40,12 +40,12 @@
private double [] vxyz = new double[3];
private double dTof;
private double dL ;
-
+
double []xy = new double[3];
private double []dpF = new double[3];
private double []dpM = new double[3];
private double []vstep = new double[3];
-
+
// Status Flags
private boolean stopTkElow = false;
private int stepFlags = 0 ;
@@ -53,18 +53,18 @@
private int atR = 0 ;
// Material Information
private int numSteps_i;
- protected double BField_i =-10. ;
+ protected double BField_i =-10. ;
private double materDeDx_i ;
-
+
// private Detector detNow ;
// Multiple Local use
private double tmpxyz ;
private double tmpp ;
private double mass_i ;
private double particleR;
- private double particlePt;
+ private double particlePt;
private double particlePabs;
-
+
/** Default Constructor */
public TrackStepper()
{
@@ -93,33 +93,33 @@
// 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
-
+
rpn=trakPs.getRp(); // get phaseSpace at entry
-
+
double a_r = r;
double a_zmax=zmax;
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
numSteps_i =(int) stpcd[0]; BField_i = stpcd[1];
double d=Math.abs(Math.sqrt(rpn[0]*rpn[0]+rpn[1]*rpn[1])-a_r) ;
// Check if distance cm or distance mm !!!!
materDeDx_i= stpcd[2]/d ;
-
+
// Call stepBy to do the stepping && update array[] rp
- double a_dc = d/numSteps_i ;
+ double a_dc = d/numSteps_i ;
System.out.println("Before "+numSteps_i+" steps particle Radius="+partR()+" Pabs="+ partPabs()); ;
-
+
for( int i=0; (i<numSteps_i);i++)
{
rpn = stepBy(a_dc,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
@@ -147,7 +147,7 @@
dL = getPathLength();
dpF = dpFromField(dc);
dpM = dpFromMatter(dc,matDeDx);
-
+
//System.out.println(" !!!=> px,dpF0,dpM0="+rpn[3]+" "+dpF[0]+" "+dpM[0]);
//System.out.println("!!!=> py,dpF1,dpM1="+rpn[4]+" "+dpF[1]+" "+dpM[1]);
vstep = getVelocity();
@@ -168,26 +168,26 @@
// Histogram.find("py vs px").fill(rpn[3],rpn[4]);
return trakPs.upDateRp(rpn);
}
-
+
/** Return dpx dpy due to the Field Bz during this step
* has been specifically adapted to the unique LCD tracking coordinate
* system. (Inverted sign of curvatures.) by changing the sign of the
* magnetic field to -Bz
*/
-
+
private double[] dpFromField(double dc)
{
- double alpha; double delta; double adjust; double mainTerm;
+ double alpha; double delta; double adjust; double mainTerm;
double pabs2 = rpn[3]*rpn[3]+rpn[4]*rpn[4]+rpn[5]*rpn[5];
mass_i = rpn[6];
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;
// assert BField_i > 4.0 : "BField too small!";
-
+
alpha = a_q*0.3*(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;
@@ -203,23 +203,23 @@
*/
private double[] dpFromMatter(double dc,double materMeanDeDx)
{
-
+
double pabs2= rpn[3]*rpn[3]+rpn[4]*rpn[4]+rpn[5]*rpn[5];
materDeDx_i = materMeanDeDx;
double tempdpoverdl;
double mass_i = rpn[6];
double a_En = Math.sqrt(pabs2 + mass_i*mass_i);
-
+
for (int i =0; i<3;i++)
{
// Change in Momentum due to material
tempdpoverdl = (materDeDx_i*a_En*(rpn[i+3]/pabs2));
dpFordL[i] = tempdpoverdl*dL;
}
-
+
return (dpFordL);
}
-
+
/**
* Return vector rp(x,y,z,px,py,pz,mass_i=MuonMass,charge,BField)
* from the information at the Interaction Point given by TrkParams.
@@ -236,7 +236,7 @@
trakPs.trkRpFromParams(tpar,inM, BField_i);
return trakPs.getRp();
}
-
+
/**
* Return vector rp(x,y,z,px,py,pz,mass_i=inMass,charge,BField)
* from the information at the Interaction Point given by TrkParams.
@@ -253,7 +253,7 @@
trakPs.trkRpFromParams(tpar,inMass, BField_i);
return trakPs.getRp();
}
-
+
/** set velocity vx,vy,vz of the particle this step
*/
public void setVelocity()
@@ -267,7 +267,7 @@
tempV = (rpn[i+3]/En)*clight*1000. ;
vxyz[i]=tempV ;
}
- }
+ }
/**
* Return velocity vx,vy,vz of the particle this step
*/
@@ -275,14 +275,14 @@
{
return vxyz;
}
-
+
/**
* Set dTimeOfFlight for the step dc
* @param dc step
*/
public void setDTOF(double dc)
{
-
+
double a=0.;double b=0.; double c=0.;
double a_dc = dc ;
xy[0]=rpn[0]; xy[1]=rpn[1];xy[2]=rpn[2];
@@ -312,7 +312,7 @@
double vabs = Math.sqrt(vstep[0]*vstep[0]+vstep[1]*vstep[1]);
dL = dTof * vabs;
}
-
+
/** @return change in pathLength this step
*/
public double getPathLength()
@@ -323,35 +323,35 @@
*/
public double[] getRAfterStep()
{
-
+
return r_ints =trakPs.getPosition();
}
/** @return Momentum (px,py,pz) after step
*/
public double[] getPAfterStep()
- {
+ {
return p_ints=trakPs.getMomentum();
}
/** @return phase-space point after step
*/
public double [] getNewRp()
{
-
+
return rpn= trakPs.getRp();
}
- /** @return Radius at which the particle is
+ /** @return Radius at which the particle is
*/
public double partR()
{
particleR= rpn[0]*rpn[0]+rpn[1]*rpn[1];
if(particleR!=0.)particleR = Math.sqrt(particleR);
return particleR;
- }
+ }
/** @return The Pt of the particle
*/
public double partPt()
{
-
+
particlePt= rpn[3]*rpn[3]+rpn[4]*rpn[4];
if(particlePt!=0.) particlePt=Math.sqrt(particlePt);
return particlePt;
@@ -360,7 +360,7 @@
*/
public double partPabs()
{
-
+
particlePabs= rpn[3]*rpn[3]+rpn[4]*rpn[4]+rpn[5]*rpn[5];
if(particlePabs!=0.) particlePabs= Math.sqrt(particlePabs);
return particlePabs;
@@ -372,10 +372,10 @@
{
resetStopTkELow();
resetAtZmax() ;
-
+
int stepFlags = 0;
}
-
+
/** Return the Stepping flags Indicator of curlBack Left for the step
*/
public int getStepFlags()
@@ -425,8 +425,8 @@
*/
public int getAtR()
{return atR;}
-
-
+
+
/** Reset- clear by calling clearTrackStepper()
*/
public void clear()
@@ -438,7 +438,7 @@
public void clearTrackStepper()
{
trakPs.clear();
- BField_i = -10.;
+ BField_i = -10.;
setDTOF(0.);
setPathLength(0.);
dpFromField(0.);
@@ -453,6 +453,6 @@
xy[i]=0.;
}
}
-
+
}